Particle-in-cell Simulation of the Neutrino Fast Flavor Instability

Neutrinos drive core-collapse supernovae, launch outflows from neutron star merger accretion disks, and set the ratio of protons to neutrons in ejecta from both systems that generate heavy elements in the universe. Neutrinos of different flavors interact with matter differently, and much recent work has suggested that fast flavor instabilities are likely ubiquitous in both systems, but the final flavor content after the instability saturates has not been well understood. In this work we present particle-in-cell calculations which follow the evolution of all flavors of neutrinos and antineutrinos through saturation and kinematic decoherence. We conduct one-dimensional three-flavor simulations of neutrino quantum kinetics to demonstrate the outcome of this instability in a few example cases. We demonstrate the growth of both axially symmetric and asymmetric modes whose wavelength and growth rate match predictions from linear stability analysis. Finally, we vary the number density, flux magnitude, and flux direction of the neutrinos and antineutrinos and demonstrate that these factors modify both the growth rate and post-saturation neutrino flavor abundances. Weak electron lepton number (ELN) crossings in these simulations produce both slow growth of the instability and little difference between the flavor abundances in the initial and final states. In all of these calculations the same number of neutrinos and antineutrinos change flavor, making the least abundant between them the limiting factor for post-saturation flavor change. Many more simulations and multi-dimensional simulations are needed to fully probe the parameter space of the initial conditions.


I. INTRODUCTION
Supernovae represent the explosive birth and neutron star mergers the cataclysmic destruction of neutron stars (see [1][2][3][4][5][6][7] for recent reviews). Neutrinos are the dominant couriers of energy and lepton number in both cases and will be observable from the next galactic event. They carry the energy that explodes the star in a core-collapse supernova, drive outflows, and modify the ratio of neutrons to protons available for nucleosynthesis, thereby determining the abundances of elements that enrich the nearby universe. These heavy elements mix with ambient hydrogen that forms later generations of stars and planets.
Although state of the art models in many cases yield qualitatively correct explosion energies, ejecta masses, neutrino signals, electromagnetic transients, and gravitational waves (e.g., [8][9][10][11][12]), there remain many holes in the details. Simulations are growing increasingly sophisticated, but still suffer from low resolution (e.g., [13,14]) and uncertain initial conditions (e.g., [15][16][17]). Simulations also require as input an equation of state for matter beyond nuclear densities (e.g., [15,[18][19][20]) and the nuclear reaction rates of heavy and very unstable elements (e.g., [21]), both of which are poorly constrained. The final dominant model uncertainty is the treatment of neutrinos. The weak interactions between neutrinos and matter and other neutrinos cause them to be largely out * srichers@berkeley.edu of equilibrium, necessitating an expensive kinetic treatment. There is a natural trade-off between inexpensive approximate methods that allow for a larger number of simulations or a focus on specific aspects of the problem and high-accuracy methods that result in different explosion dynamics, ejecta properties, and neutrino signals [9,12,[22][23][24][25][26][27].
In addition, no global simulation dynamically treats the full effects of neutrino flavor transformation. Individual neutrinos exist in a quantum superposition of flavor states, so neutrino masses and potentials due to interactions with matter and other neutrinos can drive rapid and large transformations between quantum states. The neutrino-neutrino interaction term makes the evolution equations nonlinear, resulting in an exponentially difficult many body problem [28]. Because of this complexity the majority of the work in the field has been done under the mean field approximation (e.g., [29,30]), though the history and phenomenology is quite rich even with this approximation.
Vacuum and matter-induced oscillations were applied to core-collapse supernovae shortly after they were proposed as a solution to the solar neutrino problem [31][32][33] and are still used to map supernova simulation output spectra to estimate signals detectable at earth (e.g., [34]). So-called collective oscillations occur when the potential neutrinos feel from interactions with other neutrinos is comparable to the neutrino mass energy scale and has a long history in the core-collapse supernova context (see [35,36] for recent reviews). State of the art bulb-model simulations indicate that this transformation arXiv:2101.02745v1 [astro-ph.HE] 7 Jan 2021 mode is likely unimportant for the explosion mechanism (e.g., [37]), but is still impactful for ejecta nucleosynthesis (e.g., [38][39][40]). The matter-neutrino resonance can occur under certain conditions when the neutrino potential is comparable to the matter potential and may impact neutrino signals and nucleosynthesis from both supernovae [41] and neutron star mergers [42]. In addition, a small fraction of neutrinos emitted from an ongoing core-collapse supernova explosion will scatter off of heavy nuclei in the collapsing stellar envelope far outside of the shock front, forming a diffuse halo of scattered neutrinos, some of which are moving inward. These diffuse neutrinos can drive significant flavor changes in the much larger number of outgoing neutrinos [43][44][45].
Although proposed more than a decade ago [46], it was only appreciated in recent years that a new flavor instability could drive neutrino flavor transformation orders of magnitude more quickly than the previous transformation mechanisms [47][48][49][50][51][52][53]. While the magnitude of the effect of this fast flavor instability (FFI) on the abundances of each neutrino flavor is still uncertain, straightforward arguments suggest that the FFI should be a ubiquitous feature of both CCSNe and NSMs [45,54]. Although the associated small length and timescales make simulations of the full quantum kinetic equations (QKE) presently impossible in global multidimensional simulations, neutrino transport without flavor transformation is much more tractable. Many authors have used the neutrino distributions in these multidimensional neutrino transport simulations to diagnose the presence or lack of the FFI in CCSNe [40,45,[55][56][57][58][59][60][61][62][63] and NSMs [38,54,64]. Even if flavor transformation does not affect the dynamics, it has been suggested that it can significantly modify the elements formed in the ejecta, hampering the production of heavy elements in NSMs [38,64] and enhancing the production of light-p nuclei in neutrino-driven winds from CCSNe [40].
Motivated by the potential impact of the FFI and the hope that a general understanding of the final state of an unstable distribution can eventually be applied to global simulations, several authors have begun working on direct simulations of the instability. The majority of the simulations so far assume homogeneity ("one-zone" models) and are discretized in only angle (or angular moments) and time. These simulations have demonstrated a consistency between linear stability analysis and direct evolution of the nonlinear equations [65][66][67][68] and have shed some insight into the late-time angular turbulence and kinematic decoherence [69][70][71]. Other simulations of the FFI have included inhomogeneity [72,73], a simplified treatment of collisional processes [74], or both [50].
The capabilities of existing simulation methods are currently limited by various imposed symmetries, a reduced number of neutrino species, and/or a requirement for a very large number of grid cells or basis elements. To unify and expand on these models, a more general framework for simulating kinetics is needed. Neutrino transport methods such as moment methods [75][76][77][78][79][80][81][82], dis-crete ordinates [83,84], and Monte Carlo [85][86][87][88] could in principle all be extended to treat coherent flavor effects [89,90], but the strong dependence of the evolution of each neutrino on integrals of nearby neutrino distributions is challenging to implement efficiently.
In this paper, we utilize technology from the plasma physics community and describe a particle-in-cell method that formally solves the mean field quantum kinetic equations in an efficient, scalable manner. In Section II we outline a particle-in-cell implementation of the neutrino quantum kinetic equations. In Section III we demonstrate the exponential growth, saturation, and kinematic decoherence of a toy neutrino distribution on a onedimensional mesh. Finally, in Section IV we vary the neutrino distributions to begin the parameter study needed to build a sub-grid model of the FFI. We provide some concluding remarks in Section V.
We have developed the new code Emu [91] to implement this PIC method for solving the QKEs. Emu is fully open-source and is available at https://github. com/AMReX-Astro/Emu. All parameter files and select data from this study are publicly available [92] and further data is available upon request.

II. EMU: PIC NEUTRINO FLAVOR KINETICS
Emu solves the mean-field quantum kinetic equations without collisions (Section II A) by evolving the position and quantum state of a collection of computational particles moving through a background grid (Section II B). During each simulation timestep the particles aggregate their quantum states to construct a distribution within each grid cell (Section II B 1). The neutrino and background matter distribution are interpolated to each particle's position in order to construct time derivatives of the position and quantum state (Section II B 2). All particles are then integrated forward in time using a high-order integrator and a performance portable domain decomposed parallelization scheme using the AMReX framework (Section II C). This is repeated until the simulation is evolved for the desired amount of time. We discuss each of these steps in more detail below.

A. Quantum Kinetic Equations
The quantum kinetic equations that describe the transport of relativistic quantum particles read [29,30] ∂f ab ∂t (1) where N F is the number of neutrino flavors. Throughout this paper we use the convention that a, b, c, d are flavor indices, and are ∈ {e, µ, τ } for quantities in the flavor basis or ∈ {1, 2, 3} for quantities in the mass basis. The diagonals of f ab represent the occupation probability for neutrinos of flavor a located at position x and time t, moving with momentum Ω = p/|p| at approximately the speed of light c. The neutrino energy is 2 = p 2 c 2 + m 2 c 4 ≈ p 2 c 2 with energy hν. One can write a similar equation for the antineutrino distributionf ab that involves a distinct collision integralC ab and HamiltonianH ab . Throughout this work we neglect the collision terms C ab andC ab , which are weak enough to not significantly affect the distributions on the timescales simulated here. Implicit in this form of the equations is the assumption that there is no spin coherence, there are no right-handed neutrinos or left-handed antineutrinos, spacetime is flat, and the neutrino momentum is a constant of motion.
The Hamiltonian term H ab is also a N F × N F Hermitian matrix that encodes potential energy in the form of mass and interactions with other particles. It is usually broken down into a sum of the vacuum potential (due to the neutrino mass), the matter potential (due to interactions with non-neutrino particles), and the selfinteraction potential (due to interactions with other neutrinos). In the convention we use these can be written as a c 4 δ ab /2|p|. The first term is the same for all flavors and therefore cancels out in the commutator in Equation 1, so it can be ignored here. The magnitude of the momentum vector is |p| ≈ hν, so the vacuum Hamiltonian is usually written as H vacuum,ab ≈ m 2 a c 4 /2hν. The number density and flux are angular integrals over the neutrino distribution: ( Note that we denote a scalar lepton number density with n a and a Hermitian matrix neutrino number density with n ab . f ab is the (vector of Hermitian matrices) neutrino number flux density, not to be confused with the distribution function f ab or the lepton flux f a .n a andf a are real, but we leave the complex conjugation on them in Equation 2 for comparison to the neutrino potential. For antineutrinos,H vacuum,ab = H * vacuum,ab ,H matter,ab = −H * matter,ab , andH neutrino,ab = −H * neutrino,ab . Note also that there are other conventions in the literature that affect the form of the antineutrino Hamiltonians depending on whether one writes the evolution equation forf ab or f * ab .

B. PIC Algorithm
The left hand side of Equation 1 is a comoving derivative moving along the neutrino trajectory. Therefore rather than evaluating the derivative as fluxes between adjacent grid cells as is done in a finite volume scheme, we can simulate the distribution as a series of particles such that the advection term is accounted for by moving the particle's location. A computational particle carries two scalars N andN that represent the number of physical neutrinos and antineutrinos contained in the computational particle, along with two N F × N F Hermitian unittrace quantum density matrices ρ andρ that represent the quantum state of each neutrino and antineutrino. In practice, we only store and evolve the real and imaginary components of each element in the upper right half and the diagonal of these Hermitian matrices in the flavor basis. Wherever the lower left corners of Hermitian matrices are required, we express them in terms of the components of the upper right half in our implementation.
For a given particle labeled by index p, the equations of motion can then be expressed as Formally, these equations are only true in the limit of a large neutrino energy compared to the neutrino mass and potentials, but the same is true of Equation 1. In both forms, more terms would need to be included to faithfully treat low-energy neutrinos (e.g. [29]).
Because the neutrino self-interaction couples the flavor evolution of neutrinos and antineutrinos to their local distributions, we integrate Equation 4 for all particles together in step. To compute the right hand side of Equation 4 for each particle, we need to evaluate the Hamiltonian terms in Equation 2. The matter and selfinteraction potentials respectively depend on the local lepton and neutrino density and flux. Because evaluating the local neutrino density and flux requires determining the local neutrino distributions represented by the particles, we evaluate the Hamiltonian in two steps: deposition and interpolation. The background neutrino density stored in each grid cell must reflect the particle distributions themselves. A straightforward way would be to simply add up all of the particles within a cell, which would assume that each particle is shaped as a delta function in position and only interacts with particles in the same grid cell. However, we found that this causes the particle evolution to be unstable to artificially growing perturbations with wavelengths on the length scale of the grid resolution. We eliminated this numerical instability by instead modelling each particle with an extended shape as is usually done in PIC methods in order to smooth particle-mesh deposition and interpolation operations with higher order stencils.
The particle shape routines in Emu are copied directly from the open-source plasma physics code WarpX 1 [93] and include options from a delta function (zeroth order interpolant) to a cubic distribution (third order interpolant). We generally use the quadratic shape, since the cubic shape requires more ghost zones and higher-order interpolators do not necessarily correspond to higherorder convergence in PIC codes [94,95]. We integrate the total neutrino density and neutrino number flux at grid cell j as where w pj is the fraction of particle p that contributes to grid cell j according to the particle's shape function. N particles is the total number of particles in the simulation, N p is the number of physical neutrinos that particle p represents, and dV is the volume of the grid cell. In practice, we do not include the first term in order to minimize subtractive cancellation errors, since that term cancels out in the the commutator in Equation 4. We carry out the deposition in Equation 5 by first summing particle contributions to neutrino density and flux into local grids, including neighboring ghost cells. We next communicate ghost cell values to the MPI ranks containing the corresponding valid cells and sum the ghost cell contributions into the valid cells. We finally copy the grid data from the valid cells to all corresponding ghost cells. We apply periodic boundary conditions to ghost cells at the domain boundaries, mapping them to valid cells on the other side of the domain. Now that the lepton and neutrino number density and number flux are known at the center of each grid cell, we interpolate them from the grid to each particle's position x p using the particle shape function as in Equation 6: n a,p = Nzones j=1 w pj n a,j , n ab,p = Nzones j=1 w pj n ab,j , This yields the lepton density (n a,p ), the neutrino densities (n ab,p andn ab,p ) and fluxes (f ab,p andf ab,p ). We then obtain the neutrino and antineutrino Hamiltonians H ab,p andH ab,p for each particle by applying the particle state and these interpolated quantities to Equation 2. We finally use Equation 4 to calculate the particle's density matrix time derivatives ∂ρ ab,p /∂t and ∂ρ ab,p /∂t.

C. PIC Implementation
We implement this algorithm in modern C++ using the exascale computing framework AMReX [96], which provides high level abstractions for domain decomposition and both distributed and shared-memory parallelization. We thus represent the domain as a union of nonoverlapping rectangular blocks, each containing particles and Cartesian grid data. These blocks are distributed across MPI ranks, which allocate memory for their local grid and particle data. Because higher order particle shape functions require particles to deposit and interpolate from neighboring cells, each block possesses a set of ghost cells bounding the valid cells in its interior. We choose the number of ghost cells consistent with the stencil width of the particle shape function desired (e.g. 1 ghost cell for quadratic shape functions with a stencil size of 3).
We use a fourth-order explicit Runge-Kutta method to integrate Equation 4 for all particles. We integrate particle positions concurrently with their flavor states, so we do not split the advection and flavor-changing operators in Equation 1. Thus, when constructing the state at which to evaluate a Runge-Kutta stage right hand side, we first calculate particle positions along with their density matrices at that time. As a result, a particle may move from one block on the domain to another, which may reside on a different MPI rank. Such a particle will therefore need access to grid data on a different block before we can evaluate Equation 4 for this Runge-Kutta stage. Thus, immediately after calculating the new particle positions and density matrices corresponding to a Runge-Kutta stage, if any particle position now resides within a different block we redistribute its data for all Runge-Kutta stages to its new block and MPI rank. We then resume calculating the right hand side defined by Equation 4 for the next Runge-Kutta stage using the deposition and interpolation operations discussed in Section II B 1 and Section II B 2. We repeat this redistribution procedure at all Runge-Kutta stages and at the end of the time step. This ensures that we are able to evaluate all Runge-Kutta stage right hand sides for the particles even if previous stage right hand side terms were originally evaluated on different MPI ranks.
Local particle arithmetic operations (deposition, interpolation, evolution) constitute the vast majority of the computational cost. We implement these kernels using AMReX's performance-portable strategy for targeting either CPU threading (with OpenMP) or GPUs (with CUDA). The quantities evolved in Equation 4 are arranged in a struct for each particle and each block stores an array of these particle structs. For CPUs, we distribute blocks with particles to threads using OpenMP. Each thread loops over particles within the block, though we find in practice that increasing the number of MPI ranks can give better performance on CPUs than OpenMP threading. For GPUs, we assign one MPI rank to each GPU and allocate local grid and particle memory using CUDA Unified Memory. Each MPI rank loops over its local blocks and asynchronously launches a CUDA kernel parallelizing the loop over particles within the block. We then synchronize GPU kernels only after looping over all local blocks to ensure high kernel occupancy.

III. RESULTS
In this work we restrict ourselves to simulations performed in one spatial dimension and two angular dimensions. That is, we arbitrarily assume all neutrinos have an energy of 50 MeV, though this only influences the vacuum Hamiltonian, which has negligible effect on the timescale of these simulations. In this section we examine in detail only a single choice of initial conditions. This simplified setup will serve as a basis for comparison with other simulations in this and future work. For our fiducial simulation we choose initial conditions that exhibit a strong electron lepton number (ELN) crossing (i.e., a direction along which there are equal numbers of electron neutrinos and antineutrinos), characterized by electron neutrino and electron antineutrino distributions with a flux factor of 1/3 in opposite directions. The distribution is thus quite unstable to the FFI. We randomly perturb the initial conditions to seed all modes, as we will describe in Section III A. The amplitudes of the unstable modes grow exponentially precisely in the way predicted by linear stability analysis as shown in Section III B. The instability then saturates, diffusing neutrino flavor fluctuations away from the fastest growing wavelength into modes of all wavelengths as described in Section III C. Select data and all parameter files from this study are publicly available [92]. Additional data is available upon request.

A. Simulation Setup
In our fiducial simulation we use a simulation domain 64 cm in length with a grid resolution of 1 × 1 × 1024. The unit size in the x and y directions makes the simulation behave in a one-dimensional manner. That is, the neutrino distributions are homogeneous in the x and y directions, but are not restricted to isotropy or axial symmetry. We place many particles at the center of each grid cell with directions distributed approximately isotropically using the "Polar Coordinate Subdivision" method of Ref.
[97], except that we specify the number of points along the equator rather than the total number of points. Requesting 16 directions in thex −ŷ plane results in 92 particles per grid cell. In contrast, PIC codes often distribute particles initially randomly. We elect to distribute particles regularly so we do not have to consider Monte Carlo numerical noise when interpreting the results.
We arbitrarily choose an initial neutrino distribution with a number density of n = 4.89 × 10 32 cm −3 of each electron neutrinos and electron antineutrinos, and zero density for heavy lepton neutrinos. The fastest growing mode of a two-beam model with these neutrino densities has a wavelength of 1 cm, allowing a more straightforward comparison to the two-beam test in Appendix A 4. To create the desired anisotropic neutrino distribution using an isotropic distribution of computational particles, we give each particle weights (i.e. the number of physical neutrinos or antineutrinos that a computational particle represents) of where n (n) is the input number density of (anti-)neutrinos and f (f ) is the input number flux of (anti-) neutrinos. We use f = nẑ/3 andf = −nẑ/3 for this fiducial simulation. This definesẑ as our axial symmetry axis and results in a healthy ELN crossing in thex −ŷ plane. The initial state of each particle is set to be in a nearly pure electron flavor state with a random perturbation in the off-diagonal elements. That is, U ∈ [−1, 1] is a uniform random number generated individually each time it appears and α = 10 −6 is the strength of the random perturbation. µ and τ are determined after the random numbers are generated in order to ensure unit trace and a flavor vector length of one.
Since the initial distributions have almost no ν τ or ν µ content, we set the ρ µτ off-diagonal components to zero.
The left column of Figures 1-3 show these initial conditions at t = 0. The top left panel of Figure 1 shows the initial number density of each neutrino species plotted over 16 cm of the 64 cm domain divided by the trace of the number density matrix. All neutrinos (solid curves) and antineutrinos (dashed curves) begin in the electron flavor state (blue), and hence the electron neutrinos and antineutrinos constitute essentially all of the number density. All neutrino flux is in the ±ẑ direction, so the x component of the flux shown in the bottom panel at t = 0 is zero everywhere. The left column of Figures 2 and 3 show the magnitude and complex phase, respectively, of the eµ (purple), eτ (green), and µτ (orange) components of the neutrino density. The perturbations are too small to be visible in the top left panel of Figure 2. Although the perturbations of the density matrix of each particle individually is O(10 −6 ), they average together to make the perturbation of the x component of the flux factor (bottom panel) much smaller than 10 −6 . The left panels of Figure 3 simply show a random distribution of phases of the eµ and eτ components of both the number density and flux flavor matrices. The complex phase of the µτ component is shown as zero because it has exactly zero magnitude.

B. Linear Growth Phase
Before discussing the simulation results, we briefly go through the standard linear stability analysis to point out what modes should be unstable and the rate at which they are expected to grow. Following [49], we calculate the imaginary component of the frequency for modes of the form N ex = αe −i(ωt−kz) described by a frequency ω, wavenumber k, and initial perturbation amplitude α. x here represents either µ or τ . This is a two-flavor analysis, so the results are applicable to e − µ and e − τ flavor transformation but does not describe µ − τ transformation [98]. The linear dependence of the distribution on the polar angle makes computing the integrals in the polarization tensor analytic, though a numerical root find still needs to be performed to extract the growth rate of each mode. We plot the results in Figure 4. Since the results scale with the neutrino density (assuming the flux factors and directions are held constant), we plot both the growth rate and wavenumber scaled by √ 2G F n νe as a metric of the self-interaction potential. There are two branches, one of which has a larger growth rate than the other at all wavenumbers. The faster growing mode (blue) has a peak growth rate of Im(ω) = 6.50 × 10 10 s −1 and a wavelength of 2.20 cm. This is still somewhat slower and longer wavelength than the fastest growing mode in a two-beam model with equivalent densities (Appendix A 4). The slower growing mode has a maximum growth rate of Im(ω) = 1.06 × 10 10 s −1 at a wavelength of 4.45 cm. We will later show that the blue curve rep-resents the axial symmetry preserving mode and the red curve represents the axial symmetry breaking mode.
Note that some of the unstable modes have wavelengths longer than can fit on our domain (shaded gray region), but the growth rate of these modes is significantly smaller than the peak growth rate. In addition, the finite number of grid points limits the shortest possible wavelength to 0.125 cm, or a maximum wavenumber k/ √ 2G F n νe ≈ 16, so that all modes shown in Figure 4 outside the gray shaded region are within the grid's Nyquist limit. More importantly, our finite resolution limits the number of grid cells supporting each wavelength, and we find that to resolve the instability we need many grid cells per wavelength of the fastest growing mode. We explore the resolution requirements to obtain a converged instability growth rate in Appendix B and verify that the quoted results of our simulation are robust against changes to the domain size and grid resolution.
The simulation immediately exhibits rapid growth of the flavor off-diagonal components of the neutrino density matrix. As a metric of the net growth of all modes, we calculate the magnitude of the off-diagonal components of the neutrino density and the x component of the flux averaged over the entire domain as |n ab | probes the presence of modes that preserve axial symmetry, while |f ab | probes modes that break axial symmetry because |f (x) ab | = 0 would result from axially symmetric modes alone. The purple and green solid curves in Figure 5 show growth of the eµ and eτ components of the neutrino density. For the first ∼ 0.1 ns they both exhibit rapid, but not precisely exponential growth. In this phase, there are several competing modes growing at different rates, all of which are initially large enough to contribute to |n eµ | and |n eτ | . After 0.1 ns the average amplitudes grow exponentially until they saturate at ∼ 0.27 ns. After saturation, the average flavor offdiagonal densities hover around the saturation level of a few tenths of the total neutrino density.
The top panel of the center column of Figure 1 shows the flavor diagonal components of the neutrino density at the end of the linear phase. Electron neutrinos are still dominant, though there is now a significant fraction of ν µ and ν τ , showing that the instability indeed leads to significant flavor transformation. The amplitude of the x-component of the flux has grown as well, though it remains at only ∼ 10 −7 of the total neutrino density. Although during the linear growth phase the anti-neutrino curves (dashed) match the neutrino curves (solid), as the simulation approaches saturation the nonlinear interaction between modes drives chaotic behavior and the neutrino and antineutrino distributions begin to show slight  The blue curve corresponds to modes that do not break axial symmetry, while the red curve corresponds to modes that do. Both are normalized by a characteristic selfinteraction potential so the diagram is independent of neutrino density as long as the self-interaction potential is much larger than the vacuum potential. Unstable modes exist for wavelengths λ > ∼ 0.75 cm (k/ √ 2GF nν e < ∼ 2.66), though the growth rate for the homogeneous mode (k = 0) is also 0. The gray shadow shows the region of wavenumbers where the corresponding mode cannot exist on our domain because of our limited domain size of 64 cm.
Time evolution of domain-averaged absolute-value of the neutrino number density (solid) and number flux in thex direction (dashed). The labeled line segments show the slopes corresponding to the growth rates of the fastest growing axially symmetric (blue) and axially asymmetric modes (red) shown in Figure 4. The eµ (purple) and eτ (green) components of the number density grow at approximately the axially symmetric maximum growth rate, and the number fluxes grow at approximately the axially asymmetric growth rate.
The µτ (yellow) components grow parasitically at the same rate as the µµ (brown) and τ τ (salmon) components.
in Figure 4. The antineutrinos (dashed) are similar, except that they have a negative wavenumber, resulting in phase decreasing with increasing z. The bottom panel shows the complex phase of the flavor off-diagonal x components of the neutrino flux. Unlike the number density, this quantity is sensitive to axial symmetry breaking and is not sensitive to modes that preserve axial symmetry. The wavelength is consistent with the 4.45 cm prediction from the peak of the red curve in Figure 4. Thus, we expect that the blue curve in Figure 4 represents the symmetry-preserving mode and the red curve represents the symmetry-breaking mode. The mode identification is further corroborated by the growth rates of both modes. The growth rate of the flavor off-diagonal number density (solid purple and green curves) and x-flux (dashed purple and green curves) in Figure 5 are consistent with the growth rates predicted by the peaks of the blue and red curves, respectively, shown in Figure 4. The simulations do, however, show a slightly slower growth rate than predicted by the peak of the dispersion relation. This may be partially due to the fact that the domain-averaged quantities reflect many modes, all but the fastest of which grow slower than the fastest growing mode. However, the simulated growth rate also does depend somewhat on numerical considerations. We perform a convergence study in Appendix B to show this dependence.
The µ−τ mixing is more difficult to interpret. The n µτ component (solid yellow curve in Figure 5) grows with twice the growth rate of n eµ and n eτ . We expect that this simply follows the growth of the n µµ (brown) and n τ τ (salmon) elements, each of which also rise with this doubled growth rate. This can be understood intuitively in terms of a rotating flavor vector ρ in SU(3) flavor space. The eight components of this vector are where λ ab are Gell-Mann matrices and i is the vector index. Flavor transformation rotates this vector, preserving Tr(ρ ab ) = 1 with constant | ρ|. In this case, | ρ| = 1/ √ 3 because we start in a nearly pure electron flavor state in Equation 8. Combining these two conditions yields a relationship If, as in the linear growth phase of the simulation, the neutrinos are only slightly perturbed from the electron flavor state, the second term on the right hand side is much smaller than the first. Additionally, ρ µτ is much smaller than ρ eµ and ρ eτ throughout the linear growth phase. Combining these approximations, where ρ 2 ex = ρ 2 eµ + ρ 2 eτ and ρ xx = ρ µµ + ρ τ τ . Thus, if ρ ex ∝ exp(Im(ω)t) (where x refers to either µ or τ ) then ρ xx ∝ exp(2 · Im(ω)t). Since the growth of ρ µτ requires non-zero difference between ρ µµ and ρ τ τ , the growth of the initial random differences between ρ µµ and ρ τ τ give ρ µτ the room to grow at the same rate as ρ µµ and ρ τ τ . Note that this shows that the growth of n µτ can be up to twice that of n eµ and n eτ but does not show that it must be so. It also does not explain the rapid growth of all fluxes except but the f eµ and f eτ . More work is required for a complete description of the growth of all modes and flavors (though see [98] for work in this direction).
The differing growth rates of modes of different wavelengths can be seen in the evolution of the Fourier transforms of neutrino distribution shown in Figure 6. Each curve in the plot is the Fourier transform of a grid quantity time-averaged over a 0.1 ns interval to reduce noise in the plot. The color indicates the end of the time interval for which the transform was evaluated. Beginning with the middle panel, the bottom-most curve shows the Fourier transform of n eµ from 0 − 0.1 ns and the next curve (salmon) shows the same for the interval 0.1 − 0.2 ns. We use salmon for the second curve to highlight the time interval squarely within the linear growth phase. The wavenumber of the peak, the minimum unstable wavenumber, and the maximum unstable wavenumber are at the locations predicted by the blue curve in Figure 4, as both of these intervals lie within the linear growth phase of the instability, . The third curve shows the time interval 0.2 − 0.3 ns, during which the instability saturates. The peak location is consistent with the wavenumber of the fastest growing modes, but additional features develop away from the peak that we will discuss in Section III C. The top panel of the figure shows the Fourier transform of n ee on the same time intervals. Since this component is purely real, the Fourier transform is symmetric in k, but still reflects the characteristic wavelengths in the unstable modes. The bottom panel shows the Fourier transform of f (x) eµ . Once again, the location and width of the peak in the salmon curve matches the expectations from the symmetry-breaking (red) branch of Figure 4.
The symmetry breaking and preserving modes are also apparent in the angular structure of the neutrino flavor. The top panel of Figure 7 shows a Mollweide projection of the directions of all neutrinos located within a single grid cell for a simulation with high angular resolution. Each point on the plot represents the direction of an individual computational particle. In this high-resolution simulation we request 128 directions in thex −ŷ plane, visible as 128 points along the equator in the plot, for a total of 6022 particles per cell. Our fiducial simulation uses only 16 points along the equator, but we use the high angular resolution data for better visualization. We show in Appendix B that the growth rates and saturation states are the same in both simulations. The colors of the points show the real part of ρ eµ for each packet at t = 0.27 ns. Although the FFI has greatly amplified these numbers from their initially random values of O(10 −6 ), the flavor structure is very axially symmetric. If we subtract the axial average (bottom panel), more of the axial symmetry breaking mode is so much slower. By the time the fastest growing axially symmetric mode saturates, the axially asymmetric mode has only grown by a factor of exp(Im(ω max )∆t) ≈ 20. Overall, the simulations confirm the predictions of the linear growth phase from linear stability analysis. Modes of the expected wavenumbers all grow at the expected rates. The modes with the fastest growth rates appear to be the most important, as they saturate and lead to nonlinear evolution long before the slower modes can build significant amplitude.

C. Saturated Phase
At t ≈ 0.27 ns the exponential growth of both the axially symmetric and asymmetric modes shown in Figure 5 cease as the instability saturates. After saturation, the number density of all three neutrino flavors tends toward equipartition at late times. The right column of Figure 1 shows a snapshot of the flavor diagonal elements of n ab (top) and f (x) ab (bottom) at t = 5 ns. In both cases, the flavor content is widely variable in z, but when averaged over the domain, the neutrino density and flux are very close to flavor equipartition. The right column of Figure 2 shows that the flavor off-diagonal components of the density and flux remain large compared to the flavor diagonal components. Finally, the right column of Figure 3 shows that the phase is effectively randomized in the saturation phase, indicating that the coherent modes that grew in the linear phase are all but disrupted. Figures 1-3 also show variation on much shorter wavelengths (higher wave numbers) than the fastest growing modes that dominated the linear growth phase. This evolution to short wavelengths can be seen in Figure 6, which shows Fourier transforms of n ee , n eµ and f (x) eµ from the N z = 2048 simulation in Appendix B. We display this higher-resolution calculation because it allows us to track the evolution out to larger wavenumbers than our fiducial simulation, even though the instability growth rates and final averaged neutrino flavor abundances are the same as in the fiducial simulation. We described the amplification of the unstable modes shown in the earliest two curves (the lowest purple curve and the salmon curve) in Section III B. By 0.3 ns (third curve from the bottom of each panel, though it is not clearly visible in the lowest panel) the amplitudes of the fastest growing modes have reached their maximal values and power in all three quantities begins to spread to larger and smaller wavenumbers. Though the noise in the flux spectra (bottom panel) makes them difficult to interpret, each curve of | n ee | and | n eµ | varies linearly with |k| 1 ns/t in the logarithmic plot, indicating that | n ee | ∼ | n eµ | ∼ exp −|k| t 10 −9 s −1/2 .
By t = 5 ns the base of the wings in Figure 6 have reached a wavenumber of |k| ≈ 100 cm −1 , which is the maximum representable wavenumber on our grid of 2048 spatial points. Power begins to artificially build up at these higher wavenumbers, resulting in the slight upswing in the red curves at k 1 ns/t ≈ 40 for the final (red) curve. In the lower resolution fiducial simulation this upswing occurs after t ≈ 2 ns when the wings reach that grid's maximum k value of 50 cm −1 . This suggests that in nature a true steady state solution will only be achieved by waiting long enough for the power to diffuse from the unstable wavelength to the inter-particle distance (representing the maximum physical wavenumber) multiple times, such that the distribution in wavenumber can equilibrate. However, the process of kinematic decoherence likely behaves differently in nature's three dimensions than in our simulation's single spatial dimension, so the decoherence rates demonstrated in this simulation are not necessarily realistic.
Finally, we show the angular structure of the postsaturation distribution in Figure 8 corresponding to t = 4.95 ns. Similarly to Figure 7, the top panel of Figure 8 shows Re(ρ eµ ) for each particle in a single grid cell, where particle directions are plotted using a Mollweide projection. Amazingly, the saturated distribution remainsx Although these modes are initially unstable, the axially symmetric modes are much faster and saturate while the asymmetric modes still have small amplitude. This saturation seems not only to stop the growth of the symmetric modes, but also the asymmetric modes. However, this picture also likely changes in three spatial dimensions, so three dimensional simulations will be required for a more complete understanding of axial symmetry breaking. In addition, because the initial conditions in this simulation are themselves axially symmetric, the presence of asymmetries in the initial conditions would likely make this process rich in phenomena not expressed in our fiducial simulation.

IV. PARAMETER STUDY
Our final goal in this work is to provide a first look into the dependence of the FFI on the initial matter and neutrino distributions. The initial conditions described in Section III A rely on nine input parameters: the number density of neutrinos n and antineutrinosn, the number flux vector of neutrinos f and antineutrinosf (three components each), and the background electron density n e . Although we have shown that Emu simulates the FFI with high fidelity, it is impractical to densely sample the full parameter space, even under this limited subset of possible initial conditions. Instead, we scratch the surface here by performing a parameter sweep in a single input variable at a time, holding the other input parameters the same as in the fiducial simulation.
The results of these simulations are summarized in Figure 9, where we plot the time evolution of the domain average abundance of electron neutrinos. Each curve is a distinct simulation with the input parameter varied from the fiducial simulation indicated by the color. In the top panel, we show the results of seven simulations identical to the fiducial simulation except that we have included a nonzero background electron density n e . The fiducial simulation has n e = 0 and is plotted as the darkest curve, though it is occulted by the other curves. As described in Section III, the instability saturates after ∼ 0.3 ns, resulting in a significant decrease in the number of electron neutrinos. At late times, the number density of electron neutrinos (solid curve) and antineutrinos (dashed curve) approaches one third of the number density (green line) of all neutrinos or antineutrinos, respectively, indicating flavor equilibration. As expected from the stability analysis, the electrons modify the real component of the mode frequencies but do not affect the growth rate of the unstable modes. Thus, the simulations varying n e proceed very similarly, diverging at later times due to the chaotic nature of the saturation phase.
In the second panel from the top we show a series of eight simulations where we hold the number densities and flux directions of neutrinos and antineutrinos constant, but vary the magnitude of the antineutrino flux from 0 (isotropic, black curve) ton/3 (fiducial simulation, light orange curve). In the third panel from the top we show results from a series of seven simulations where we hold the neutrino densities and flux magnitudes constant, but change the direction of the antineutrino flux from aligned with the neutrino flux (black curve) to opposed to the neutrino flux (fiducial simulation, light orange curve). In both cases, the parameter variation reduces the growth rate of the instability, increasing the time when the instability saturates. However, these parameters do not appear to change the final abundances of neutrino flavors. All simulations approach flavor equilibration except for the case where the antineutrino and neutrino fluxes align and no flavor transformation occurs. Also note that when we vary either the magnitude or direction of antineutrino flux, the neutrinos and antineutrinos follow the same tra- In black, we show the results of simulations where we set f andf and varyn/n from 0−1 (the same data as in the bottom panel of Figure 9). The fiducial simulation is plotted as a purple point, and the horizontal green line marks complete flavor equilibration ( nee = Tr(n ab )/3). In both parameter sweeps, the antineutrinos consistently experience a greater degree of relative flavor transformation.
rapid rate than antineutrinos. When there are no antineutrinos present (uppermost black line) the distribution shows no flavor transformation. However, the fact that only the ratio of antineutrino to neutrino density changes the final flavor abundances is an artifact of our choice of initial conditions for our fiducial simulation. We performed another set of simulations with initial parametersn/n = 2/3, f = 0, and f = −ẑ and varied |f | from 0 ton/3. Using these parameters, there is an ELN crossing only for 3|f |/n ≥ 1/2. We run each simulation for 25 ns to ensure that all unstable simulations had sufficient time to reach an approximately stationary final state. In Figure 10 we plot time averages of the electron neutrino (solid) and antineutrino (dashed) abundances in blue. As expected, when there is no ELN crossing (3|f |/n < 1/2) the distribution exhibits no flavor transformation. Beyond this point, increasing the flux factor increases the size of the ELN crossing, increases the growth rate of the instability, and increases the amount of eventual flavor transformation. For none of our simulations in this parameter sweep do we see complete flavor equilibration ( n ee = Tr(n ab )/3), though this may occur for larger antineutrino flux factors. For comparison, we also plot the same quantities for then/n parameter sweep (i.e., calculated from the results shown in the bottom panel of Figure 9) in black and label the fiducial simulation in purple. In both parameter sweeps the antineutrinos consistently experience a greater degree of relative flavor transformation. This is simply a reflection of the fact that in all of our simulations, the same number of neutrinos and antineutrinos oscillate. That is, we find that in our simulations n − n ee =n − n ee .
As elsewhere in the paper, all variations of n in the equation are number densities with units cm −3 . Since there are fewer antineutrinos than neutrinos, a larger fraction must change flavor to match the total number of oscillating neutrinos. This is not entirely surprising. Since the neutrino self-interaction potential dominates the total potential, neutrinos and antineutrinos experience the same potential (up to a sign) and the flavor vectors of an individual particle rotate the same amount for neutrinos and antineutrinos, though in opposite directions. Equation 14 does not predict the amount of neutrinos that oscillate, but if this result holds in general for fast flavor transformations, it does say that the net electron flavor content (i.e. number of electron neutrinos minus the number of antineutrinos) does not change under the fast flavor instability even in saturation. In addition, it provides a bound on the amount of neutrinos that change flavor. In the limit of zero heavy lepton neutrinos, n − n ee = min(n ee ,n ee ). We can extend this to three flavors, noting that the amplitude of possible flavor transformation is given by √ 3l, where l is the length of the SU(3) flavor isospin vector (analogous to Equation 10 for n ab instead of ρ ab ). That is, n − n ee ≤ √ 3 min(l,l) (15) and similarly for antineutrinos. Explicitly, when the flavor off-diagonal elements are zero, these can be written l 2 = 1 6 (n νe − n νµ ) 2 + (n νe − n ντ ) 2 + (n νµ − n ντ ) 2 (16) and similarly for antineutrinos. If the antineutrino density flavor vector length is smaller than that of the neutrinos, the flavor rotation of the antineutrinos is not bound, but the flavor rotation of neutrinos must match. Complete flavor equilibration of both neutrinos and antineutrinos is possible only if both flavor vector lengths are equal. Once again, these bounds do not actually predict the final abundances; the abundances of electron neutrinos in Figure 10 show an actual amount of flavor change that is significantly smaller than the bound. For instance, the bound for the fiducial simulation only says that the n ee ≥ 0.
In summary, we show that the details of the neutrino distribution affect both the growth rate of the FFI and the final abundances of each neutrino flavor. Although we do not have a general way to map initial conditions to final neutrino abundances, we show that a small ELN crossing can produce minimal flavor change. Clearly, the amount of flavor transformation present in the final saturated state must be mapped out by a more comprehensive parameter sweep. We save this for future work.

V. CONCLUSIONS
We present the new open-source particle-in-cell (PIC) neutrino quantum kinetics code Emu. The PIC algorithm evolves the quantum states of individual computational particles by accumulating number density and flux on a background grid and then interpolating the net number density and flux from the background grid to evaluate the potential at each individual neutrino's location. The code is based on the AMReX framework and is performance portable to CPU and GPU hardware. We perform a series of test problems and a resolution study in the appendix.
We use Emu to perform simulations of the neutrino fast flavor instability (FFI) in one spatial dimension and two momentum (direction) dimensions. We simulate and analyze in detail a distribution of initially electron neutrinos and antineutrinos with number fluxes in opposite directions and a flux factor of 1/3 (Figures 1-3). We demonstrate the simultaneous evolution of the axial symmetrypreserving and symmetry-breaking modes with growth rates ( Figure 5) and wavelengths (Figure 3) that match predictions from linear stability theory (Figure 4). The symmetry-preserving mode has the faster growth rate, resulting in a high degree of axial symmetry during instability growth (Figure 7) that persists into the postsaturation phase (Figure 8). After the instability saturates, spatial fluctuations in neutrino flavor cascade from the length scale corresponding to the wavelength of the fastest growing axial symmetry mode to high wavenumbers ( Figure 6) following a simple analytic relationship (Equation 13). The final abundances of neutrinos and antineutrinos in this simulation approaches flavor equipartion at late times.
We then vary the initial distribution of antineutrinos from the fiducial simulation in a series of simulations. We show that that the antineutrino number density, flux magnitude, and flux direction all modify the growth rate of the instability (Figure 9). In addition, we show that varying antineutrino flux factor and number density result in different final abundances of each neutrino flavor after the instability saturates ( Figure 10). We furthermore show that a small ELN crossing leads to slower instability and little flavor change between the initial postsaturation states, and that the growth rate and flavor change increase with larger ELN crossings ( Figure 10). Although we do not have a way of mapping the initial distributions to final flavor abundances, we do demonstrate that in all of our simulations the same number of neutrinos and antineutrinos change flavor. That is, the final abundances of neutrino flavors imply the final abundances of antineutrino flavors (Equation 14). This relationship results in a bound on the amount of neutrino and antineutrino flavor change (Equation 15), though the simulation results never come close to this bound.
Our initial parameter study based on our highly symmetric fiducial simulation (Figure 9) demonstrated how a sparse exploration of the parameter space can lead to incomplete conclusions. Exploring more variations of the initial conditions (blue curve in Figure 10) revealed that the antineutrino flux factor affects the final abundances of neutrino flavor, even though the parameter sweep in Figure 9 showed all antineutrino flux factors resulting in flavor equilibration. We expect that it will be possible to build a mapping from initial distributions to final neutrino abundances, but this will require many more simulations. One major shortcoming of this work is that the simulations in this work were performed in a single spatial dimension. Although a dense covering of the parameter space using three-dimensional simulations is not currently feasible, it will be important to anchor the 1D simulations with select 3D ones. This, too, we leave to future work.
We also expect that the one-dimensional nature of the simulations in this work will significantly affect the process of kinematic decoherence and cascade of power to small scales due to the artificially reduced number of degrees of freedom. In addition, all of the neutrinos in these simulations were of a single energy. Although we include neutrino mass terms, the neutrino potential in these simulations is much larger than the vacuum potential, so all of these simulations are well outside the realm of collective oscillations. Many of these conclusions likely do not apply when the vacuum potential becomes significant. This is yet another avenue for future study.
With these prospects in mind, we note that the Emu implementation of modern and established numerical methods from the plasma physics community enables efficient numerical predictions with even modest computational resources. For scale, each of the simulations in the parameter survey in Figure 9  To test our implementation of the vacuum and matter potentials we follow [90] and perform a simple two-flavor one-zone test of neutrino oscillations in a regime where the neutrino potential is negligible and the matter density is tuned to put the neutrino oscillations in resonance. Specifically, we set up a 1 cm cube in a 1 × 1 × 1 domain filled with matter at a density of 1.35 × 10 9 g cm −3 and an artificial electron fraction of Y e = 1. The neutrino mass eigenstates have masses of m 1 = 8.60×10 −3 eV and m 2 = 0 eV, and the mixing angle is set to θ 12 = 33.82 degrees. We only initialize two particles, each of which represents a single neutrino and a single antineutrino at an  [36]) are overplotted as dotted black curves. The antineutrino values (green and red) undergo complete transformation, but since the neutrinos (blue and orange) are not in resonance they do not completely swap. The maximum difference between the analytic and simulated quantities is 3.7 × 10 −4 when using a timestep of 6.5 × 10 −13 s.

Bipolar Oscillations
To test our implementation of the isotropic components of the self-interaction potentials we follow [90] and perform an isotropic bipolar oscillation test. Once again, we use only one particle moving in the +z direction and one moving in the −z direction, since this exactly replicates the behavior of an isotropic distribution in the context of this test. We choose particle weights to ensure the neutrino and antineutrino number densities are n =n = 10(m 2 − m 1 ) 2 c 4 /(2 √ 2G F E) with energy E = 50 MeV. We use the same neutrino masses as in Section A 1, but we instead choose a mixing angle of θ 12 = 0.01 = 0.573 degrees and a background mass density of zero. This corresponds to the [99] and [90]. Note that unlike [90] we use an inverted mass hierarchy, so we start in an electron flavor eigenstate instead of a muon flavor eigenstate.

Homogeneous Two Beam Fast Flavor Instability
To test the anisotropic terms in the self-interaction potential, we perform a two-flavor two-beam fast flavor conversion test. The growth rate of the fast flavor instability in this simplified setup is easily calculated analytically as described by [48]. We use the same neutrino masses as in Section A 1, but choose a mixing angle of θ 12 = 10 −6 degrees. We once again use a single grid cell containing one particle moving in the +z direction representing only neutrinos, one moving in the −z direction representing only antineutrinos, and a vacuum matter background. The number of neutrinos each particle represents is set such that the total number density of each of neutrinos and antineutrinos is n = (m fastest growing instability given the neutrino parameters according to Equation 2.10 in Ref. [48]. The results are shown in Figure 13. The blue curve on top shows the evolution of the electron diagonal component of the density matrix for neutrinos (solid) and antineutrinos (dashed). As time progresses, the neutrinos mix into the the muon flavor state (red) and a quantum superposition of muon and electron flavor states (purple). Linear stability theory predicts an exponential increase in the off-diagonal component with a growth rate of Im(ω) = (m 2 2 − m 2 1 )c 4 /(2hhν) = 1123 s −1 . The fitted growth rates between t = 8.4 ms and 11.7 ms match the theoretical rate with a relative error of at most 3.2×10 −5 .

Inhomogeneous Fast Flavor Instability
Finally, we perform a simulation of the inhomogeneous two-beam two-flavor fast flavor instability to test all of the previous components of the code together with the spatial advection of neutrinos. We choose the same neutrino masses and mixing angle as in Section A 3. In order to get the wavelength of the fastest growing mode to match the size of our box (i.e. k = 2π/(1 cm)), we set each of the neutrino and antineutrino number densities to n = (m 2 2 − m 2 1 )c 4 /(2hν) +hck 2 √ 2G F = 4.89 × 10 32 cm −3 (A1) such that the fastest growing mode has a wavelength of 1 cm for a neutrino energy of 50 MeV (see [48] Equa-tion 2.7). To isolate the fastest growing mode, we also perturb the initial conditions using Re(f eµ ) = Re(f eµ ) = 10 −6 sin(kz 0 ), where z 0 is the initial position of the particle. The results are shown in Figure 14. As for the homogenous case in Section A 3, we see neutrinos (solid) and antineutrinos (dashed) mix from the nearly pure initial electron flavor state (blue) into the muon flavor state (red) and a quantum superposition of muon and electron flavor states (purple). We see the growth of the off-diagonal component is exponential, with a growth rate matching the prediction of linear stability theory of Im(ω) = (m 2 2 − m 2 1 )c 4 /(2hhν) + ck = 1.88 × 10 11 s −1 between t = 33 ps and 53 ps with a relative error of at most 4.3 × 10 −3 .

Appendix B: Convergence
The calculations we perform here are on a very small length and time scale (typically 64 cm domain size) compared to a supernova or neutron star merger. Thus, the quantity we are primarily interested in is the net amount of electron neutrinos that remain after the instability saturates, averaged over the entire domain and over time. We compute the spatial-temporal average of the neutrino and antineutrino density matrices N ab and N ab as where i labels a grid zone and j labels the simulation snapshot. j start is the first snapshot after t = 3 ns. This allows us to average only over times unaffected by the initial transitory phase of the instability growing and saturating. After saturation, the system is chaotic and thus impossible to demonstrate exact convergence. However, since the quantity significant for core-collapse and merger simulations is precisely N ab , this is the quantity we ensure is converged. There are three parameters dictating the fidelity of these simulations. The first is the domain size of the simulation. Larger domains can host modes with longer wavelengths, which even if they are not the fastest growing unstable modes, they may still be unstable or interact with shorter-wavelength modes after saturation. The second is the number of grid cells. The background grid discretizes the neutrino potentials, so finer grid cells are required to allow modes with shorter wavelengths to be represented in the background potential. Finally, there is the number of computational particles per cell. More particles means that a larger number of unique particle directions are represented on the grid and modes with finer angular structure can be expressed by the simulation. We use only one particle per direction per cell, since we have found that placing multiple particles at different locations in each cell but moving in the same direction does not result in a solution significantly different from the one calculated using only one particle per direction per cell. In the simulations discussed below, we initialize the neutrino distribution as described in Section III A to assess what simulation parameters lead to a numerically converged result.
The top left panel of Figure 15 shows the space-and time-averaged fraction of ν e in the post-saturation state under variations in the number of grid zones for a fixed domain size of 64 cm. Using eight or fewer grid zones causes the instability to never develop, leaving the final fraction of electron neutrinos equal to the initial conditions. As the number of grid zones increases beyond that, the saturated state appears to move increasingly toward an evenly mixed flavor distribution until 32 grid zones. Then the separation between the flavors increases slightly until 1024 grid zones, beyond which the results are mostly constant. Note that the numbers shown here are the values averaged over a simulation with only a duration of 5 × 10 −9 s, so this gap between electron and heavy lepton neutrinos is at least in part a result of not allowing the simulations to extend until momentum-space decoherence is complete. In the bottom left panel, we plot the fitted growth rate of the axially symmetric (blue) and asymmetric (red) modes along with the corresponding theoretical results from a linear dispersion analysis (gray horizontal lines). The growth rate of the symmetric mode appears to asymptote to the theoretical value. The growth rate of the asymmetric mode approaches the theoretical value until n z = 256 and then appears to slowly decrease. We expect that this is due in part to contamination from the underlying random perturbations, since the amplitude of the asymmetric mode is never large, and in part due to the limited resolution. The asymmetric mode does not significantly contribute to the postsaturation state, so we leave a more fine simulation of the asymmetric mode to future work.
A similar progression occurs for changes in the domain size while keeping the size of each grid cell constant (middle panels in Figure 15). If the domain has an extent of only 1cm (4 grid cells) only a weak instability occurs, since the grid cannot contain the fastest growing mode. The amount of transformation over the course of the duration of the simulation varies up until a domain size of 54 cm, after which it remains roughly constant. The growth rate of the symmetric mode is nearly correct at a 2 cm domain size, since it can nearly fit the fastest growing mode with a wavelength of 2.2 cm. There are once again small variations until a domain size of 64 cm, after which the growth rate remains relatively constant. The growth rate of the symmetric mode once again proves to be difficult, peaking at a domain size of 16 cm and then slowly decreasing with increasing domain size. We expect this to be the result of the same reasons as for the behavior in the n z convergence test above.
Finally, we test the robustness of the results while varying the number of particles per cell, which also allows the particles to probe a larger number of directions. The leftmost point in the right panels of Figure 15 represents a two-beam calculation. In each grid cell one particle is moving in the +z direction and one particle is moving in the −z direction. For the rest of the points, the x-axis describes the number of particle directions along the equatorial (x−y) plane. The remainder of the points particles are given directions according to Section II. 4,8,16,32,64,128, and 256 equatorial directions correspond to 6, 24, 92, 378, 1056, 6022, and 24088 particles per cell, respectively. The amount of post-saturation flavor transformation (top panel) does do not appear to be sensitive to the angular resolution of the simulation beyond 8 equatorial directions. The symmetric mode growth rate appears to vary little after 16 equatorial directions. Once again, the asymmetric mode growth rate is more difficult to resolve, but approaches the theoretical line with increasing angular resolution. It is interesting that such sparse sampling of the angular distribution can yield such accurate simulations of the FFI. This is in contrast to simulations of collective oscillations that require thousands of angular bins [37].
We also investigate our integration method's order of space-time convergence. We use the data set shown in the leftmost column of Figure 15 where we kept the extent of the domain constant while varying the number of grid zones, thus varying the grid resolution ∆z. We measure convergence by measuring the instability growth rate during the linear phase in the time interval [0.15, 0.25] ns, firmly centered in the linear regime as shown in Figure 5. We take as a reference point for computing the convergence error the highest resolution simulation shown in Figure 15 corresponding to n x = 2048 grid cells in the domain. We denote the grid resolution of our n x = 2048 simulation as ∆z 2048 and each other grid resolution where we measured convergence error as ∆z. We then denote the convergence error of the growth rates in the lowerresolution simulations as δω = |Im(ω − ω 2048 )| and compute this error for the same neutrino densities n ab and fluxes f (x) ab shown growing in Figure 5. In the top panel of Figure 16, we shift all logarithmic growth rate errors of each component of n ab by the minimum logarithmic error for that component across all resolutions. This permits us to clearly compare our achieved convergence order to representative lines (dashed) for first (red), second (magenta), third (cyan), and fourth (blue) order convergence rates. We see that towards high resolutions, density components n eµ (purple), n µτ (yellow), and n µµ (brown) all converge to second order or better. Exceptions are n eτ (green) and n τ τ (salmon), which converge between first and second order. We note that although we are using a fourth order Runge-Kutta method as our time integrator, the PIC deposition and interpolation shape functions we use reduce the achievable space-time convergence to second order, as we discuss in Section II B 1.
In the bottom panel of Figure 16 we next show the unscaled convergence error δω for the axially asymmetric flux f der, although f (x) eµ and f (x) eτ do not. We note nevertheless that at the fiducial resolution of n x = 1024 (the leftmost set of points in the lower panel of Figure 16), the size of the error δω corresponds to ≈ 10% error compared to the axially asymmetric growth rate controlling f These results justify our choice of fiducial simulation parameters. However, it should still be noted that these simulations are too short to probe the long-term behavior of the post-saturation nonlinear evolution.