{'fraction_non_alphanumeric': 0.029695144096220812, 'fraction_numerical': 0.010596252857356492, 'mean_word_length': 4.801933564700049, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Evaluating Approximate Flavor Instability Metrics in Neutron Star Mergers

Neutrinos can rapidly change flavor in the inner dense regions of core-collapse supernovae and neutron star mergers due to the neutrino fast flavor instability. If the amount of flavor transformation is significant, the FFI could significantly affect how supernovae explode and how supernovae and mergers enrich the universe with heavy elements. Since many state of the art supernova and merger simulations rely on neutrino transport algorithms based on angular moments of the radiation field, there is incomplete information with which to determine if the distributions are unstable to the FFI. In this work we test the performance of several proposed moment-based instability tests in the literature. We perform time-independent general relativistic neutrino transport on a snapshot of a 3D neutron star merger simulation to generate reasonable neutrino distributions and check where each of these criteria correctly predict instability. In addition, we offer a new ``maximum entropy'' instability test that is somewhat more complex, but offers more detailed (though still approximate) estimates of ELN crossing width and depth. We find that this maximum entropy test and the resonant trajectory test are particularly accurate at predicting instability in this snapshot, though all tests predict instability where significant flavor transformation is most likely.


I. INTRODUCTION
Neutrinos come in one of three known flavors associated with the three known charged leptons. In the canonical theory for core-collapse supernovae (CCSNe), neutrinos are the dominant means by which energy is transported outward, enabling the collapse of the stellar core to result in the explosion of the rest of the star [1][2][3]. In neutron star mergers (NSMs), the neutrino losses determine the thermal evolution of the disk [4,5], and in both cases neutrino irradiation determines the nuclear composition of the ejecta that enriches the universe with heavy elements (e.g., [6,7]). Electron neutrinos play a unique role because they participate in charged-current reactions that transform neutrons into protons and vise versa. Because of this, models of the dynamics and nucleosynthesis in these systems are sensitive to the generation, movement, and absorption of neutrinos of different flavors [4,8].
The propensity of neutrinos to change flavor in flight thus poses a significant challenge to these models (see [9,10] for reviews). In CCSNe and NSMs, neutrinos can be sufficiently dense that neutrino-neutrino interactions significantly modify flavor transformations in a nonlinear manner that has yielded a rich phenomenology, including collective oscillations [9], the neutrino halo effect [11], and the matter-neutrino resonance [12]. More recently, the neutrino fast flavor instability (FFI) [13,14] was shown to occur nearly ubiquitously in CCSNe (e.g., [15][16][17][18][19]) and NSMs [20][21][22][23] in a way that may significantly modify the nuclear composition of ejected matter. The FFI can occur in regions inaccessible to other flavor transformation mechanisms, but the short timescales and lengthscales associated with the FFI preclude a direct treatment in * srichers@berkeley.edu global simulations. Because of this, insight is needed to predict where the instability occurs and the net effect it produces.
Although global simulations of flavor transformation in CCSNe and NSMs are not yet possible, local dynamical simulations of the FFI (e.g., [24][25][26][27][28][29][30][31][32][33][34][35][36][37]) and analytic calculations [33,[38][39][40] are able to predict the post-instability equilibrium with an increasing realism, but a general solution is still lacking. Fortunately, linear stability analysis can be used to predict where in a CCSN or NSM the FFI occurs, even if it cannot predict the nonlinear behavior of the instability after the instability saturates. Following an extensive history in application to other collective neutrino instabilities (e.g., [14,17,[41][42][43][44][45][46][47]) linear stability analysis has led to a straightforward criterion for instability: a neutrino distribution is unstable to the FFI at any location where there is a propagation direction along which there is an equal number of neutrinos and antineutrinos [48][49][50]. This simple concept is very amenable to post-processing of global simulations that do not include flavor transformation (e.g., [16][17][18][19][51][52][53][54]). However, in many of these calculations, the full neutrino distribution is not calculated, as only angular moments of the neutrino field are simulated [55,56] to reduce the computational cost. A number of methods have been proposed that use the limited information present in these moments to predict whether an ELN crossing, and thus the FFI, is present for a particular combination of angular moments. These include the k 0 test [46], the polynomial test [52], the α = 1 test [57,58], the unstable pendulum and resonant trajectory tests [26], and a fitting method to more exact calculations [54]. In this work, we also provide analytic expressions for a general maximum entropy test, another ELN crossing test based on the shape of the distribution used in the popular maximum entropy closure [26,59,60].
The k 0 and polynomial tests have been dynamically incorporated into global simulations of NSMs [22,23], and similar simulations of CCSNe are likely underway (though [61] incorporate flavor transformation into models of neutrino-driven wind from a protoneutron star). The ability of some of the tests to accurately predict instability has been tested in the context of simulations of one-dimensional (spherically symmetric) CCSNe with Newtonian gravity [62][63][64], in which hand-tuned criteria can perform quite well, but more conservative criteria naturally tend to under-predict instability. It is also not clear how these results extend to . In this work, we take each of these tests and assess how well they perform in the context of a three-dimensional NSM with a general spacetime metric. With this information, we hope to provide some insight into the biases associated with each test such that flavor transformation can be more realistically incorporated into simulations of NSMs and so we can better interpret the results of such simulations.
We begin in Section II by reviewing the definition of ELN crossings, taking care to discuss classical neutrino distributions without reference to quantum kinetics. In Section III we describe the time-independent Monte Carlo radiation transport method we use to calculate the full neutrino distribution information and review how analytic closures are used to determine higher angular moments of the radiation field from the energy density and flux. In Section IV, we describe the structure of the resulting full radiation field, derive our new maximum entropy crossing test, and demonstrate the ability of each of the proposed tests to accurately detect crossings in the ELN distributions. Finally, we provide some concluding remarks in Section V.

II. ELN CROSSINGS AS AN INDICATION OF THE FAST FLAVOR INSTABILITY
In this section, we briefly review the conditions for the growth of the neutrino fast flavor instability. For the sake of simplicity, we make no reference to the quantum kinetic equations in this work, and instead appeal to the equivalence between instability and crossings in the angular distribution of electron lepton number [49]. Understanding the origin of the instability and how it evolves requires a treatment of the quantum kinetic equations, but identifying instability in a distribution of neutrinos in pure flavor states requires only knowledge of each flavor's distribution. We also assume a flat spacetime in this discussion, since the fast flavor instability tends to operate on length scales much smaller than the spacetime curvature.
The distribution of each neutrino species ν a is represented by the distribution function f νa (x, Ω, , t), which for neutrinos takes on values of 0 ≤ f νa ≤ 1. The distribution function is a seven-dimensional function of the position x, direction unit vector Ω, the energy , and the time t. The number density, number flux, number "pressure tensor", and number "heat tensor" of each species are integrals of the distribution function over momentum: The neutrino lepton number distribution is defined as an energy integral of the difference between the distributions of a pair of neutrino flavors. Specifically, (2) There is a neutrino lepton number crossing, and thus flavor instability [49], at any x and t where G νaν b takes on positive values in some directions and negative values in others. It is common to assume that due to the energy scales involved in core-collapse supernovae and neutron star mergers, interactions producing heavy leptons are kinematically suppressed, and the distributions of mu and tau neutrinos and antineutrinos are all similar. Although small deviations from this assumption can be important (see [62,65]), we do assume that heavy lepton neutrinos have the same distribution for the sake of analyzing crossings in the electron flavor sector. With this assumption, the neutrino lepton number distributions become In the following, we refer simply to the electron lepton number (ELN) distribution G to mean either G νeνµ or G νeντ . Angular moments of the lepton number distribution are defined analogously to moments of the distribution function itself. That is, Note that the moment subscript also denotes the tensor rank of the moment. I 0 is a scalar, I 1 is a vector, I 2 is a rank-2 tensor, etc.
Many of the instability metrics described in this work were derived assuming the distributions of all species exhibit axial symmetry around the same axis, which is generally not true in three-dimensional systems, and especially difficult to justify in neutron star mergers [60,66]. In order to connect with work that assumes axial symmetry, we extract the component of each moment along the direction of the net ELN fluxÎ 1 . The scalarized moments are then We use Einstein summation notation in these expressions, but since everything is defined in an orthonormal tetrad, down-index quantities are identical to up-index quantities.
The rest of this work is devoted to assessing how well we can predict the presence of the fast flavor instability using only these angular moments of the ELN.

III. METHODS
We perform general-relativistic Monte Carlo neutrino radiation transport to directly solve for a realistic distribution function everywhere on the domain and determine where there are ELN crossings. We then take angular moments of the distribution in order to assess how well several moment-based tests are able to detect ELN crossings. In this section, we lay out the details of the Monte Carlo calculation that yield the full angular distribution. We then review how an analytic closure is used to estimate the pressure and heat tensors when only the number density and number flux are known.

A. Monte Carlo Radiation Transport
We use SedonuGR [60] to calculate the steady-state radiation field in a snapshot of a three-dimensional neutron star merger simulation from [66] at 10 ms after merger. SedonuGR imports the mass density ρ, electron fraction Y e , temperature T , and spacetime metric g µν at every point in space. We use only one refinement level spanning a domain of size 563 km × 563 km × 145 km (assuming reflection symmetry across z = 0) with a grid size of 207 × 207 × 54, corresponding to a spatial resolution of 2.7 km in all directions.
The emissivity η, absorption opacity κ abs , and elastic scattering opacity κ scat for each neutrino species is determined by NuLib [67], including charged current absorption/emission on nucleons and nuclei, elastic scattering on nucleons, nuclei, and electrons, and neutrino pair creation and annihilation (including nucleon-nucleon Bremsstrahlung). For pair processes, the neutrino annihilation rate is determined by applying Kirchoff's law to the emissivity. We use the LS220 equation of state [68] to calculate the neutrino interaction rates, consistent with that used in the simulation that produced the background data. The steady-state approximation and the approximate treatment of the scattering and pair processes are not realistic, but they suffice to produce believable distributions of neutrino radiation that we can use to judge schemes for detecting an ELN crossing. We employ 12 energy groups with upper bounds logarithmically spaced from 4 to 150 MeV.
We outline the major features of the Monte Carlo method, but refer the reader to [60] for details. Sedonu initializes a large number of Monte Carlo particles in each grid cell and each energy bin. In this work, we create a total of 1.2 × 10 10 Monte Carlo particles for each flavor. Each particle is given a weight N (i.e., the number of physical neutinos the particle represents) according to the emissivity of each neutrino species in that spaceenergy zone. The direction of each particle is isotropically randomly sampled in the frame comoving with the fluid. The distance (again in the comoving frame) to the next scattering event is randomly sampled from an exponential distribution. The particle then propagates that distance or to the next grid cell wall (whichever is closer) according to the geodesic equation, and the comovingframe distance traversed is labeled ∆s. Throughout this step, the weight of each particle is continuously decreased according to the absorption opacity. If the distance chosen was the scattering distance, the particle is then given a new random direction in the new comoving frame, preserving the neutrino energy in that frame. In any case, all opacities and metric quantities are then re-interpolated from the background grid, a new distance is sampled, and the process repeats until the particle weight decreases below a threshold (in which case it is rouletted) or it leaves the domain of the calculation.
Each space-energy zone collects radiation information from the particles that pass through it. This aggregate radiation field is discretized into discrete direction bins, with 16 bins uniformily spaced in azimuthal angle around theẑ axis, and 8 polar bins uniformily spaced in the cosine of the angle from the same axis, all defined in a comoving orthornormal tetrad. During the step, each particle contributes a bit of energy density ∆E to the radiation field stored in the space-energy-direction zone it occupies. The energy density contribution is determined by where N is the average neutrino weight during the step (recall, it is changing due to absorption) and V is the fourvolume of the grid cell. Both the emission (which determines the initial weight N ) and the four-volume assume a particular coordinate time interval ∆t, but this arbi-trary choice cancels in the energy density accumulation. By the end of the calculation, each zone contains contributions from many separate particles. In addition, to remove Monte Carlo noise, we perform gaussian smoothing in space (x, y, z) with a width of 1 grid cell and a maximum extent of 1 grid cell. Doing the data analysis with and without this smoothing allows us to confirm that our results do not vary under differing amounts of noise.
The main output of the Monte Carlo transport is sixdimensional energy density grid E lmnpqr , where (l, m, n) are spatial grid cell indices, p is the energy bin index, q is the azimuthal angle grid index, and r is the polar angle grid index, all defined in a comoving orthonormal tetrad. We also define a unit vector pointing to the center of each direction bin Ω i qr and the bin-center neutrino energy p . The first four number density moments for each neutrino flavor can then be evaluated as straightforward sums over the energy density array:

B. Maximum Entropy Closure
We briefly review the classical maximum entropy closure of [59], as this is currently the most popular choice of analytic closures in modern moment-based neutrino transport methods. The closure also lends itself to an approximate crossing test described in Section IV B. Maximizing the angular entropy of the energy-integrated distribution constrained to a given number density and flux of each neutrino species as described in [59] yields a functional form of the distribution at each location: Although n, F, and Z are different for each species, we drop the species subscripts for the rest of this section with the understanding that this whole process is applied separately to each species. Here, θ is the angle between Ω and the direction of the net number flux of the given neutrino species, such that cos θ = Ω·F . Z is a parameter determined by solving the transcendental equation The left-hand side of this equation is the flux factor, defined as f = |F|/n.
When using any analytic closure, all components of the pressure and heat tensor are constructed by interpolating between the optically thick and thin limits as Again taking advantage of our orthonormal tetrad for simplicity, these thick and thin limits are In the expressions for L thick , repeated indices are assumed to be the same and distinct indices are assumed to be different. All components of both tensors can be determined noting that they are symmetric upon exchange of any pair of indices.
Taking the second and third angular moments of the maximum entropy distribution along the flux direction yields the familiar closure relations [59,60] In the case of spectral transport, this process is generally applied separately to each energy bin, as we do when evaluating the "closed" moments later in this work. Also, the closure is usually, though not always, applied to energy moments rather than number moments (i.e., with one more factor of energy in the integrand). However, in the case of simulated spectral moments where one assumes for numerical purposes that all of the radiation within a bin has the same energy, the result is identical.

IV. RESULTS
In this section, we assess how well various momentbased ELN crossing detection schemes perform. In order The flux directions differ most significantly in the polar regions and just above the accretion disk. Imposing a closure does not change these quantities, since they are defined with only the first two moments.
to put those results in context, we first describe some of the prevalent features of the neutrino radiation field. We then show where ELN crossings occur in the full radiation field data in Section IV A. We introduce the generalized maximum entropy crossing test in Section IV B and demonstrate the efficacy each moment-based crossing test in Sections IV B-IV F.
We show results of the Monte Carlo radiation transport calculation in Figure 1 to demonstrate differences between the electron neutrino and antineutrino distributions that lead to ELN crossings. The top panel shows the lepton number density asymmetry, where dark red implies 100% of the neutrinos are electron neutrinos, and dark blue implies that 100% of the neutrinos are electron anti-neutrinos. Far from the merger, there is an overall over-abundance of electron anti-neutrinos, a reflection of the fact that the neutron star matter is by and large increasing its electron fraction and emitting antineutrinos. There is an over-abundance of electron neutrinos in the hot and dense parts of the inner accretion disk where the electron fraction is somewhat higher (up to about 0.25), since electron antineutrinos are able to escape more easily. Finally, there is a significant over-abundance of electron anti-neutrinos in the central hypermassive neutron star because the neutron chemical potential significantly exceeds the proton and electron chemical potentials. The black contour (repeated in all other plots in this work) shows where electron neutrinos and antineutrinos have an equal number density, guaranteeing the presence of an ELN crossing [57,58]. Even in regions where one flavor is significantly more abundant than another, an ELN crossing is possible if the fluxes of the two distributions are sufficiently different. The center panel shows the difference between the electron neutrino and antineutrino flux factors. In the central hypermassive neutron star, both flux factors are approximately 0 (hence a difference also of 0). Far from the merger, both flux factors approach 1 (also trending toward a difference of 0). In intermediate equatorial regions (10 km x 200 km), the electron antineutrino flux factor is significantly larger than the electron neutrino flux factor, a result of the fact that the electron antineutrino interaction rates are smaller, allowing them to decouple from the fluid more easily.
The lower panel shows the angle between the angle between the electron neutrino and antineutrino fluxes. At large radii, all fluxes trend toward pointing radially. We color by the logarithm of the angle in order to better show small differences between the fluxes at |x| 150. Even though the electron anti/neutrino flux factors differ significantly in these regions, the flux directions differ by at most a few degrees and do not exhibit as much structure as the flux factors.
In the following sections, we demonstrate that these differences between electron neutrino and antineutrino distributions lead to ELN crossings and assess how well these crossings are detected by various moment-based tests.

A. Direct Crossing Search
When the full distribution of neutrinos is available, one can search for ELN crossings without approximation beyond the numerical discretization. However, there is generally a trade-off in simulations between the accuracy of the radiation transport and other components of the simulation, so this is only possible for a small subset of simulations. In the language of our discrete energy density array output from the Monte Carlo calculation (see Section III A), the discrete ELN distribution is The coarseness of our angular grid prevents us from detecting ELN crossings smaller than the angular grid cell size of ∆φ = 22.5 • . Small and shallow crossings seem to lead to minimal flavor transformation [35], so although finer angular resolution may reveal slightly more volume with an ELN crossing, we do not expect this to significantly influence implications for flavor transformation, mass ejection, and nucleosynthesis. In addition, we verified that we use sufficient angular resolution by checking that all results are unchanged resulting from another Monte Carlo calculation in which we collect the radiation field directly into angular moments (see [60] for details).
The regions containing an ELN crossing are shown as green in Figure 2. There are no ELN crossings within the hypermassive neutron star because the distributions of both electron neutrinos and antineutrinos is nearly isotropic and there is a strong over-abundance of electron antineutrinos. However, in the rest of the domain, the relative amounts of each species are much more similar (top panel of Figure 1), and variations in the angular distributions are sufficient to induce crossings in 98% of the domain. Extrapolating beyond the calculation domain, collisional processes are very weak due to low densities. Trajectories that have equal numbers of neutrinos and antineutrinos will remain so, implying that crossings should remain present at larger distances, in agreement with [21][22][23]. In following subsections, we will try to reproduce these results with a variety of moment-based tests.

B. Maximum Entropy Test
We provide analytic expressions for a generalized ELN crossing test based on the assumption that the energyintegrated neutrino distributions follow the form assumed in deriving the maximum entropy closure (Equation 8). This is in general not a valid assumption, since even if each of several neutrino energy bins follows a maximum entropy distribution, the sum of the distributions from those energy bins (each with different flux factors and directions) is not a maximum entropy distribution. However, we will see that it is nevertheless useful for estimating other properties of the distributions, although its applicability is limited in large-scale simulations due to the need to iteratively solve a transcendental equation. Note that Johns and Nagakura 26 use this concept The radius represents the differential number density of electron neutrinos (blue) and antineutrinos (red) propagating in the direction given by the angle θ from some arbitrary direction θ = 0. In most directions, there are more antineutrinos than neutrinos (shaded red), but the directions between the ELN crossings (black points) are dominated by electron neutrinos (shaded blue).
to analyze spherically symmetric neutrino distributions, but the approach presented here allows the neutrino and antineutrino fluxes to point in arbitrary directions.
There is an ELN crossing at any direction where G = 0. We can make intuitive sense of the crossings by taking a cross-section of the distributions in momentum space as in Figure 3, plotting the differential number density of each neutrinos (blue) and antineutrinos (red) in each direction θ as the radial coordinate of the curve. In this example, the two distributions cross at the black points, which are part of a continuous loop passing through the plane. Given the number densities and number fluxes of each distribution (and thus also Z from Equation 9), it is straightforward to determine the directions in this plane where the distributions cross by solving This can be expressed more simply as where α =Z sin θF − Z sin θ F , γ =Z cos θF , and η = ln(nZ sinhZ/nZ sinh Z). Further defining the variable θ = tan −1 (α/γ) allows a simple expression for the angles  Equation 19). The ME test correctly predicts wide crossings out to r 200 km, but falsely predicts radial structures with alternating wide and no crossings outside this region.
at which a crossing occurs A crossing exists if θ is real. Therefore, the condition for instability is Nonlinear flavor transformation simulations are still needed to precisely predict the implications for the eventual flavor transformation, but previous work has indicated that wide crossings are favorable for more significant flavor transformation (e.g., [28,33,35]). In addition to predicting the presence of a crossing, we can follow the analysis behind the maximum entropy test further to estimate the properties of the crossing. The inverse cosine in Equation 17 yields two results, and we can use the difference between them to estimate the angular width of the crossing as ∆θ = 2 cos −1 η The bottom panel of Figure 4 shows the angular width of the crossing as determined from Equation 19 (bottom panel), along with an estimate of the same quantity extracted directly from the Monte Carlo results. We also estimate the angular width of the crossing from the Monte Carlo data using ∆θ MC ≈ 2 cos −1 1 − ∆Ω 2π (20) where ∆Ω is the solid angle occupied by the inverted portion of the ELN. The results are displayed in the top panel of Figure 4. In both cases, there is a broad region out to r ≈ 200 km that exhibits wide crossings (green ). The ME test qualitatively predicts the crossing depth in the disk, but predicts an artificially larg crossing depth in the polar regions. and yellow). However, there are significant differences in the structure. The ME estimation of the width in the polar region is significantly larger than the MC estimate. This is unsurprising given that analytic moment closures are known to perform poorly in polar regions. Outside of r ≈ 200 km, the ME test seems to indicate significantly wider crossings than those present in the MC data along radial structures. This structure reflects the structure apparent in the differences between flux factors between neutrinos and antineutrinos shown in the center panel of Figure 1, falsely correlating large electron neutrino flux factors with wide crossings. The growth rate of the FFI is sensitive to the depth of the crossing. Following [16], the growth rate (ω) can be estimated to scale with the "crossed" and "uncrossed" ELN densities I + and I − , defined as We calculate the "crossed" and "un-crossed" ELN densities from the discrete Monte Carlo data as where Θ is the Heaviside theta function. This is directly evaluated from the Monte Carlo data and displayed in the top panel of Figure 5. The deepest crossings are present in the dense part of the accretion disk, but crossings are present almost everywhere in the domain. As already described in several previous works, even the regions with a small crossing depth have growth rates that are much faster than the relevant advection or collisional timescales. While it is in general possible to integrate I + and I − for ME distributions, we instead approximate Equation 21 in a way that is more straightforward to evaluate in the context of a global two-moment radiation hydrodynamics simulation. We evaluate the ME distribution in the direction of the net ELN flux (purple vector in Figure 3) and in the opposite direction. Specifically, if I 1 is oriented with angle θ I1 , we evaluate Using, δf ME = f ME νe − f MĒ νe , We can then approximate the crossing depth as This is plotted in the bottom panel of Figure 5. Once again, the ME test reproduces the actual crossing depth rather well within the accretion disk. However, the ME test over-predicts the crossing depth in the polar regions, which is again expected due to the known problems of analytic closures in this region. However, since the growth rate everywhere is faster than other timescales, the particular growth rate is not as important is the presence of instability and the net flavor change it produces.
Finally, the total amount of eventual flavor change is related to the relative amount of electron neutrino and antineutrino excess (i.e. the relative sizes of I + and I − ). That is, if either the blue or red shaded regions in Figure 3 is small, there is not significant freedom for the fast flavor instability to transform overall flavor [28,35,69].
We evaluate the crossing ratio as These are plotted on the top and bottom panels, respectively, in Figure 6. For the case of the direct MC data (top panel), the ratio is close to unity (indicating possible significant flavor transformation) near the contour of n νe − nν e = 0 (black curve). The ME test predicts a large crossing ratio at similar locations, but at smaller radii. This estimate would likely be improved by a full angular integral of the ME test, but such an approach is likely too expensive to implement in global simulations of neutron star mergers. Note that the inability of much of the domain to undergo significant changes in flavor is consistent with the calculations by [70,71] using a related "instability parameter" and a toy model of flavor transformation in the merger system. The main results for this test are shown in the top left panel of Figure 7. The maximum entropy test predicts a crossing almost everywhere that one exists (green), never showing a false positive and predicting no crossing where one exists in only a few percent of the domain. Overall, the ME closure is quite good at predicting where ELN crossings are present and offers the ability to estimate qualitative details of the crossings. These details (crossing width, crossing depth, and crossing ratio) are approximately correct within the disk out to ∼ 200 km where most of the flavor transformation is expected to occur, but rather inaccurate in polar regions and at large radii.

C. Polynomial Test
The polynomial test of [52] states that if two different positive-weighted angular integrals of the ELN distribution G(Ω) have opposite sign, then the ELN itself must carry positive and negative values. This test has the advantage that if the full tower of angular moments is known, it can exactly predict the presence of an ELN crossing. In general, it also requires a sweep over parameter space, but we will show here that one can judiciously choose parameters to maximize the ability of polynomial tests to capture a crossing without requiring a parameter sweep. Unfortunately, polynomial tests also do not offer insight into the properties (wavenumber, growth rate) of unstable modes.
Specifically, for two different positive semidefinite functions F ± (Ω), the integrals are Note that the definitions of I ± are different from those ). The maximum entropy and resonant trajectory tests agree with the largest volume of Monte Carlo data, but all of the tests exhibit instability in the regions near the black contour where significant flavor transformation is possible (Figure 6). The tests that use information from the pressure and heat tensors are sensitive to the choice of closure. in Section II. The distribution is unstable if The weighting functions can then be chosen such that the integrals I ± are combinations of the known moments by using the form resulting in corresponding integrals of the form Following this process, one must search through all coefficients a, b i , etc. that make F ± positive for all Ω to see if any two combinations yield I ± of different signs. With moments up to rank 3, this is a 40-dimensional parameter space that is not practical to fully search. Instead, in the following we try to simplify the approach to extract straightforward polynomial expressions that can be applied without a numerical search procedure. We can capture the most significant effects by noting that, intuitively, the ELN is likely to have values of opposite sign in the directions along and opposite the direction of the net ELN flux I 1 . Thus, we can create a smaller class of weighting functions F ± (µ) = a + bµ + cµ 2 + dµ 3 + ...
with corresponding integrals where µ = Ω ·Î 1 and the scalaraized ELN moments are defined in Equation 5.
If we only use the number densities, the only meaningful weighting function is F + = F − = 1 (and positive scalar multiples), resulting in I + = I − = I 0 . Thus ELN crossings are guaranteed if Equation 27 is satisfied with which is equivalent to the α = 1 test of [57,58] and is plotted as a black curve in all cross-sectional figures.
As we explain in Section IV B in reference to Figure 6, this traces the regions most likely to experience significant flavor transformation. In core-collapse supernovae, this condition is satisfied in regions inside the protoneutron star where neutrinos are nearly in equilibrium with zero chemical potential, implying that even if the distributions are unstable to flavor mixing, the flavors are already effectively fully mixed. This is not the case in neutron star mergers, where there are significantly fewer heavy lepton neutrinos, leaving a great deal of room for significant flavor transformation.
If we now also use information from the number flux, we can choose weighting functions in the form of Equation 30 with only a and b nonzero. The impact of the flux information is maximized if we choose F + (−1) = 0 and F + (1) = 1 (for a "forward-weighted" integral), and (I 0 ± I * 1 ) .
The results of this test are plotted in the "Polynomial 1" panel of Figure 7. The region of predicted instability now obtains some spatial extent, largely encompassing the regions of large R crossing shown in Figure 6. Thus, although this method does not predict crossings everywhere they are present in the Monte Carlo data, it likely predicts crossings in the regions most important for flavor transformation, though this must be tested by nonlinear simulations. This treatment of the order-1 polynomial test is also exactly equivalent to a complete search, so there is no possibility that a different choice of coefficients of the energy density and fluxes would yield a larger unstable region. In addition, this test does not depend on the choice of closure, since it does not use information from the pressure or heat tensors.
If we now also allow use of the pressure tensor, c can be nonzero. The weighting function is maximally weighted to one side if we require F + (−1) = F − (1) = 0, F + (−1) = F − (1) = 0, and F + (1) = F − (−1) = 1. This requires that a ± = 1/4, b ± = ±1/2, and c ± = 1/4, and the resulting functions are plotted in green in Figure 8. The distribution is unstable if Equation 27 is satisfied with The results are displayed in the top panel of the "Polynomial 2" plot in Figure 7. The method is able to capture a bit more volume than the order-1 polynomial test, but not significantly so. In the bottom panel of the "Polynomial 2" plot, we replace I * 2 with that determined by applying the maximum entropy closure (Equation 12) to the number density and number flux of each species separately for each energy bin. This causes the crossing test to detect crossings in wings above and below the disk where the antineutrino flux factors are significantly larger than the neutrino flux factors (center panel of Figure 1). This indicates that the crossing test can be significantly influenced by the choice of closure.
We could also try to maximally weight the ELN distribution at µ = 0 (the "equatorial" region around I 1 ) by requiring F + (0) = F − (−1) = F − (1) = 0 and F + (−1) = F + (1) = F − (0) = 1. This constrains the coefficients to a ± = (1 ∓ 1)/2 and b ± = ±1. The corresponding weighting functions are plotted in green in Figure 8. The distribution is unstable if Equation 27 is satisfied with The results are shown in the top panel of the "Polynomial 2 (Equatorial)" plot in Figure 7. Applying the ME closure to the second moment (bottom panel) slightly further shrinks the volume of detected crossings. This method detects instability in only a subset of the region detected by other polynomial methods, and does not appear to be representative of any qualitative features of the FFI-unstable regions. However, it does affirm our choice to use polynomials that maximally emphasize the distribution along or opposite I 1 , as in the other polynomial tests. Finally, if we now allow use of the heat tensor, d can be nonzero. The weighting function is maximally weighted to one side if we require F + (−1) = F + (−1) = F + (−1) = F − (1) = F − (1) = F − (1) = 0 and F + (1) = F − (−1) = 1. This constrains the coefficients to a ± = 1/8, b ± = ±3/8, c ± = 3/8, and d ± = ±1/8. The corresponding weighting functions are plotted in gold in Figure 8. The distribution is unstable if Equation 27 is satisfied with The results are shown in the top panel of the "Polynomial 3" plot in Figure 7. As expected, the extra information allows the method to detect a slightly broader region of instability, the most significant differences being inside of the I 0 = 0 contour. For larger radii, the top panel of Figure 4 shows that the crossing angular width can be quite small, suggesting the need for very high-order polynomials to be able to detect them. It thus seems unlikely that polynomial crossing tests will be able to detect the full range of crossings in neutron star merger simulations. However, [22,23] already show that the FFI is present in a large fraction of the domain during simulations that dynamically include FFI-inspired flavor mixing. In addition, we note that all of the polynomial tests encompass the region where I + ≈ I − in which significant flavor transformation is possible, so the inability of the polynomial method to detect crossings at large radii does not necessarily make it significantly less realistic. An order-3 polynomial maximally weighting µ = 0 can be derived by requiring F ± (0) = 0, but the result is identical to F (2,eq) ± . D. k0 Test Dasgupta et al. 46 note that there is always a wave number k 0 for which the dispersion relation can be expressed as a function of only moments up to rank 2. This is a conservative test in that it cannot yield false positives, it provides insight into the growth rate of the k 0 mode, but it has the disadvantage that it cannot detect instability of any other mode. It requires finding roots of a matrix's characteristic polynomial, compared to transcendental equation solving in the Maximum Entropy test or parameter sweeping in some Polynomial tests, but the brevity of this section attests to the simplicity of the idea and the lack of free parameters to tune.
While we do not reproduce the derivation here, the frequency ω of this special mode can be determined by solving where V tt = I 0 , V ti = V it = I i 1 , and V ij = I ij 2 . The mode with wavenumber k 0 is unstable if there is a solution ω with nonzero imaginary component. The results of the test are shown in the top half of the k 0 panel of Figure 7, where it appears to perform very comparably to the Polynomial 2 test. This is perhaps not unexpected, as the two tests use information from moments up to the pressure tensor, even though they are sensitive to different unstable modes. The bottom half of the panel shows the results when the ME closure is used to provide the pressure tensor. Again, the results are very similar to the Polynomial 2 test that uses the ME closure. Although this suggests that simulations performed with different stability metrics may be consistent, this similarity might not be present in other realizations of NSM simulations or other systems (e.g., CCSNe). [26,72] showed that for isotropic modes (i.e. k = 0) and with axially symmetric distributions, if there is a direction that satisfies a particular resonance condition, then the distribution is unstable to the FFI. We do not repeat the derivation here, and the general distributions in this work are not axially symmetric. However, we nevertheless use the scalarized moments defined in Equation 5 to test how well this test performs in a realistic environment. This criterion states that the distribution is unstable to the FFI if

E. Resonant Trajectory Test
The above condition reduces to the instability criterion of [26] if axial symmetry is restored. This is an approximate test, since breaking the axial symmetry assumption allows for false positives, but this test performs surprisingly well. Despite its approximate nature, it misses only a small part of the domain known to be unstable and over-predicts instability in a few isolated pixels (top panel of the Resonant Trajectory plot in Figure 7). When I * 2 is determined using the maximum entropy closure, the results do not significantly change (bottom panel of the Resonant Trajectory plot in Figure 7). More testing would be required to determine whether this holds up in general at different points in time in the merger and in other situations like core-collapse supernovae, but its simplicity and apparent robustness against the choice of closure would make it very attractive to include in dynamical simulations.

F. Unstable Pendulum Test
Similar to the resonant trajectory test, [26,72] appeal to the pendulum-like nature of the evolution of angular moments of the neutrino radiation field, again under the assumptions of homogeneity and axial symmetry. This test suggests instability if Once again, this can be considered an approximate instability criterion, as there is no guarantee against false positives when the assumption of axial symmetry is broken. Like the resonant trajectory test, the unstable pendulum test performs better than expected (top panel of the Unstable Pendulum plot in Figure 7), although it misses a significant amount of instability in the polar regions. In addition, it appears to be more sensitive to the choice of closure, such that applying a closure (bottom panel) causes the test to falsely determine that disk some disk regions outside of the black contour are not unstable. We expect the resonant trajectory test to be more representative of instability than the unstable pendulum test, but more testing in other scenarios is required to see if this holds true in general.

V. CONCLUSIONS
We calculate a representative neutrino radiation field in a snapshot of a neutron star merger simulation using time-inedependent Monte Carlo radiation transport ( Figure 1). We use the results to show that there are electron lepton number crossings, and hence flavor instability, everywhere on the domtain except in the central hypermassive neutron star (Figure 2). We then take angular moments of this radiation field and assess how well a number of proposed tests are able to correctly determine the presence of ELN crossings using only these moments ( Figure 7). All of the methods predicted instability near the regions where significant flavor transformation is likely (Figure 6). The resonant trajectory test and the generalized maximum entropy test derived in this work predicte instability in almost all locations where ELN crossings are present in the full Monte Carlo data. Many of the tests showed significant dependence on the choice of closure, but the resonant trajectory showed remarkably little dependence, and the maximum entropy tests and order-1 polynomial tests are independent of the closure choice by construction. We note that each of these tests has particular advantages, including simplicity of implementation, guarantees to not over-predict instability, and insight into the growing modes of the distribution, and the optimal test to use varies by the need for each of these. In addition, while we chose a challenging and rich environment in which to test these instability metrics, the reader should be cautioned that the successes and similarities of the metrics may not carry through to other realizations of NSMs or other systems like CCSNe.
We generalized the maximum entropy test mentioned above in order to qualitatively predict the width (Figure 4), depth ( Figure 5), and relative size ( Figure 6) of the ELN crossing in addition to a binary determination of the presence of an ELN crossing. While the maximum entropy test is not able to quantitatively reproduce these quantities, much of the qualitative structure is reproduced. However, the need to solve a transcendental equation iteratively may make it too expensive to include in dynamical simulations.
Although many of the crossing tests are able to predict where instability occurs, they all remain rather agnostic about the amount of flavor transformation. More work in nonlinear simulation and analytic estimation of the final fixed point to which neutrinos relax after the instability (including coherent flavor waves [73]) are needed to address this deficiency. Finally, non-local effects that follow from neutrinos experiencing flavor instabilities in multiple parts of the domain can only be addressed with global simulations. Although much work remains to be done, we are hopeful that the results presented here guide the use and further development of the treatment of the FFI in simulations of neutron star mergers.