Neutrino Fast Flavor Instability in three dimensions for a Neutron Star Merger

The flavor evolution of neutrinos in core collapse supernovae and neutron star mergers is a critically important unsolved problem in astrophysics. Following the electron flavor evolution of the neutrino system is essential for calculating the thermodynamics of compact objects as well as the chemical elements they produce. Accurately accounting for flavor transformation in these environments is challenging for a number of reasons, including the large number of neutrinos involved, the small spatial scale of the oscillation, and the nonlinearity of the system. We take a step in addressing these issues by presenting a method which describes the neutrino fields in terms of angular moments. We apply our moment method to neutron star merger conditions and show it simulates fast flavor neutrino transformation in a region where this phenomenon is expected to occur. By comparing with particle-in-cell calculations we show that the moment method is able to capture the three phases of growth, saturation, and decoherence, and correctly predicts the lengthscale of the fastest growing fluctuations in the neutrino field.


Introduction
Neutron Star Mergers (NSMs) are dramatic multi messenger events producing a variety of observable signals.
The first direct evidence for the synthesis of the heaviest elements -the r-process or rapid-neutron-captureprocess elements -in neutron star mergers came from the optical/infrared emission (kilonova) associated with the first neutron star merger observed through gravitational waves, GW170817 [1][2][3][4][5][6][7][8][9][10].This kilonova light curve decayed on a relatively slow timescale, requiring the presence of lanthanides in the ejected matter [11,12].This confirmed many decades of prediction that neutron star mergers would make r-process elements [13,14].While the original proposal for r-process in NSMs required tidal stripping, there are now several additional types of ejecta that are proposed to make r-process elements in NSM, including the collisional ejecta [15,16], viscous evaporation [17] and magnetically- [18,19] or neutrino-driven winds [20,21].As the neutrinos affect not only the thermal evolution of the neutron star merger, but also the number of neu-trons available for capture in the ejecta [22,23] and therefore the shape and magnitude of the kilonova signal, e.g., [24,25], the investigation of the neutrino dynamics is of paramount importance.As neutrinos propagate through the media, their numbers, energies, momenta, and flavors will change as functions of time and position.In particular, an accurate accounting of the flavor of the neutrinos as a function of position and time is essential.
In some large scale dynamical simulations of NSMs neutrinos are already included, albeit classically, as they facilitate energy transfer and charge exchange in the evolving compact object.Much work has been done to characterize the transport processes which change number, energy, and momenta.There exists an extensive literature on using different techniques to capture neutrino transport in compact objects, taking inspiration from methods developed for photon and neutron transport.Moment methods [26][27][28][29][30][31][32] approximately reproduce the behavior of neutrino transport in many cases, and their efficiency enables simulations of a much larger scale than would be possible with direct Boltzmann transport (e.g., [32][33][34][35][36][37]).However, due to their approximate nature, moment methods always require validation against more exact grid-based [38,[38][39][40][41][42] or Monte Carlo [43][44][45][46][47][48][49][50] Boltzmann transport calculations.In this manuscript, we adapt an existing two-moment neutrino transport scheme to simulate the transport of neutrino flavor.
Not only do neutrinos change flavor due to the non-zero neutrino rest mass, the oscillations are modified by coherent forward scattering on neutrons, protons and electrons, i.e., the MSW effect [51][52][53] which operates on the outskirts of an NSM.In addition, deep in the interior of an NSM, coherent forward scattering on other neutrinos, the so-called neutrino self interaction effects, can play a dominant role [54][55][56].The transformations due to the self interaction can be roughly grouped into three classes: the matter-neutrino resonance [57,58], collective oscillations [59] and Fast Flavor Oscillations (FFO) [60].The importance of fast flavor transformation in mergers has been previously pointed out using parameterized outcomes of the FFO, e.g.[35,37,61,62].Multi-angle methods have been used to simulate the FFO directly in idealized scenarios or localized domains (e.g., [63][64][65][66][67][68][69]).However multi-angle methods are computationally expensive and the expense is a barrier to studying FFO in more realistic geometries.Using moments would significantly reduce the computational expense.But FFO are sensitive to the angular distributions -something the moments integrate over -so it is not clear whether moment methods are capable of capturing FFO.
In this work, we present the first simulations of quantum moment based neutrino FFO using neutrino distributions taken directly from a general-relativistic simulation of an NSM.We find that for these conditions, a moment based approach is capable of qualitatively capturing the local phenomenon and to determine its accuracy, we compare results with particle-in-cell calculations.

Evolution Equations and Methods
Our starting point is the result of the classical global general relativistic two-moment radiation hydrodynamics simulation of Ref. [70] which provides moments of the neutrino distributions for an NSM remnant.Quantum coherence of multi-neutrino wave functions require us to treat the system as an ensemble of n particles with a corresponding density matrix.If we ignore higher-order correlations by employing a mean-field approximation, we can construct dynamical quantities which only encode flavor coherence.(For investigations of the appropriateness of the mean field approximation see [71][72][73][74][75][76]). These quantities are the neutrino generalized density matrices, ϱ and ϱ which characterize the phase-space density for neutrinos and antineutrinos at time t, 3-position x, and 3-momentum p.When integrated over phase space, their diagonal components represent the number of neutrinos of a given flavor while the off-diagonal components give a measure of the coherence between two flavors.
References [77][78][79][80] derive evolution equations for ϱ and ϱ which are commonly known as the Quantum Kinetic Equations (QKEs) and involve a Hamiltonian-like opera-tor, where the dot over a quantity represents differentiation with respect to time.The left-hand-side (lhs) of Eqs. ( 1) and (2) give the evolution of the density matrices through phase space.We will set ṗ = 0 when we consider flat Minkowski space for our calculations.The rhs of Eqs. ( 1) and ( 2) gives a generic source term which governs the evolution of the density matrices.C and C are momentumchanging collision terms which we will neglect here.The Hamiltonian-like operators are sums of three terms: a vacuum potential which mixes flavor states during neutrino propagation; a matter potential which gives rise to forward scattering of neutrinos on background electrons, positrons, and nucleons; and a self-interacting potential H ν due to forward scattering on the other neutrinos, namely where G F = 1.166 × 10 −11 MeV −2 is the Fermi constant, p is a free variable for the density matrix in the commutators of Eqs. ( 1) and ( 2), and q is a dummy variable of integration.
The QKEs may be recast as a tower of evolution equations for the angular moments of the density matrices [81][82][83][84].The first three angular moments are more familiarly known as the energy density E, the flux vector F and the pressure tensor P .The explicit expressions for these moments are where p = |p|, subscripts a, b ∈ {e, µ, τ } represent flavor indices, and superscripts j, k ∈ {x, y, z} represent space indices.The factor of 4π in these definition is a normalization convention.As the phenomenon of fast flavor transformation is driven by the self-interaction, it is sufficient to write the first two moment evolution equations as keeping only the self-interacting term of the Hamiltonian while neglecting the matter and vacuum contributions.
These equations have implicit summations over repeated indices j and k using a Euclidean metric, and where The partial derivatives with respect to space denote advection and are key to capturing inhomogeneous modes of FFO.Similar equations can be derived for the evolution of the antineutrino energy density, E and flux vector F j .
The recasting of the QKEs as evolution equations for the moments is not an approximation: the moment evolution equations are simply a series of linear combinations of the QKEs with weighting factors proportional to the components of the neutrino momentum unit vectors.However, truncating the series of moment equations does involve an approximation since a priori we do not know the exact closure -the relationship between the lowest order unevolved moment [P jk in Eq. ( 8)] and the evolved moments.We algebraically solve for P jk by stipulating a specific closure relationship: the Maximum Entropy Closure (MEC).The MEC relies on an interpolation between the optically thin and thick limits (or more properly, neutrino thin and thick limits) where |F ab | is the magnitude of the flux vector for the ab component of ϱ and δ jk is the Kronecker delta function.Analogous expressions exist for anti-neutrino moments E, F j , P jk .The quantity χ (χ) is analogous to the classical Eddington factor, and we use the flavor traced flux factor in determining this quantity for all density matrix components [85].We calculate χ based on the MEC relation.(For more detail on the MEC in the classical limit, see [86,87].)For the flavor off-diagonal components of P jk ab , we evaluate Eq. ( 11) for the real part of P ab using the real part of f ab , and separately calculate the imaginary part of P ab using the imaginary part of f ab .Future improvements will be needed to remove the basis dependence that stems from this prescription.
Many codes already exist which contain the infrastructure for organizing moment variables (including Refs.[30,[88][89][90][91][92][93], among many others) and time-evolving those variables by calculating the advection terms on the left hand side of Eqs.(7) and (8).We choose the FLASH code [94] and amend it to include neutrino flavor transformation.Modifications include adding in new variables for the off-diagonal terms of the energy density and flux moments, and adding in new subroutines to calculate the commutators in Eqs.(7) and (8).Although we use energy density in our moment scheme, we only consider monoenergetic neutrinos in this work, implying the relationship between number and energy densities is simply  (1/p)E ab (x, p).For a basis of comparison to the FLASH results, we also do comparable QKE calculations with the Particle-In-Cell (PIC) quantum kinetics code EMU [95,96].
EMU follows the evolution of neutrino quantum states along discrete directions instead of evolving angular moments.
A necessary condition for fast flavor transformation is that the lepton number flux as a function of direction changes sign [97][98][99].We look for spatial points where these crossings occur using the values of the angular moments from the merger simulation [70], by first transforming from curvilinear coordinates into a comoving orthonormal tetrad to allow us to take advantage of the local flatness of the metric, rotating the coordinates to align the net electron lepton number flux with the z-axis, and then applying the maximum entropy crossing test of Ref. [100]. Figure 1 shows a polar slice of the merger simulation 5 ms after the collision, where shaded regions indicate instability.The black cross at (x, z) ∼ (25, 20) km shows the location from which we extract the initial conditions.We deliberately choose a point close to the disk where the flavor content of the neutrino field is likely to influence the wind nucleosynthesis [20,57,101].Table 1 gives the initial conditions for the neutrino moments at the black cross.
Emerging from the dense environment of the remnant, we would expect the neutrino worldlines to follow different paths from the site of the last scattering [102,103].The vacuum term in the QKEs will act to populate the flavor-off-diagonal components with different phases for the different neutrino worldlines.That associated potential is smaller than the self-interacting potential of Eqs. ( 9) and (10), where we estimate the ratio to be 10 −6 at the conditions present in Fig. 1.To incorporate phases into the off-diagonal components of our moment simulation, we perturb the imaginary and real components of the energy density moment by using random numbers scaled by the ratio 10 −6 of the maximum value of the flavor-diagonal components where −1 < A, B < 1, and a ̸ = b.For the flux moment, we implement the same perturbation as the energy density moment using flux factors δF ab (x) = δE ab (x)(Σ c N cc f cc )/ (Σ c N cc ), implying the off-diagonal components for E and F are correlated.Analogous expressions exists for the antineutrino moments.For the PIC method, we give the density matrices ϱ nonzero off-diagonal components with the same 10 −6 scaling and randomization for each particle and each off-diagonal element in each cell.If we calculate the energy density for the PIC method using Eq. ( 4), we would expect the incoherent sum for δE ab to be reduced by √ n for n particles per cell.As a result, we expect different behavior between the moment and PIC methods due to a different initialization of δE ab .Note that we see identical behavior with smaller perturbations since the perturbations grow exponentially.

Results
To solve the QKEs in both the moment and PIC methods, we use a 3-dimensional Cartesian geometry.The system is a cube with side length L = 7.87 cm.We divide the domain into 128 cells per side and institute periodic boundary conditions.Our choice of geometry allows us to study a range of spatial modes and resolve the fastest growing one, with little interference from boundary effects of a finite domain.For the flavor evolution, the moment method uses 2 flavors, e and a heavy-lepton x.
Figure 2 gives averaged components of the neutrino number density moment N ab plotted against t−t max , where t max is the time at the point of saturation (t max ≃ 0.5 ns for EMU simulations; t max ≃ 0.2 ns for FLASH).The solid dark red curve shows the results from the FLASH simulations, whereas the solid (dashed) black curve gives the results from the 2-flavor (3-flavor) EMU calculations.To test for convergence, we also give results for other simulations where we changed the domain parameters.Instead of using a cube with side-length 7.87 cm, the additional simulations have a side-length of 3.93 cm.The number of grid points per side also differs.The gray lines for both EMU calculations, and the medium red line for FLASH have 64 grid points per side, yielding the same resolution as the baseline case.The light red line for FLASH has 128 grid points per side, implying double the resolution of the baseline case.We observe that the FLASH calculations of differing domain size and resolution begin to diverge for t − t max ≃ 0.4 ns.This is due to an accumulation of error in the finite differencing of the advective term in Eqs.(7) and (8).Therefore, we focus our discussion on times before t − t max ≃ 0.4 ns, which includes the growth and saturation phases of the FFO.
We first examine the growth of the instability which can be seen as the abrupt drop of all curves in N ee at t − t max = 0 (top panel) and the corresponding rise in |N ex | (bottom panel).|N ex | grows from the initial perturbation until saturation, where the slope of the line is the exponential growth rate.We see excellent qualitative agreement between the moment and PIC methods, with our FLASH simulations having a slightly faster growth rate (8.09 × 10 10 s −1 ) than our EMU calculations (5.58 × 10 10 s −1 for 2-flavor).
Second, we examine the point of saturation, corresponding to the peak of ⟨|N ex |⟩ in the bottom panel.The magnitude of this peak is nearly identical for both methods.EMU has a larger duration of time where ⟨|N ex |⟩ is large, and we see this reflected in the top panel where the amplitude of the first oscillation in ⟨N ee ⟩ is larger for EMU than FLASH, indicating more flavor transformation in the former.
Third, we examine the decoherence phase, which corresponds to the leveling out of the ⟨N ee ⟩ curves in the top panel.We see that the electron number density agrees to about 30% for our choice of closure.For reference, the green lines show the value of ⟨N ee ⟩ if the system had reached a 2-flavor equilibration (solid) or 3-flavor equilibration (dashed).Examining the bottom panel, we see a larger loss of coherence in the FLASH calculation as well as the resulting freeze-out of ⟨N ee ⟩ at a smaller value as compared to EMU.
To show the behavior of the phase of the off-diagonal term ϕ ex , we plot 3D contours of the phase at 2 different time slices for the FLASH calculation in Fig. 3.The left panel shows the linear growth phase, where extended wavefronts of equal phase have begun to form.These waves grow larger in amplitude until the point of the FFO saturation.Although there is a net decoherence after saturation, the characteristic wavelength of the phase variations increases (right panel).This qualitative behavior is similar to what we see in PIC simulations, where we perform the calculation with both 2 and 3 flavors.[Top] Diagonal ee component plotted against t − tmax.Solid (dashed) black lines show the results for the 2-flavor (3-flavor) EMU baseline calculations.The solid red line shows the results for the FLASH baseline calculation.Also plotted are the results for the simulations with a smaller domain but identical resolution as gray and medium red lines.In addition, we give a light red line for a calculation with a smaller domain but double the resolution of the baseline calculation.Finally, we give a solid (dashed) green line to show the value of Nee  Fast flavor transformation is a phenomenon where the rapid growth occurs only on certain fluctuation scales as illustrated in Fig. 3. To determine this scale we show a Discrete Fourier Transform of the complex N ex component normalized by the flavor trace of N ab versus the magnitude of the wavenumber in Fig. 4.Both the EMU and FLASH calculations show power peaking at nearly the same wavevector k ∼ 5 cm −1 , corresponding to the fastest growing mode on a scale λ = 2π/k ∼ 1.0 cm.Thus the moment method effectively captures the fluctuation scale as well.There are some small differences that would likely change if a different closure were applied in the moment method.The peaks are slightly offset and there is a larger amount of power in the EMU calculations at the peak.In addition, there is more power in the FLASH calculations at larger wavenumber, which we expect to have a numerical origin and be intrinsic to the moment scheme.
The moment method relies on the MEC to reconstruct the flux distribution as a function of polar and azimuthal angle at all times.To compare the evolved angular distributions between FLASH and EMU, we show reconstructions of the flavor-diagonal density matrix distributions at a time t = 2 ns using mollwiede projections in Fig. 5.The top row gives cell-averaged angular distributions for FLASH with the MEC.Columns 1-4 are for electron neutrino, x-flavor neutrino, electron antineutrino, and x-flavor antineutrino, respectively.We plot dϵ ϵ 2 ϱ aa , where ϱ aa is averaged over the simulation domain and the integral over energy ϵ serves to specify the units.that the initial ELN crossing has been eradicated and that there is no longer a possibility of a FFO in both simulations.However, the most visible difference between the two sets of calculations is indeed the ELN, as there is a larger range of ELN values in the EMU calculations.This is the result of larger ν e fluences in the z direction, i.e., a smaller amount of flavor transformation.

Conclusions
To summarize our results, we addressed the open question of whether or not a moment method would exhibit fast flavor transformation in neutron star mergers by designing a moment method appropriate for quantum neutrino transport, and applying it to classical neutron star merger conditions.We found that fast flavor transformation did indeed occur although the moment method with MEC closure predicted more conversion of the electron neutrinos to other types but a similar level of mixing at the saturation peak.Specifically, the final electron (ondiagonal) neutrino number density in the two-flavor moment method with MEC closure asymptoted to 64% of its original value as compared with the two-flavor PIC calculation where it ended up at 74% of its original value.At the saturation peak, the moment method predicted |N ex | to be 16% of the original (on-diagonal) electron neutrino number density while the PIC calculation had 18%.Additionally, when compared to PIC calculations as a baseline, the moment method captured the growth rate to 45 percent and the wavelength of the fastest growing mode in the neutrino field to 33 percent.We attribute the differences between the moment method and the PIC method to our use of a classical closure which has previously been shown to overestimate the coherence of the moments [104].
When including neutrino transport, FLASH and other comparable state of the art codes take many millions of CPU hours for 3D calculations using classical moment based methods.Not only does a quantum moment flavor transformation method offer the promise of more seamless integration, but also the computational time required to model the FFO physics is a factor of 30 less for the baseline FLASH calculation as compared to EMU, which offers hope that a robust flavor transformation procedure can be coupled with a hydrodynamic code.
For a moment method to be a viable option, we need to sample additional points that have a wide range of neutrino angular distributions in order to test the limits of the moment method [85].In addition, future work will need to address several key issues.For example, the bottom panel of Fig. 2 shows divergences between different resolutions during the latter decoherence phase.The centimeter length scales simulated here are also still orders of magnitude smaller than the kilometer scales simulated in state of the art merger simulations.In addition, an investigation into quantum closures must be undertaken and the effect of differing closures [40,87,105] on the results quantified.Finally, in order to capture the flavor physics in NSM or core collapse supernova scenarios, the momentumchanging collision term [106-110] must be included.This is very challenging for computational-resource-heavy multiangle neutrino schemes, but the moment method we have presented here is a step toward achieving this milestone.

Fig. 1 .
The first three rows show the number densities of each neutrino flavor.For clarity, the third row shows the sum of all four heavy lepton anti/neutrino densities.Three-flavor simulations assume Nµµ = N µµ = Nττ = N τ τ = ΣN (x) /4, while two-flavor simulations assume Nxx = N xx = ΣN (x) /4.The next three rows show the flux factor vectors, the norm of which are the flux factors.The last row is the flux factor vector for a heavy-lepton in either the two or three flavor simulations.

Figure 1 :
Figure 1: Locations of a Lepton-Number crossing in a neutron star merger simulation 5 ms after merger as indicated by the presence of a zero crossing in the electron lepton number current.The radiation field energy density and flux come from the simulation of [70].Neutrino distributions are assumed to follow the angular distribution described by the maximum entropy closure, consistent with the original simulation.A colored region indicates that there exists at least one direction where the flow of electron-flavor neutrinos is equal to the flow of electron-flavor antineutrinos, equivalent to a zero crossing in electron lepton number as determined using the maximum entropy crossing test in[100].White regions indicate no such direction exists.The black cross indicates the location from which the parameters for calculation were drawn.From outside to inside, the green contours indicate matter densities of {10 11 , 10 12 , 10 13 , 10 14 } g cm −3 .

Figure 2 :
Figure 2: Spatially averaged components for the number density moment plotted against time measured from the point of saturation.[Top] Diagonal ee component plotted against t − tmax.Solid (dashed) black lines show the results for the 2-flavor (3-flavor) EMU baseline calculations.The solid red line shows the results for the FLASH baseline calculation.Also plotted are the results for the simulations with a smaller domain but identical resolution as gray and medium red lines.In addition, we give a light red line for a calculation with a smaller domain but double the resolution of the baseline calculation.Finally, we give a solid (dashed) green line to show the value of Nee for 2-flavor (3-flavor) equilibration.[Bottom] Magnitude of off-diagonal ex component plotted against t − tmax.Color and linestyle conventions same as top panel.For the 3-flavor EMU calculation, ex = eµ.The growth rate for EMU (2f) is 5.58 × 10 10 s −1 and for FLASH is 8.09 × 10 10 s −1 .
Figure 2: Spatially averaged components for the number density moment plotted against time measured from the point of saturation.[Top] Diagonal ee component plotted against t − tmax.Solid (dashed) black lines show the results for the 2-flavor (3-flavor) EMU baseline calculations.The solid red line shows the results for the FLASH baseline calculation.Also plotted are the results for the simulations with a smaller domain but identical resolution as gray and medium red lines.In addition, we give a light red line for a calculation with a smaller domain but double the resolution of the baseline calculation.Finally, we give a solid (dashed) green line to show the value of Nee for 2-flavor (3-flavor) equilibration.[Bottom] Magnitude of off-diagonal ex component plotted against t − tmax.Color and linestyle conventions same as top panel.For the 3-flavor EMU calculation, ex = eµ.The growth rate for EMU (2f) is 5.58 × 10 10 s −1 and for FLASH is 8.09 × 10 10 s −1 .

Figure 3 :
Figure 3: Volume rendering of contours of ϕex (the phase of the number density moment Nex) during the linear growth phase (left) and after saturation (right).The geometry is a cube with side length L = 7.87 cm.In the growth phase, approximately planar structures show wavefronts of the fastest growing mode.After saturation the phases chaotically mix.

Figure 4 :
Figure4: Discrete Fourier transform of Nex plotted against wavenumber k. k is the magnitude of the 3-vector.Linestyles and colors are the same as Fig.2.The time is taken as t − tmax = −0.1 ns for all three curves.The mode with the peak power for EMU (2f) is 4.79 cm −1 and for FLASH is 6.39 cm −1 .The rightmost end of the curves shows power at the grid scale.Lower resolution simulations (not shown) truncate at smaller k, but have a peak in the same location.

Figure 5 :
Figure 5: Energy-integrated and space-averaged angular distributions at t = 2 ns.The top row shows the distributions implied by applying the MEC to the simulated moments from the FLASH calculation, where each moment is averaged over the spatial domain.The bottom row shows the results from the EMU calculations, where each direction is averaged over the spatial domain and plotted as a colored dot.

Table 1 :
List of simulation parameters initially taken from Ref. [70] at the black cross in