Chaos in Inhomogeneous Neutrino Fast Flavor Instability

In dense neutrino gases, the neutrino-neutrino coherent forward scattering gives rise to a complex flavor oscillation phenomenon not fully incorporated in simulations of neutron star mergers (NSM) and core collapse supernovae (CCSNe). Moreover, it has been proposed to be chaotic, potentially limiting our ability to predict neutrino flavor transformations in simulations. To address this issue, we explore how small flavor perturbations evolve within a narrow centimeter-scale region inside a NSM and a toy neutrino distribution. Our findings reveal that paths in the flavor state space of solutions with similar initial conditions diverge exponentially, exhibiting chaos. This inherent chaos makes the microscopic scales of neutrino flavor transformations unpredictable. However, the domain-averaged neutrino density matrix remains relatively stable, with chaos minimally affecting it. This particular property suggests that domain-averaged quantities remain reliable despite the exponential amplification of errors.


I. INTRODUCTION
Neutrino flavor oscillation, proposed to explain the missing flux of neutrinos coming from the Sun, has ushered in a new era of neutrino physics beyond the standard model [1,2] and a wide phenomenology of neutrino flavor transformation in astrophysics.Neutrinos propagate in a quantum superposition of weak interaction eigenstates that is modified by the forward interaction of neutrinos with charged leptons via the Mikheyev-Smirnov-Wolfenstein (MSW) effect [3][4][5].In addition, coherent neutrino-neutrino forward scattering produces a not well understood non-linear flavor oscillation in dense neutrino gases [6][7][8][9][10].
Neutrinos play a pivotal role in NSM (see [11][12][13][14]) and CCSNe (see [15][16][17][18]).These represent the sole locations in the universe, post the big bang, where neutrinos are generated at densities high enough to undergo strong interactions with each other.CCSNe arise when the core of a massive star collapses following the depletion of its nuclear fuel.The infalling material halts abruptly as the core attains nuclear density, giving rise to a bounce and a shock wave that propagates throughout the star.Neutrinos originating in the deep regions of the star drive the explosion, transporting energy from the collapsed core to the material falling below the shock front.In the delayed CCSNe scenario, this mechanism propels the stalled shock wave through the star, inducing the explosion [19][20][21][22].Numerous simulations have yielded predictions for CCSNe explosion energies, neutrino luminosities, electromagnetic signals, and ejected material [23][24][25][26][27]. Several sources of uncertainty necessitate further investigation, including a consistent treatment of quantum neutrino transport [28,29], the equation of state for matter beyond nuclear densities (e.g., [30,31]), and the nuclear reaction rates of heavy and highly unstable elements [32][33][34][35].
NSM and CCSNe are prime candidates for the production of much of the heavy nuclei of the universe (e.g., [36,37]).Neutrinos emitted from the resulting hot accretion disk and ejected material can significantly affect the production of neutrons and protons via the reactions p + νe ↔ n + e + n + ν e ↔ p + e − that seed the conditions for slow and rapid neutron capture processes.Given that heavy lepton neutrinos do not partake in the creation and annihilation of protons and neutrons, a potential transformation of electron neutrinos into heavy flavors could significantly impact the ultimate production of heavy elements [38][39][40][41][42].
Despite its importance, the neutrino flavor transformation in dense neutrino environments is not fully understood.An ensemble of neutrinos that undergoes coherent forward scattering is an interacting quantum many-body system that could develop quantum entanglement.This can leave an imprint on the evolution histories of the neutrino flavor (see [62] for a recent review).Quantum many-body systems exponentially increase complexity as the number of particles increases.Neglecting helicity and pair coherence the Hilbert space of an ensemble of N interacting neutrinos and antineutrinos with n f number of flavors has a dimension of n N f .The exact solution for quantum many-body neutrino simulations have not involved more than 20 neutrinos [62][63][64][65].Approximations of the quantum many body problem are based on compact representations of the wave function like tensor networks [66][67][68] and other product state configurations [69].This allows hundreds of neutrinos to be simulated.Neglecting quantum entanglement, i.e., approximating expectation values of operator products as a product of the expectation values of individual operators, the neutrino and antineutrino ensemble can be described in the mean-field approximation by evolving N pure quantum states as n f × n f density matrices.This reduces the Hilbert space to a dimension of N n F .However, even under the mean field approximation, a global implementation of the full quantum neutrino transport in NSM and CCSNe is not computationally possible due to the non-linear complexity of the Hamiltonian and the small spatial and timescale of the flavor transformations.
The non-linear evolution of flavor leads to intricate and seemingly chaotic dynamics, posing challenges for achieving convergence in dynamical calculations (e.g., [99]).While chaos in the context of fast flavor instabilities remains unexplored, earlier investigations have hinted at chaos in bipolar oscillations within two-beam models (i.e., two momentum states) [100,101].Reference [100] presents a method for obtaining the complete spectrum of twelve Lyapunov exponents and covariant Lyapunov vectors across various initial conditions with different symmetries that either maintain or disrupt the periodicity of bipolar oscillations.However, the long-term evolution of small perturbations and their maximum magnitudes were not explored.Reference [101] demonstrates that small perturbations in neutrino flavor states yield a distinct flavor evolution history in a model incorporating fast axial symmetry-breaking modes.In such models, flavor oscillation modes either lack regularity or exhibit features resembling periodic behavior for a brief duration.The introduction of a small perturbation does not guarantee that the system will evolve in a manner where its evolution paths closely align.
Future observations of CCSNe and NSM neutrinos will provide crucial scientific insights into the collapsing stellar core (e.g., [102]) and nucleosynthesis conditions.The exponential amplification of errors in simulations, due to chaos, may limit our ability to predict the evolution of neutrino flavor within CCSNe and NSM, hindering our capacity to compare theoretical predictions with observational data.To shed light on the problem, we extend the study of chaos to a more realistic distribution of neutrinos that exhibit fast flavor instabilities.We conduct two simulations of dense neutrino gases with distinct angular distributions, one of which involves a narrow domain a few centimeters wide, capturing a snapshot of a neutron star merger (NSM).The resolution scale is defined by a domain of 1×1×64 cm divided into 1×1×1024 cells, each accommodating 24, 088 momentum states.This setup allows for flavor anisotropies and inhomogeneities.The primary goal is to characterize the evolution of small perturbations in neutrino flavor states over time, specifically quantifying the sensitivity of the results to the initial conditions.This endeavor will shed light on how chaos affects our capability to predict the evolution of neutrino flavor in dense astrophysical neutrino gases within simulations.
In Section II, we briefly overview the standard theory of chaos in dynamical systems.Neutrino quantum kinetics is reviewed in Section III, while Section IV outlines our numerical approach to simulating neutrino flavor transformation.Evidence of chaos is presented in two extensive one-dimensional simulations in Section V. We further explore the unpredictable nature of chaos and its impact on both large and small scales of neutrino flavor transformation in the same Section V. Concluding remarks are provided in Section VI.

II. CHAOS
The discovery by E. Lorenz of an atmospheric convection model exhibiting irregular and unpredictable behavior formally marked the beginning of the study of chaotic systems [103].Chaotic systems are non-linear dynamical systems characterized by erratic and irregular complex behavior.Although these systems are fundamentally deterministic -meaning precise knowledge of initial conditions allows for the prediction of future behavior -even a minute variation in these conditions, such as measurement uncertainty, results in an entirely different prediction.This contrasts with non-chaotic systems, where the approximate present determines the fu-ture approximately.In chaotic systems, the approximate present leads to an entirely different future (see [104][105][106][107][108] for modern expositions of chaotic systems).Since Lorenz's findings, many dynamical systems have been discovered to exhibit chaos, revealing how chaos and complexity can emerge from seemingly simple dynamic systems such as the logistic equation [109], a driven damped pendulum [110], or a double pendulum [111].Chaos, rather than being an isolated behavior of specific dynamical systems, is a universal property of complexity in all of them [108].
Various methods can be employed to test whether a dynamical system exhibits chaos.One commonly used approach involves measuring how sensitively results depend on initial conditions through the Lyapunov exponent.Consider a general dynamical system described by state space coordinates r i following the time-evolution equation with an initial condition ⃗ r t0 .The solution ⃗ r t traces a path in the state space that we assume is bounded to a certain volume.Consider a nearby initial state ⃗ r t0 + ⃗ δ t0 where the magnitude of the perturbation | ⃗ δ t0 | is extremely small.If the perturbation satisfies λ represents a Lyapunov exponent of the dynamical system and ⃗ δ t a Lyapunov vector.The Lyapunov exponent characterizes the average rate of exponential divergence or convergence of nearby trajectories in state space [112].
A positive Lyapunov exponent indicates exponential divergence and chaos, while a negative exponent represents exponential convergence.A zero Lyapunov exponent implies regions in state space where ⃗ δ t0 remains constant over time.
The exponential divergence of two nearly identical states suggests that systems with subtle, hard-to-notice distinctions in the initial conditions will soon start to behave differently, and the ability to make forecasts will be rapidly reduced.The magnitude of the Lyapunov exponent indicates the rate at which the system's behavior becomes unpredictable.
Over short time intervals, perturbations may exhibit intricate behavior that is not clearly identified as exponential divergence or convergence.Chaos is identified by the average exponential increase in perturbation size, even in the presence of rapid nonlinear fluctuations on a small time scale.In bounded dynamical systems, perturbations are expected to cease growing exponentially once they reach certain magnitudes, as state space coordinates cannot exceed values within a specific domain.
In an n + 1-dimensional state space, there will be n + 1 Lyapunov exponents.A practical representation is as a small n-sphere of initial conditions ⃗ r t0 + ⃗ δ k t0 centered at a given point ⃗ r t0 .k = 1, 2, ... n + 1 represent each state space dimension of the dynamical system.As time progresses, this n-sphere evolves into an n-ellipsoid, with each dimension stretching and contracting at different rates according to its Lyapunov exponent λ k as If the perturbations evolve for an extended period, the largest Lyapunov exponent will dominate the shape of the n-sphere, and the perturbations will follow the direction of maximum divergence [104,108].

III. NEUTRINO QUANTUM KINETICS
A proper treatment of a flavor neutrino field in a dense environment involves a seven-dimensional distribution function for neutrinos f ab (t, x, p) and antineutrinos fab (t, x, p) (n F × n F matrices, where n F is the number of neutrino and antineutrino flavors) whose elements are the expectation values of bilinear creation and annihilation operators < a † a a b > and < b † a b b > in the mean field approach.The diagonal terms (a = b) of f ab (t, x, p) and fab (t, x, p) represent the neutrino and antineutrino occupation numbers for each flavor (a, b ∈ {e, µ, τ }), and the off-diagonal term (a ̸ = b) represents the flavor correlation [8,10].The time evolution of the neutrino distribution function is governed by the quantum kinetic equation (QKE) Barred quantities ( fab , Cab and H) need to be substituted to arrive at the equation for antineutrino distributions.The matrix C ab ( Cab ) is the collision term which accounts for the non-forward scattering of neutrinos(antineutrinos).The first term on the left accounts for the explicit time dependence of the distribution function and the second term accounts for the drift that is caused by the particles that are streaming freely.We ignore other external forces like gravity and hypothetical electromagnetic interactions.The Hamiltonian H ab accounts for the flavor transformation due to mixing between mass and flavor eigenstates H vacuum ab , the coherent forward scattering of neutrinos with background matter (other than neutrinos) H matter ab , and the coherent forward scattering of neutrinos with others neutrinos and antineutrinos H ν−ν ab .These are combined as .449× 10 −3 eV 2 under the neutrino normal mass ordering as in reference [113], although the calculations are in a regime where these terms are of negligible import.p is the neutrino momentum.n a (n a ) denote the scalar number density and f a ( fa ) denote the number flux density for leptons(anti-leptons) of flavor a.The Hermitian neutrino number density n ab and number density flux f ab matrices come from the neutrino distribution function as follows: The antineutrino equations are given by substituting barred quantities as In this work we suppress the matter term in the Hamiltonian to focus on the flavor transformation produced by the neutrino-neutrino coherent forward scattering mechanism.

IV. METHODS
We employ mean-field neutrino quantum kinetics simulation in a 1+2D setup, assuming homogeneity in two directions and a single neutrino energy.We use the multidimensional neutrino quantum kinetics code Emu, which we describe in Section IV A. The pair of initial conditions for the neutrino systems under study is detailed in Section IV B.

A. Emu
Emu [78] is a particle-in-cell code designed for simulating neutrino flavor transformations.It solves the QKE without considering collisions under the mean-field approximation.Emu utilizes computational particles to discretize the neutrino field into distinct packages and tracks the evolution of each computational particle's position on a background grid.This approach efficiently accounts for the advection term in the Liouville operator v • ∇ x by moving the computational particles in accordance with the equation where r is the position vector of the particle and p is the momentum vector of the particle.This implies that there is no momentum exchange due to incoherent scattering, creation and annihilation of neutrinos and antineutrinos and energy exchange C ab = Cab = 0 so that each computational particle in EMU satisfies where N and N are the number of physical neutrinos and antineutrinos that a computational particle represents.Each computational particle carries two quantum states defined by the Hermitian density matrices ρ and ρ, representing the physical neutrinos and antineutrinos carried by the particles.The flavor quantum state evolves in time according to the Schrödinger equation To compute the Hamiltonian based on a discrete approximation of the neutrino field, EMU employs a deposition and interpolation algorithm to estimate the number density of neutrinos and antineutrinos (n ab and nab ) and the number density flux (f ab and fab ) at every position of the particle.During each simulation time step, EMU models each computational particle with an extended shape centered on the particle position, with an extent comparable to the size of the grid cell.This is used to combine the quantum states of the computational particles that overlap over a single cell based on the parameterization of a second-order shape function.This combination attributes a distribution of neutrinos and antineutrinos from each particle to each cell of the grid.Subsequently, the neutrino distribution moments are interpolated from the grid cell distribution to the location of each computational particle, using the same second-order shape function.The quantum state is then integrated using a fourth-order Runge-Kutta method in an efficient and scalable manner based on the AMReX framework, ensuring performance portability to both CPU and GPU hardware.

B. Initial conditions
We study the presence of chaos in two separate simulations.In the first (Section IV B 1) we look at a wellunderstood toy model, and in the second (Section IV B 2) we extract the neutrino distribution from a multidimensional NSM simulation.In both cases, we construct a simulation domain that is one-dimensional in space (i.e., assuming homogeneity in the x and y directions) using a one-dimensional array of 1 × 1 × 1024 cells embedded in a domain of 1 × 1 × 64 cm with periodic boundary conditions.The final flavor content resulting from the flavor instabilities in one-dimensional simulations provides a reasonably accurate approximation compared to computationally expensive three-dimensional simulations [79].We simulate 24, 088 computational particles initialized at the center of each cell, with isotropically distributed directions ensuring that each computational particle represents the same solid angle.This significantly larger particle count per cell is chosen as we observe that the properties of chaos are notably more sensitive to numerical errors than the global distribution averages investigated in previous works [78,79].

Fiducial simulation
In this simulation, the neutrino and antineutrino distributions have fluxes in opposite directions with a flux factor of 1/3 and number densities of 4.89 × 10 32 cm −3 .The initial angular distributions satisfy where n να (n να ) is the neutrino(antineutrino) number density of flavor α. dΩ is a differential solid angle in momentum space centered on the momentum direction p.The direction z defines the direction from which θ is defined zero.Figure 1 shows the initial angular distribution in the xz plane.The linear angular distribution exhibits a strong lepton-number crossing in the xy plane.Previous work shows that this distribution exhibit fast flavor instabilities that match the unstable modes predicted by linear stability analysis [78].The input parameters of this simulation are publicly available in the EMU code [114].
To seed the instability and set a small initial amplitude for fast unstable flavor modes, the quantum state ẑ x ν e νe FIG. 1.Initial angular distribution for neutrinos and antineutrinos in the one-dimensional Fiducial simulation.A rotation of the curves around the ẑ direction (ẑ axial symmetry) reproduces the three-dimensional neutrino angular distribution.The neutrino and antineutrino number density emitted per solid angle is proportional to the distance of the curves from the origin.
of each computational particle is initialized in nearly the pure electron flavor state, with a small random perturbation on the non-diagonal electron-muon and electron-tau components: Here, Q is a random number generated between −1 and 1 each time it appears, and α is the strength of the random perturbation, which we set as 10 −6 .ϵ µ and ϵ τ are computed once the random non-diagonal components are generated in order to satisfy the conservation of the length of the flavor polarization vector and the unit trace of the density matrix.

NSM snapshot simulation
We evolve the quantum flavor states of neutrinos in a NSM snapshot after 5 ms post-merger from the classical global general relativistic two-moment radiation hydrodynamics simulation in [115], at a location approximately 40 • from the accretion disk plane and at 30 km from the compact object center as in [96].The total number densities and fluxes of neutrinos and antineutrinos is shown in Table I.To maintain consistency with the assumptions made in the original NSM simulation, each computational particle carries a number of neutrinos and antineutrinos that satisfy the angular distribution determined by the maximum entropy closure [116].This maximizes the angular entropy of the energy-integrated distribution while conforming to the total neutrino and antineutrino number densities and fluxes outlined in Table I.Specifically, the initial distribution function is given by where n aa is the number density of flavor a, and θ is the angle between the momentum direction and the flux direction.Barred quantities ( faa and naa ) need to be substituted for antineutrino distributions.The factor Z is numerically calculated to yield the expected total number density flux outlined in Table I.This distribution exhibits a angular directions of equal fluxes of neutrinos and antineutrinos [96], indicating its unstable to fast flavor transformations (e.g., [84]).

V. RESULTS
We explore how small perturbations in the flavor state of neutrinos and antineutrinos evolve over time using Lyapunov exponents (see Section II), on both macroscopic and microscopic levels of the neutrino flavor transformation.The macroscopic scales are represented by the domain-averaged density matrix, while the microscopic scales are computational particles that represent the smallest spatial resolution of the neutrino distributions in the simulations.
We first simulate the neutrino flavor transformation in a domain a few centimeters wide, located above the accretion disk 5 ms after the merger of two neutron stars [115] (see Section IV B 2 for a description of the simulation setup).This is a noteworthy area for the generation of r-process nuclei [52,53] and the neutrinos are unstable to the fast flavor instability, making this region an optimal testing ground for researching chaotic properties in neutrino flavor evolution.We also simulate a well-studied artificial neutrino distribution in which one third of the total number of neutrinos and antineutrinos travel in opposite directions (see Section IV B 1 for a description of the simulation setup).This is an extreme and unique initial condition that is useful for understanding of the patterns of chaos.Although, the three-dimensional neutrino flavor transformation for these systems have been reported in [79,96], we conduct high-resolution one-dimensional simulations to precisely compute the Lyapunov exponents.
Essential data have been archived and are publicly accessible at [117], with additional data available upon request.

A. Overall Dynamics
To start, we need to first examine the impact of the fast flavor instability on the flavor distribution before delving into the additional effects introduced by chaos.This preliminary analysis will aid in interpreting the results presented in the subsequent sections.Figure 2 shows the neutrino density matrix averaged across the spatial domain, denoted as ⟨ ρ ⟩.For flavor components i, j ∈ { e, µ, τ }, the domain-averaged density matrix is defined as a sum over the computational particles according to where Here, N par represents the number of computational particles in the simulation, N k is the number of physical neutrinos represented by the computational particle k, and ρ k ij denotes the neutrino density matrix of particle k.
As previously reported in [78,79,96], the evolution of neutrino flavor exhibits three distinct phases: linear growth, saturation, and decoherence of flavor oscillation modes.During the linear growth phase, observed between 0.10 and 0.25 ns in both panels of Figure 2, the off-diagonal components of the domain-averaged density matrix are predominantly influenced by highly unstable flavor oscillation modes that exhibit exponential growth as ⟨ ρ ⟩ ij = Ae −iωt with Im(ω) ̸ = 0 and A a small initial amplitude.The imaginary frequency of the dominant unstable modes can be read from the slopes of the off-diagonal curves during the linear growth phase.In After the linear growth phase, the flavor oscillation modes that begin the simulation with a negligible amplitude become noticeable, provoking saturation of the flavor oscillation modes.This saturation phase occurs between 0.25 and 0.30 ns in both simulations.Significant flavor conversion does not take place until the conclusion of the linear growth phase and the initiation of the saturation phase.After the saturation phase, all the flavor oscillation modes start evolving in an incoherent mixture.This is the flavor decoherence phase.Coherent modes that grew in the linear growth phase are disrupted, and the phase of the oscillation modes is randomized, although long-lived modes can exist after the saturation phase [118].In the decoherence phase, neutrinos and antineutrinos start evolving in a complex non-linear relationship, archiving huge amounts of flavor conversion.After about 0.3 ns, the domain-averaged neutrino density matrix in the NSM simulation reaches an incomplete flavor mixing state as an equilibrium point and fluctuates slightly around that state.The domainaveraged antineutrino density matrix follows the same pattern.The dynamics of the Fiducial simulation are remarkably similar to that of the NSM simulation, except that the equilibrium distribution of flavors is very close to an even mixture.

B. Chaos
To quantify the evolution of the small perturbation ⃗ δ t over time, we introduce a mathematical vector henceforth referred to as the flavor vector ⃗ r, with components given by Here, k ∈ {0, N par ) denotes the particle index, l ∈ [0, 7] represents the Gell-Mann vector component of a single particle's neutrino flavor vector, and l ∈ [8,15] designates the Gell-Mann component for a particle's antineutrino flavor vector.P k l = 1 2 Tr(ρ k G l ) is the neutrino flavor polarization vector for the computational particle k, where G l refers to the Gell-Mann matrices (with the eight Gell-Mann matrices being repeated for l ∈ [8,15]).Given that l ∈ [0, 15], the flavor vector has 16N par components.The flavor vector completely determines the flavor content of our simulations and its evolution can be thought of as tracing a path through a 16N par -dimensional state space.Given that the length of the flavor polarization vector of each particle is a constant of motion, the trajectories of the flavor vector in the state space is bounded to the surface of a 16N par -dimensional sphere of radius the length of the flavor vector.
The evolution over time of small perturbations can be obtained by simulating the paths in the state space of two flavor vectors with initial conditions ⃗ r bas (t 0 ) = ⃗ r t0 and ⃗ r per (t 0 ) = ⃗ r t0 + ⃗ δ t0 .Here, ⃗ r t0 is the initial flavor vector of the NSM and Fiducial simulations described in Sections IV B 2 and IV B 1, respectively.The perturbation ⃗ δ t can be obtained by computing The first simulation, described by ⃗ r bas serves as the baseline simulation for comparison.The other, described by ⃗ r per , is randomly perturbed from the baseline simulation.Specifically, we randomly choose the initial magnitude of the perturbation so that This choice of perturbation is sufficiently large to avoid numerical errors (see Appendix A).The mathematical maximum magnitude of the perturbation at any point in time is 2|⃗ r t | (i.e., when ⃗ r per is oriented opposite to ⃗ r bas ).It is important to note that the "chaos" perturbation is independent of the perturbation used to seed the instability (see Section IV), which is present in both the baseline and perturbed simulations.In the following discussion, the perturbations referred to are the "chaos" perturbations and not the initial instability seed.
If the neutrino flavor evolution is chaotic, the time evolution of small perturbations follows an approximately exponential trend (see Equation 1) with λ > 0. To test whether perturbations exhibit this exponential behavior, we performed two simulations in which we apply perturbations at t 0 = 0 and t 0 = 2.65 ns in both the NSM and Fiducial simulations.The first case represents a perturbation applied before the linear growth phase, and the second a perturbation applied in the decoherence phase.
The blue curve on the left panel of Figure 3 shows the evolution of a perturbation applied at t 0 = 0 in the NSM simulation.During the linear growth phase, between 0.10 and 0.25 ns, the perturbation follows an exponential trend characterized by λ ≈ 59 ns −1 .This trend is mainly driven by fast flavor unstable modes, and does not represent a signature of chaos.Additionally, as previously exposed, there is no significant flavor conversion at this time.Towards the end of the linear growth phase, the perturbation's exponential trend abruptly transitions to a slower rate (λ ≈ 2.6 ns −1 ).At this point, as observed in Figure 2, unstable modes have vanished and the domain-averaged density matrix evolves into a complex, incoherent mixture of flavor oscillation modes that has attained an equilibrium distribution, closely resembling flavor equipartition.This suggests that the exponential growth of the perturbation in decoherence phase is not driven by fast flavor unstable modes and is a signature of chaotic flavor evolution in the flavor vector.The Fiducial simulation (orange) shows a very similar trend for the same reasons.The growth rates before and after saturation in this case are approximately 58 ns −1 and 2.6 ns −1 , respectively.
To isolate the chaotic behavior from the fast flavor instability, we perform another simulation where we apply the perturbation at t 0 = 2.65 ns during the decoherence phase (right panel of Figure 3).In the NSM simulation (blue), the perturbation amplitude exhibits an approximately exponential growth with λ = 0.44 ns −1 , and in the Fiducial simulation (orange), the growth rate is λ = 0.32 ns −1 .This observation highlights that the paths of similar flavor vectors in the state space diverge exponentially, indicating chaotic behavior.The Lyapunov exponent is one order of magnitude smaller than the growth rate of the dominant flavor unstable modes present in the linear growth phase.
It is also curious that the growth rates for perturbations applied before the linear growth phase are so different from those applied in the decoherence phase.The Lyapunov exponent depends on the position in the state space of both baseline and perturbed flavor vectors (i.e., the direction of the 16N par -dimensional perturbation vector ⃗ δ t ), so as the system evolves λ can change, and a single simulation is not able to fully explore such a large parameter space.To illustrate this, we run a third set of simulations in which we periodically normalize the perturbation maintaining its direction every 0.9 ns.The upper panel of Figure 4 shows the time evolution of the perturbation after each normalization, while the lower panel provides an approximation of the Lyapunov exponents calculated over each interval between normalizations.Given the significant differences in the initial conditions, it is remarkable how similar chaos manifests in the NSM and Fiducial simulations.The Lyapunov exponents emerging from the pre-saturation perturbations are larger because the fast flavor instability causes the perturbation to grow in a very specific direction that is closer to directions that exhibit the strongest chaos.Randomly perturbing after the instability saturates simply sets the stage with a random perturbation in a less chaotic direction.In bipolar simulations involving only two beams it is possible to obtain the full spectrum of Lyapunov exponents [100], but that is infeasible with these large-scale simulations.However, we do find that random perturbations produce consistent Lyapunov exponents across multiple runs, leading us to believe that the numbers reported for the post-saturation perturbations represent a lower bound for the chaoticity of the system.
Although the exponential trend of the perturbations is predominantly driven by positive Lyapunov exponents (λ > 0), the upper panel reveals the existence of directions with negative Lyapunov exponents (λ < 0).One such direction is marked by each a blue and orange star.Stable directions with negative Lyapunov exponents are swiftly suppressed by divergent directions.In other words, perturbations that initiate evolution in a stable direction will eventually transition to a direction characterized by a positive Lyapunov exponent.This is another illustration of the dependence of the Lyapunov exponent on the state itself and the direction of the perturbation.
Negative Lyapunov exponents arise from the conservation of the state space volume of the canonical variables, as stated in the Liouville theorem.The QKE can be transformed into a conservative classical Hamiltonian system with the canonical coordinates and momenta given by any function of the coordinates used to describe the flavor polarization vector (see Section II A of [100]).This implies that for each direction in the state space with exponential divergence, there will be a direction of exponential convergence at the same rate, ensuring the conservation of the state space volume of the canonical variables, even if its shape is deformed.The exponential divergence of small perturbations also has consequences for the integrability of the QKE.In Hamiltonian systems, for each conserved quantity, there exists one direction in the state space with a zero Lyapunov exponent.Since In the left panel, the perturbation is applied at t0 = 0 ns (before the linear growth phase).Between 0.10 and 0.25 ns, the perturbation growth is exponentially driven mainly by the fast flavor instability.In the right panel, to avoid the presence of unstable fast flavor modes, the perturbation is applied in the decoherence phase at t0 = 2.65 ns.The magnitude of the perturbation grows, following an approximate exponential trend.This demonstrates that paths of similar flavor vectors in the state space diverge exponentially, illustrating the chaotic evolution of neutrino flavor.As the simulation concludes, the magnitude of the perturbations approaches a constant value larger than the magnitude of the flavor vector |⃗ rt| (dotted lines).
integrable systems have conserved quantities as degrees of freedom, the Lyapunov exponent is zero for every direction in the state space.In other words, the shape of the state space volume of the canonical variables remains unchanged.Even though our state space does not correspond to the state space of the canonical variables, our finding that close flavor vector paths in state space diverge exponentially suggests that the QKE is a nonintegrable system.
Towards the end of the simulation, the magnitude of the perturbation stabilizes, reaching a constant value (t ≳ 4 ns in the left panel and t ≳ 60 ns in the right panel).In both cases, the perturbation's magnitude asymptotically settles to a value less than the maximum possible value of 2|⃗ r t | but greater than |⃗ r t |.If we interpret ⃗ δ t as the uncertainty in the flavor vector ⃗ r t ± ⃗ δ t , the relative error in ⃗ r t can be up to 106% for the NSM simulation, and 142% for the Fiducial simulation.This implies that even small uncertainties are exponentially amplified making the flavor vector unpredictable on a time scale of the inverse of the Lyapunov exponent.The Fiducial simulation maximum perturbation amplitude is considerably larger than in the NSM simulation because the Fiducial simulation experiences a higher amount of flavor conversion in all directions.All particles in the Fiducial simulation change flavor from a pure electron state to a flavor-mixing equilibrium.

C. Individual Particles
While the domain-averaged density matrices tend to a flavor mixing equilibrium, this stationary state comprises numerous particles whose flavor states undergo random and incoherent fluctuations at small scales.The fine details of these fluctuations, observable at the level of single computational particles, remain highly chaotic.This can be seen in figure 5, which shows the neutrino flavor transformation of a single computational particle in the EMU code for the NSM simulation.The upper panel shows the magnitude of the off-diagonal components of the neutrino density matrix, while the lower panel displays the diagonal components.Similar to Figure 2, the same linear growth, saturation, and decoherence phases are discernable.Flavor transformations become evident toward the end of the linear growth phase, around 0.25 ns.The amount of flavor conversion that a computational particle experiences is related to the direction of its momentum.The flavor transformation follows a complex non-linear trend with no discernible pattern.
The Lyapunov exponent extracted from individual particles is similar to those obtained from the full state vector in Section V B. The blue (NSM) and orange (Fiducial) curves in the upper panel of Figure 6 show the difference between the flavor state of a single particle extracted from the baseline simulations and the simulations perturbed at 2.65 ns.To avoid focusing on a single flavor, The blue curve corresponds to the NSM simulation, while the orange curve represents the Fiducial simulation.The lower panel shows the Lyapunov exponent calculated over each interval between normalizations.The Lyapunov exponent of a perturbation in the state space varies based on the direction of the perturbation.The presence of exponential divergence (λ > 0) and convergence (λ < 0 e.g.below the blue and orange stars) is an expected result of the Liouville theorem of conservation of the state space volume of canonical variables.
we sum over flavors and plot This can be interpreted as the combined uncertainty in the average electron, muon, and tau flavors generated by the perturbation in the flavor vector, whose maximum value is two.This difference follows a similar exponential trend (λ ≈ 0.41 ns −1 for NSM and 0.32 ns −1 for Fiducial) as the full flavor vector (right panel of Figure 3).Toward the end of the simulation, this difference ceases exponential growth and instead shows irregular and large-amplitude fluctuations.The maximum combined uncertainty reaches values of up to 0.17 for the NSM particle and 1.92 for the Fiducial particle.By tracking particles propagating in other directions (not shown), we note that computational particles experiencing high-amplitude flavor oscillation modes reach values of Tr|∆ρ| close to two (similar to the Fiducial particle), while particles undergoing low-amplitude flavor oscillation modes exhibit small values of Tr|∆ρ| (like the NSM particle).

D. Domain-Averaged Quantities
The flavor vector in Equation ( 4) fully characterizes a many-body flavor quantum state (which by construction never builds multi-particle entanglement).However, in astrophysical applications, such as CCSNe and NSM, this level of detail is unnecessary and computationally infeasible.A more relevant quantity is the domain-averaged density matrix over a spatial domain greater than the flavor oscillation length scale.The quantification of the FIG. 5.The upper panel shows the magnitude of the nondiagonal components of the neutrino density matrix for a single computational particle in EMU, obtained from the NSM snapshot simulation.The lower panel shows the diagonal components of the neutrinos density matrix for the same computational particle.The flavor transformation becomes evident at the end of the linear growth phase and evolves with large amplitude, in contrast to the domain-averaged quantities in Figure 2.
impact of chaos on these quantities is important for a reliable implementation of the neutrino flavor transformation in CCSNe and NSM simulations, where even a small numerical error could propagate exponentially.
To investigate how small perturbations propagate on macroscopic scales, we compute the domain-averaged density matrices defined in Equation ( 3) for the baseline simulations and the simulations perturbed at 2.65 ns used in the right panel of Figure 3.The lower panel of Figure 6 shows the time evolution of the trace of the difference between the baseline and perturbed neutrino density matrices averaged in the domain, defined as in Equation ( 5) but using the domain-averaged density matrix ⟨ρ⟩.In the NSM simulation the maximum value it reaches represents 0.07% of the maximum possible value of two, while in the fiducial simulation it represents 1.3% of the maximum.This is considerably smaller than the relative error of the flavor vector due to the same perturbation, implying that small perturbations do not escalate  5) of single computational particles extracted from the baseline and postsaturation perturbed simulations.The lower panel shows the same quantity, but using the domain-averaged density matrix ⟨ρ⟩.The blue and orange curves denote the NSM and Fiducial simulations, respectively.The black dotted lines show the maximum possible value.Both panels exhibit a comparable exponential trend, akin to the perturbations observed in the flavor vectors (see the right panel of Figure 3), but late-time uncertainties in the domain-averaged quantities are much smaller than in individual particles.
to large magnitudes in domain-averaged quantities.The diagonal components of the domain-averaged density matrix have a high degree of predictability compared to the flavor vector even amidst chaotic flavor evolution.Due to this distinctive feature, (i.e., small initial perturbations do not reach large magnitudes) it is reasonable to use the domain-averaged density matrix as a key variable in simulations of neutrino flavor transformation in CCSNe and NSM.This is also encouraging for thermodynamic theories of flavor transformation (e.g., [119]).

VI. CONCLUSIONS
To inform simulations of neutrino flavor transformation in CCSNe and NSM, we delve into the chaotic na-ture of neutrino flavor transformations in two distinct one-dimensional dense neutrino gases.In the first simulation (Section IV B 2), we extract the neutrino distribution from a region above the accretion disk of a multidimensional NSM simulation, and in the second simulation (Section IV B 1), we examine a well-understood toy model.Both distributions undergo significant fast flavor instabilities (see Section V A).
Using Lyapunov exponents, we study the stability of paths in the state space of similar flavor vectors (Equation 4) with similar initial conditions.These paths diverge exponentially (see Section V B) showing that the neutrino flavor transformation is chaotic.Interpreting the perturbation as uncertainty results in a final relative error of over 100 percent in the position of the flavor vector in the state space that grows with a Lyapunov exponent with a lower limit of 0.44 ns −1 (NSM) or 0.32 ns −1 (Fiducial), as shown in Figure 3. Since the flavor vector completely characterizes the flavor content in the simulations, we conclude that the long-term evolution of neutrino and antineutrino flavors at this level of detail is unpredictable.
We analyze how the Lyapunov exponents depend on the direction in the state space of the perturbation (see Figure 4).In both simulations, we find that the time evolution of perturbations is mainly driven by directions of positive Lyapunov exponents.A perturbation applied in a stable direction (λ < 0) will eventually move to an unstable direction in state space (λ > 0).
We identify stable directions in state space in which the perturbations briefly converge exponentially below the stars in the upper panel of Figure 4. Since the QKE can be transformed into a classical Hamiltonian system (as in Equation 2 of [100]) several Hamiltonian mechanic theorems can be applied.The Liouville theorem claims that the state space volume of the canonical coordinates is a constant of motion, i.e., the spectrum of the Lyapunov exponents is symmetric (λ 1 , λ 2 , λ 3 , ... , −λ 3 , −λ 2 , −λ 1 ) , meaning that for each direction in the state space with exponential divergence, there exists another direction of exponential convergence at the same rate, so the volume of the state space is conserved even though its shape is deformed.One Lyapunov exponent is zero for each conserved quantity, and all of them are zero for a stable or integrable system where there exist conserved quantities as degrees of freedom.The presence of a direction of exponential divergence and convergence in our simulations suggests that the QKE forms a non-integrable system.
The unpredictable nature of chaos manifests differently on the macroscopic and microscopic scales of neutrino flavor transformation.The macroscopic scale is represented by the domain-averaged density matrix, while the microscopic scale is the quantum state of single computational particles that represents the smallest spatial resolution of flavor in the simulations.On single computational particles, small perturbations grow exponentially (see Sec-tion V C), reaching magnitudes of order unity for particles moving in angular directions with a large amount of flavor conversion.This implies a relative error of 100% on the quantum states of neutrinos and antineutrinos, meaning that slightly different initial flavor vectors reproduce a different flavor evolution in single computational particles.Individual computational particles destroy information at a rate of λ log 2 e ≈ 0.6 bits/ns, implying that initial variables specified with 64-bit floating-point precision (of which only 53 bits are used to store significant digits) cannot be accurately predicted after approximately 100 ns.
In the domain-averaged density matrix, small perturbations in the flavor vector grow exponentially reaching small magnitudes compared to the maximum possible value.In both simulations, this produces a combined maximum uncertainty in the diagonal components of the domain-averaged density matrix of less than 1.3% (see Section V D).This suggests that CCSNe and NSM simulations could safely rely on domain-averaged quantities even if microscopic details are unpredictable.
In the following sections, we examine the impact of the domain size, the number of cells, the number of particles, and the magnitude of the perturbation on the Lyapunov exponent.This establishes the parameters for which the simulation accurately reproduces the Lyapunov exponent, minimizing numerical errors.To approximate the Lyapunov exponent in Equation ( 1

Number of particles
To determine whether the Lyapunov exponent is affected by the number of particles, we conducted five simulations wherein neutrinos are emitted in a roughly isotropic distribution from the center of each cell.We considered 92, 378, 6, 022, 24, 088, and 54, 202 particles per cell for each case.The simulation domain is 1×1×64 cm, divided into 1024 cells in the ẑ direction.
Figure 7 shows the Lyapunov exponents for the five simulations in this test.As the number of particles (and hence the number of directions) increases, the Lyapunov exponent asymptotes to 0.6 ns −1 .Around 24, 088 directions, convergence is achieved.This represents a highprecision simulation compared to the 92 particles needed to achieve convergence in the long-term domain-averaged density matrix [78].

Perturbation magnitude
To evaluate how the initial magnitude of the perturbation affects the Lyapunov exponent, we perform three sets of five simulations with initial perturbation magnitudes | ⃗ δ t0 |/|⃗ r t0 | on the order of 10 −14 , 10 −12 and 10 −10 .The simulation domain is 1 × 1 × 64 cm divided by 1024 cells.Each cell is initialized with 24, 088 particles at the center of the cell and the perturbation is applied at t = 2.65 ns.We terminate the simulation at 20 ns or until | ⃗ δ t |/|⃗ r t | = 10 −6 .
Figure 8 show the Lyapunov exponent (blue, orange and green) for three sets of five simulations with distinct initial perturbation magnitudes (10 −14 ,10 −12 , and 10 −10 ).The blue, orange, and green shaded areas represent a spread of 2σ in the Lyapunov exponents of the same color.Larger perturbation magnitudes result in less dispersion in the Lyapunov exponent without significantly changing the mean (0.61 ns −1 ), confirming their independence from the perturbation magnitude.At a  FIG. 9. Test for convergence of the Lyapunov exponents with the simulation domain size.The vertical axis shows the Lyapunov exponents for three simulations with domain sizes of 16, 32, and 64 cm in the ẑ direction and 1 cm in the x and ŷ directions.The ẑ direction has 16 cells per centimeter, while the x and ŷ directions each have one cell per centimeter.Each cell contains 24, 088 particles.The Lyapunov exponents remain relatively stable as the simulation size increases.

Domain size
To investigate the impact of domain size on the Lyapunov exponent, we conducted three simulations with dimensions of 16, 32, and 64 cm in the ẑ direction and 1 cm in the x and ŷ directions.There are 16 cells per centimeter in the ẑ direction and one in the x and ŷ directions.In each cell, 24, 088 particles are emitted in an approximately isotropic distribution from the center.A perturbation with | ⃗ δ t0 |/|⃗ r t0 | ∼ 10 −10 is applied at t = 2.65 ns long after the saturation of flavor instability.We allowed the perturbation to evolve until the simulation reached 20 ns or until | ⃗ δ t |/|⃗ r t | ∼ 10 −6 .
Figure 9 shows the dependence of Lyapunov exponents on the size of the simulation domain.The Lyapunov exponents do not exhibit significant variation and remain approximately constant as the simulation domain increases.We consider convergence achieved for a simu-lation domain of 1 × 1 × 64 cm.

Cell size
To investigate the impact of the cells size on the Lyapunov exponents, we perform four simulations with 128, 256, 512, and 1024 cells in the ẑ direction and one in the x and ŷ directions.The domain size is 1 × 1 × 64 cm, and there are 24, 088 particles per cell.We Figure 10 shows that the Lyapunov exponents are slightly affected by the cell size.We consider convergence achieved for 1, 024 cells within a domain size of 1 × 1 × 64 cm.

FIG. 2 .
FIG.2.Density matrix averaged over the spatial domain of the NSM snapshot (upper panel) and Fiducial (lower panel) simulations.Individual components are visualized in distinct colors.In the off-diagonal components, flavor oscillation modes that start the simulation with small amplitudes grow, creating a combination of modes that stop the exponential trend due to unstable modes.The diagonal components reach a final state near complete flavor mixing (indicated by the black dotted line) and slightly fluctuate around it.The domain-averaged antineutrino density matrix follows the same trend.
FIG.3.Time evolution of small perturbations | ⃗ δt| in the NSM (blue) and Fiducial (orange) simulations.In the left panel, the perturbation is applied at t0 = 0 ns (before the linear growth phase).Between 0.10 and 0.25 ns, the perturbation growth is exponentially driven mainly by the fast flavor instability.In the right panel, to avoid the presence of unstable fast flavor modes, the perturbation is applied in the decoherence phase at t0 = 2.65 ns.The magnitude of the perturbation grows, following an approximate exponential trend.This demonstrates that paths of similar flavor vectors in the state space diverge exponentially, illustrating the chaotic evolution of neutrino flavor.As the simulation concludes, the magnitude of the perturbations approaches a constant value larger than the magnitude of the flavor vector |⃗ rt| (dotted lines).

FIG. 4 .
FIG.4.The upper panel shows the time evolution of perturbations periodically normalized to keep the perturbation small.The blue curve corresponds to the NSM simulation, while the orange curve represents the Fiducial simulation.The lower panel shows the Lyapunov exponent calculated over each interval between normalizations.The Lyapunov exponent of a perturbation in the state space varies based on the direction of the perturbation.The presence of exponential divergence (λ > 0) and convergence (λ < 0 e.g.below the blue and orange stars) is an expected result of the Liouville theorem of conservation of the state space volume of canonical variables.

FIG. 6 .
FIG.6.The upper panel shows the flavor-traced difference between the neutrino density matrices (Equation5) of single computational particles extracted from the baseline and postsaturation perturbed simulations.The lower panel shows the same quantity, but using the domain-averaged density matrix ⟨ρ⟩.The blue and orange curves denote the NSM and Fiducial simulations, respectively.The black dotted lines show the maximum possible value.Both panels exhibit a comparable exponential trend, akin to the perturbations observed in the flavor vectors (see the right panel of Figure3), but late-time uncertainties in the domain-averaged quantities are much smaller than in individual particles.
), we employed the least squares linear fit method for y = mx + b using the numpy.polyfit()function.Here, y = ln | ⃗ δ t |, m = λ, x = t and b = ln | ⃗ δ t0 |.The numpy.polyfit() function returns values of m and b that minimize the sum of the squares of the residuals.

FIG. 7 .FIG. 8 .
FIG. 7. Test for convergence of the Lyapunov exponent with the number of particles in the EMU simulations.The vertical axis shows the Lyapunov exponent for five Fiducial simulations.The simulation domain is 1 × 1 × 64 cm divided into 1024 cells.Convergence is attained for a number of directions greater than 24, 088.

FIG. 10 .
FIG.10.Test for convergence of the Lyapunov exponents with the cell size.The vertical axis shows the Lyapunov exponents for four simulations with 128, 256, 512, and 1024 cells in ẑ and one in x and ŷ directions.The domain size is 1 × 1 × 64 cm, and there are 24, 088 particles in each cell.The number of cells in the simulation has a slight effect on the Lyapunov exponents.In this work, we consider convergence achieved for 1, 024 cells in a domain size of 1 × 1 × 64 cm.