Equilibration of quantum many-body fast neutrino flavor oscillations

Neutrino gases are expected to form in high density astrophysical environments, and accurately modeling their flavor evolution is critical to understanding such environments. In this work we study a simplified model of such a dense neutrino gas in the regime for which neutrino-neutrino coherent forward scattering is the dominant mechanism contributing to the flavor evolution. We show evidence that the generic potential induced by this effect is non-integrable and that the statistics of its energy level spaces are in good agreement with the Wigner surmise. We also find that individual neutrinos rapidly entangle with all of the others present which results in an equilibration of the flavor content of individual neutrinos. We show that the average neutrino flavor content can be predicted utilizing a thermodynamic partition function. A random phase approximation to the evolution gives a simple picture of this equilibration. In the case of neutrinos and antineutrinos, processes like $\nu_e {\bar{\nu}}_e \leftrightarrows \nu_\mu {\bar{\nu}_\mu} $ yield a rapid equilibrium satisfying $n( \nu_e) n({\bar \nu}_e) = n( \nu_\mu) n({\bar \nu}_\mu) = n( \nu_\tau) n({\bar \nu}_\tau)$ in addition to the standard lepton number conservation in regimes where off-diagonal vacuum oscillations are small compared to $\nu-\nu$ interactions.


I. INTRODUCTION
In hot and dense astrophysical environments such as core collapse supernovae and binary neutron star mergers, neutrinos are emitted in such large number fluxes that the average flavor content is important to the dynamic and chemical evolution.Through weak interactions with local nucleons electron flavor neutrinos and antineutrinos can alter the local proton-to-neutron ratio in these environments thereby affecting, for example, r-process nucleosynthesis [1][2][3][4][5][6][7].It is therefore crucial to understand the flavor content of the neutrinos if one wishes to perform detailed studies of the evolution of these systems [8][9][10][11][12][13][14][15].
In regions of high density, neutrinos can undergo coherent forward scattering which is sensitive to the quantum mechanical flavor state of the neutrino.When scattering on charged leptons, a flavor dependent relative phase can develop between components of the flavor state.The neutrinos can also exchange flavor content with other neutrinos through neutral current coherent forward scattering.When the number density of neutrinos is sufficiently large, this flavor exchange effect is expected to dominate the flavor evolution of the neutrino gas, and novel coherent effects such as flavor spectrum splits and swaps may occur [2,[8][9][10][11][12][16][17][18].
In this work, we will study the flavor evolution of such a neutrino gas in the limit that the potential generated by coherent ν − ν forward scattering is the dominant effect; this is a regime known in the literature as the "fast" flavor oscillation limit [14,[19][20][21][22][23][24][25][26][27][28][29][30][31][32][33].In this regime, the redistribution of flavor content among the neutrinos is expected to occur on length scales of approximately a meter.We will study this regime in a quantum manybody formalism with a parametrization of the coherent forward scattering potential which does not impose any simplifying symmetries on the momenta of the neutrinos.
In a core-collapse supernovae the neutrino density is governed by n ν = L/( Ē 4πR 2 ).Taking a total luminosity L = 10 53 erg/sec, R = 50 km, and an average energy Ē = 10 MeV gives a total density n ν ≈ 6.6 × 10 −7 fm −3 .For a degenerate relativistic Fermi gas, the average energy is given by Ē = 3E F /4 where E F is the Fermi energy.For a two species Fermi gas (e and µ flavor neutrinos) the number density is given by n = E 3 F /(3π 2 ).Equating the average neutrino energy to the Fermi energy, we find that the density of a degenerate, zero temperature twocomponent Fermi gas is approximately 15 times greater than the above estimated neutrino density.For 6 neutrino species (3 neutrino and 3 antineutrino species) with roughly equal fractions, the density is approximately 45 times lower than a degenerate gas of the same average energy.As we estimate the effects of degeneracy to be minimal, we use Boltzmann statistics for both the neutrinos and antineutrinos.We will also work in the twoflavor approximation such that a single neutrino's flavor quantum state can be written as an SU(2) spinor.
The ν − ν coherent forward scattering potential takes the form of an all-to-all Heisenberg-like interaction, and such Hamiltonians are generically expected to be nonintegrable except in special cases.We will provide evidence that this interaction is non-integrable in the absence of simplifying symmetries, and that we observe a characteristic signature of (non)integrability when simplifying symmetries are imposed (relaxed).Furthermore, non-integrable Hamiltonians are hypothesized to "thermalize" in the sense that expectation values of few-body operators are expected to equilibrate to a value which can be predicted from an appropriate thermodynamic partition function.For the generically parameterized potential we observe strong agreement between the one-body flavor expectation values obtained in the exact many-body evolution of the system and those predicted from a grand-canonical Boltzmann distribution.

II. ENERGY SPECTRUM
The general Hamiltonian governing the coherent evolution of the neutrino flavor content is composed of three parts (see e.g.[34]).The first is the vacuum potential which stems from the fact that the neutrino mass states are linear combinations of flavor states.The second is the matter potential [35,36] and it is generated by coherent forward scattering through charged current interactions with local charged leptons in the environment.Finally is the ν − ν potential generated through neutral current coherent forward scattering between neutrinos [8,15].For the rest of this work, we will assume that the vacuum oscillation potential is negligibly small relative to the other two potentials, and that its primary effect is to provide perturbation to an otherwise pure flavor product state initial condition.We also assume that the matter profile is uniform in the regions of the environment under consideration, and as such we may consider a co-rotating frame to remove its overall effect on each neutrino.
After these modifications, only the the ν − ν coherent forward scattering Hamiltonian remains.For N neutrinos it has the form [15] Here ⃗ σ is the usual vector of Pauli operators, and µ = √ 2G F n ν is the scale of the ν − ν interaction.Given a typical core-collapse supernovae flux at a radius of 50 km a simple estimate gives µ ≈ 1 cm −1 .As it is the only dimensionful parameter in the problem, we hereafter measure all distances and times in units of µ and set µ = 1.
We perform simple numerical simulations for small N in periodic boundary conditions.This is a toy model because an estimate using a typical flux for 20 neutrinos at this density gives a box size of order 300 fm, much too small to be realistic.This is a standard choice in the literature, [13,14,31,[37][38][39][40][41][42][43][44][45][46][47][48][49][50][51][52], but imposes a great deal more symmetry than a more realistic case.Relaxing this condition would only induce more decoherence in the simulations, a critical ingredient as we discuss below.(For more discussion on the role of decoherence see e.g.[37,38,46,53].) We will address two commonly employed simplifying symmetries in the literature.The first is the uniform coupling symmetry, under which it is assumed that the v i • v j term averages to zero under the action of the allpairs summation in Eq. (1).For this case, the total H νν is proportional to J 2 and is straightforwardly diagonalized in the absence of the other two relevant potentials.For some choices of one-body potentials, the uniform couplings Hamiltonian may be diagonalized with an alge-braic Bethe ansatz [41,51].This simplification has been well studied so we will not consider it further.
Another simplifying assumption is to impose an axial symmetry to the velocity coupling term such that the velocity components orthogonal to the momentum symmetry axis (denoted z) average to zero.This is equivalent to making the substitution It was recently shown that this symmetry imposition leads to an integrable Hamiltonian which can also be diagonalized with an algebraic Bethe ansatz [52].
There is no known method for systematically diagonalizing the generic Hamiltonian of Eq. ( 1) through the use of an extensive set of nontrivial conserved charges.A "nontrivial" operator is an operator which is not simply proportional to a projection operator into an eigenstate of the Hamiltonian.The only obvious operators which commute with Ĥνν are the projections of the total angular momentum (where a ∈ [1,3] is the vector index of the Pauli matrices) and the square of the total angular momentum ( Ĵ2 = 3 a=1 Ĵa Ĵa ).Because Ĥνν commutes with Ĵ2 , Ĵ3 , and the ladder operators Ĵ± , the energy eigenstates can only be linear combinations of |j, m⟩ in these degenerate subspaces.In order to obtain a the full energy spectrum, one can "brute force" diagonalize Ĥνν in the lowest |m| subspace for each j.Because Ĥνν commutes with the angular momentum ladder operators, the energy spectrum is 2j + 1-fold degenerate for each j, and eigenstates with higher values of |m| can be systematically constructed from the minimal |m| states with suitable applications of the ladder operators.
In the following, we construct the Hamiltonian by first drawing the z components of the velocities on the interval v z ∈ [−1, 1] randomly from a Gaussian-like distribution centered on v z = 1 with a standard deviation of 1/2 such that P (v z ) ∝ e −2(vz−1) 2 .We choose this distribution to approximate the forward-peaked-ness of the neutrino momenta emitted from the proto-neutron star surface.For the axially symmetric case, the x and y components are chosen to be identically zero.For Ĥνν , which retains all three components of the unit velocities, we choose where ϕ is a chosen to be a random number uniformly drawn in the interval [0, 2π].
Finally, we will solely treat the evolution of initial product states.This necessarily assumes that the neutrinos entering a given region of space have developed a negligible correlation in their flavor evolution prior to the interaction we herein consider.We make this approximation only in the interest of analyzing the salient dynamics of the evolution under Ĥνν .More careful treatment of the initial state will require a detailed study of the neutrino-neutrino self interaction in the presence of non-negligible momentum changing scattering processes in the dense matter of the proto-neutron star core.

A. Level spacings
We next consider the structure of the spectrum of Ĥνν in its m and j symmetry sectors.In each sector we sequentially order the energy eigenvalues and investigate the differences between the sequential energies which we denote s α ≡ (E α+1 − E α )/s where s is chosen such that the average over the energy differences is unity.
Integrable systems are characterized by an extensive set of (nontrivial) operators which commute with the Hamiltonian, and these conserved quantities permit substantial degeneracy in the energy spectrum.When considering the sequential energy differences for an integrable system, the probability of two sequential energy eigenvalues having spacing s is a Poisson distribution, P (s) = e −s .In contrast, non-integrable Hamiltonians demonstrate repulsion between energy levels, and the probability distribution for spacing s is, in accordance with the Wigner surmise, of the form P (s) = π 2 se −πs 2 /4 for a Gaussian orthogonal ensemble (GOE) [54].
When extracting the distribution of energy differences from a particular Hamiltonian, an unfolding procedure must be performed in order to approximately remove the effect arising from the fact that the Hamiltonian has a finite density of states which may be sparser in some energy regions than others [55].Instead of considering directly the distribution of level differences s i , one can consider the probability distribution for the ratio of sequential level differences, defined as r α ≡ sα+1 sα which is insensitive to the local density of states.The new probability distribution, P (r), for the integrable (Poisson) case is While for non-integrable (GOE) cases, P (r) is . ( For the GOE case, the expression in Eq. ( 5) was derived as an exact expression for 3 × 3 random matrices drawn from a Gaussian Orthogonal Ensemble.It is in good agreement with numerical studies of larger random matrices, and the correction derived from fits to numerical distributions extracted from larger random matrices can be found in [55].The approximate correction for the large matrix case is O(10 −2 ) and is not visible on the scale plotted in Fig. 1.
In Fig. 1 we compare the distribution of r i for the two Hamiltonians we have discussed.The Hamiltonian, Ĥνν , implements generic unit magnitude three velocities, while the axially symmetric (abbreviated "ax") Hamiltonian, Ĥax νν , implements the identical v z,i values as Ĥνν , but takes the x and y components of the velocities equal to zero by construction.As previously mentioned, Ĥax νν is known to be integrable, and therefore the level-spacings in its spectrum are expected to obey Poisson-like statistics.For both cases, we use 16 neutrinos in the construction of the Hilbert space.
We diagonalize Ĥνν and Ĥax νν in the j = 2 subspace and plot the normalized histogram of ratios of level spacings (r α ).We also compare directly with the expected universal probability distributions for the integrable and nonintegrable cases.This figure provides evidence that the generically parameterized Hamiltonian does not have an extensive set of "hidden" conserved charges corresponding to some hitherto unrecognized symmetry.
To validate that the level spacing statistics in the other j subspaces behave similarly, we compute the average rα with respect to the extracted probability distribution and show the results for j ≤ 4 in table I.The distribution of ratios of level spacings should not depend on the subspace, so this average value should take a universal value.For the GOE case, it can be shown that the value is rα ≈ 1.7789 (see table 1 of [55]).However, for the Poisson case the average diverges as r α → ∞, so for any finite distribution of level spacings the average can take arbitrary and unbounded values.

III. ENERGY DEPHASING
The distribution of the ratio of level spacings is in good agreement with the probability distribution one would expect if the Hamiltonian under consideration was a ran- dom matrix drawn from a GOE.The level repulsion resulting from a lack of an extensive set of symmetries implies that energy differences should not be expected to vanish for the bulk of states in the spectrum.Furthermore the effective random nature of the Hamiltonian implies that differences in energies should not generically be simple rational numbers.As such, the time-dependent phases of off-diagonal matrix elements of operators in the energy basis should be expected to individually average to zero over arbitrary time windows after the state has evolved for a sufficiently long time.Thus, under a time average over an arbitrary interval at late time, the expectation value of σ3,i should be in good agreement with the value predicted by computing the average with respect to the incoherent energy diagonal distribution which preserves the energy state probabilities of the initial condition.Thus we will compute and compare the approximation where |ψ 0 ⟩ represents the initial quantum state and t eq > 1/µ is the (as yet unknown) time it takes for the system to reach its equilibrium value.The time window [t eq , t eq + ∆] must also be selected such that the oscillations due to the coherent evolution of the quantum system about the average expectation value can be fully captured.We will refer the value t eq as the equilibration timescale, and we will discuss our observations of it in the next section.We refer to the incoherent probability weighted sum of operator expectation values in the final line of Eq. ( 6) as the energy mixed state distribution (EMSD).

A. Thermalization
Non-integrable many-body Hamiltonians are widely expected to obey the so-called "Eigenstate Thermalization Hypothesis" (ETH) [56][57][58], see [59] for a review.Many numerical studies of ETH have been conducted, we refer to Refs.[60][61][62] as a very incomplete list of examples for when the system has an underlying lattice structure.For systems with all-to-all interactions or no discernible lattice structure, there is a less extensive literature, and we refer to Refs.[63][64][65] as examples.Besides integrable systems, for a generic Hamiltonian, ETH may only hold for the majority of the eigenstates.Certain eigenstates (quantum scars) embedded in the spectrum may have additional symmetries beyond that of the Hamiltonian, invalidating ETH for that state, see for instance [66,67].However, these states are expected to have measure zero in the density of states in the thermodynamic limit, unless there is special symmetry protection.Loading the Hamiltonian with additional symmetry, while not inducing integrability, results in the fracturing of the Hilbert space into many block-diagonal subsectors, with each subspace being simply too small for ETH [68].
The ETH is fundamentally a statement about the structure of "few-body" operators that act on a subextensive number of the local Hilbert spaces that tensor together to form the full many-body Hilbert space.For a "few-body" operator Ô, the matrix element between two eigenstates |E α ⟩, |E β ⟩ of the Hamiltonian is hypothesized to be given by the ansatz: where , and ω αβ = E α − E β .Ô(E ave ) is the microcanonical average of the operator Ô over states near E ave , and S is the suitable microcanonical entropy, counting the number of eigenstates near E ave within a small window given by ϵ.R αβ is an independent random number for each α, β with average zero, but variance 1.When the dimension of the total Hilbert space is large, we can expect S to be large even for very narrow windows, so that in the thermodynamic limit, we can take ϵ → 0 and have the matrix element of the operator in eigenstates be dominated by the microcanonical average.Finally, the function f Ô is related to the linear response correlation function of the thermodynamic system when perturbed by the operator Ô [69] and, for fixed E ave , decreases exponentially in the energy difference ω αβ for such chaotic systems (see e.g.[70,71]).
In general we expect the ETH ansatz (Eq.( 7)) for the matrix elements to hold for a chaotic system but not for an integrable one.As shown above, the neutrino Hamiltonian Ĥνν in Eq. (1) we use can reproduce both behaviors depending on the choice of the velocities.To test the validity of ETH for our problem we consider matrix elements of σ3,i on the first neutrino in both the integrable and non-integrable regime.We first show the diagonal matrix elements ⟨E α | Ô|E α ⟩ in Fig. 2 for both the integrable (red data on right panel) and non integrable case (blue data on left panel).One can see that in the former (integrable) case the diagonal matrix element has large fluctations over nearby frequencies throughout the whole spectrum while for the non-integrable case the expectation value in the bulk of the spectrum is approximately a function of energy alone as expected from the ansatz in Eq. (7).Next, in Fig. 3 we show the magnitude of the off-diagonal matrix elements of σ3,i on the first neutrino as a function of the energy difference |ω αβ | for all states in a narrow energy window of size ϵ = 0.005[µ] around an average energy E ave in the middle of the spectrum: for the results shown here we took E ave = 2.4[µ].One can clearly see an exponential decay for large energy difference in the non-integrable case (left panel) while for the integrable model (right panel) the size of the off-diagonal matrix elements fluctuates by more than six orders of magnitude for every value of |ω αβ |.To better visualize the importance of a large number of outliers in the integrable case, which are instead not present in the nonintegrable system on the left, we also report in Fig. 3 the running average of the matrix element size obtained by taking the median over a window of 50 matrix elements.This is a direct estimate of the absolute value of the function f Ô (E, ω) in Eq. ( 7) for E in the middle of the spectrum (cf.[59,72]).
All results shown in Fig. 2 and Fig. 3 were obtained for a system with N = 16 neutrinos and for states restricted to j = 2 and m = 2 subspace.
The basic intuition behind eigenstate thermalization is that few-body operators are unable to distinguish which eigenstate the system is in, with nearby (in energy) eigenstates displaying similar few-body behavior.When the dimension of the Hilbert space becomes large, the eigenstate thermalization is often described using canonical or grand-canonical ensembles, since the expectation values are expected to be effectively equal in either ensemble.When using the (grand)-canonical ensemble, the temperature (or other chemical potentials) are tuned to match the energy and other quantum numbers of the desired eigenstate of the system.After this, the matrix elements of generic few-body operators can be computed in the appropriate thermal ensembles.When ETH holds, we can effectively claim that up to small corrections, the time-averaged spectrum of few-body operators "thermalize" even in a single eigenstate, and are thus computable from thermal ensembles.One can generalize the eigenstate thermalization to any state of the system, provided expectation value of the variance of the energy for that state is appropriately small.
If our system "thermalizes" in this sense, then we can predict the (time-averaged) expectation value of the individual neutrinos at large times utilizing a grandcanonical statistical distribution.As the system is time independent, and Ĵ3 and Ĵ2 commute with the Hamiltonian, we construct a partition function using one temperature parameter (β) and two chemical potentials (µ 3 and µ 2 ).For a given initial condition, we find these parameters by fitting the expectation values of the relevant operators calculated using the partition function to the invariant expectation values calculated with respect to the initial condition.This amounts to finding β, µ 3 , and µ 2 such that for the conserved quantities Ô = Ĥνν , Ĵ3 , Ĵ2 and Z = Tr e −β Ĥνν +µ3 Ĵ3+µ2 Ĵ2 is the partition function.Once the temperature and chemical potentials have been determined, the (time-averaged) expectation value of the flavor content for a given neutrino (i) can be predicted by computing ⟨σ 3,i ⟩ ≈ Tr σ3,i e −β Ĥνν +µ3 Ĵ3+µ2 Ĵ2 .We next investigate the late time behavior of an example product state initial condition for 16 neutrinos.We employ the same couplings in Ĥνν as in the previous sections, and we order them from lowest to highest value of v z,i in increasing values of the particle index i ∈ [1,16].We next choose the first i ∈ [1, ⌊N/2⌋ + 1] neutrinos to be electron flavor (corresponding to the lowest 9 v z values for 16 neutrinos), and the remaining (highest v z ) neutrinos to be τ flavor.Finally, to mimic the effects of the neglected one body contributions to the Hamiltonian, we perturb the initial flavor configuration by a small random polar rotation away from initial pure flavor states, as well as a random rotation in the x−y plane of the flavor Bloch sphere.The first rotation mimics the effect of the small effective mixing angle induced by noncommutation of the dense matter and vacuum potentials, while the second rotation mimics the phase accumulation from the rapid rotations about the flavor axis induced by the matter potential.
Once the couplings and initial conditions are specified, we then evolve the 16 neutrino quantum state by numerically solving the many-body Schrödinger equation in the interval t = 0, 10 3 µ −1 .Using the time evolved quantum state, we compute the average of ⟨σ 3,i ⟩ over the interval t = [10 1 , 10 3 ]µ −1 .We show this time averaged quantity in the top panel of Fig. 4 and we compare time averaged numerical predictions for ⟨σ 3,i ⟩ with the mixed state predictions for each spin.The red squares represent the time-averaged values of the evolved quantum state, and the green circles (blue diamonds) show the expectation value estimated using the EMSD (grandcanonical) approximation of Eq. ( 6) (Eq.( 9)).The bottom panel shows the difference between the late time averaged many-body solution and the two energy diagonal predictions, while the filled red region indicates the rootmean-square (RMS) deviation of the time oscillations of the expectation value of the many-body solution about its mean (as seen in Fig. 5).We observe that both of the energy diagonal approximations for the estimation of ⟨σ 3,i ⟩ are fully within the RMS oscillations of the the many-body solution about its mean, excepting one outlying point estimated using the grand-canonical ensemble.
Its important to note that while the interaction is SU (2) invariant, we choose a particular spin direction along which to quantize (which is informed by the initial polarization) and denote it as ⃗ e 3 .Because Ĥνν commutes with Ĵ+ and Ĵ− their expectation values are also conserved, but these operators cannot be simultaneously diagonalized with Ĵ3 and Ĵ2 .Therefore no mixed state which is energy diagonal will be able to accurately predict their expectation values.For the initial conditions we consider here, ⟨ Ĵ+ ⟩ ≈ ⟨ Ĵ− ⟩ ≈ 0. As such, the energydiagonal mixed states we consider here are adequate for describing the polarization configurations.However in the presence of substantial orthogonal polarization (i.e.substantially nonzero values of ⟨ Ĵ± ⟩), a more careful treatment of the non-abelian conserved charges would be necessary [73].

B. Approach to equilibrium
In the previous sections we have argued that in the absence of simplifying symmetries the ν − ν interaction Hamiltonian in its various symmetry sectors has a level spacing distribution which is in good agreement with that of a random Hamiltonian drawn from the GOE.We have also argued that this implies generic product state initial conditions should display substantial dephasing in energy such that one-body expectation values can be computed from an energy-diagonal distribution.Once dephased, the one-body expectation values are expected to achieve an equilibrium value from which they stray only transiently and with small amplitude.
To investigate the approach to this equilibrium, we evolved seven different systems each with N spins from N = 10 to N = 16 in integer steps.For each system size we independently selected velocities and initial conditions in the same manner as described for the N = 16 case shown in detail above.We then considered three simple measures of the speed of the evolution.First, for each system size (N ), we consider the one-body von Neumann entropy (S i = −Tr(ρ i log 2 (ρ i ))) of each neutrino, and we find the earliest time for which S i > 0.95.We then take the average of these times over all the spins for a given system size, we denote this average large entropy time as T S , and plot T S as black diamonds in Fig. 6.We also computed the standard deviation of T S over the spins for each system size, however in each case the error bar would be too small to resolve on the plotted scale of Fig. 6.Next, we note the time dependent behavior Fig. 5 of P 3 ≡ ⟨σ 3 ⟩.After evolving away from the initial value, the expectation value crosses the thermal prediction before turning and approaching it again.We observe this behavior across all of the chosen system sizes and for each spin.We therefore consider, for each system size, the time at which the value ⟨σ 3 ⟩ for each spin reaches a turning point at which its first time derivative vanishes.We then average these turning-point times over all the spins, denote the average time as T P3 , and plot it as black circles in Fig. 6, and we indicate the standard deviation from the average value over the spins in a given system size with error bars.
Finally, we consider the Loschmidt echo L(t) = |⟨ψ 0 |e −i Ĥνν t |ψ 0 ⟩| 2 which quantifies the probability of measuring the time evolved state in the initial configuration.The t = 0 curvature of the Loschmidt echo is given by the variance of the Hamiltonian, and the earliest time at which it (the Loschmidt echo) can vanish is bounded by the quantum speed limit [74].For chaotic systems, the Loschmidt echo is expected to saturate to an equilibrium value which scales inversely proportionally to the size of the system's Hilbert space [75], a scaling behavior that we have verified for Ĥνν .The dynamics of the echo at intermediate times is an active area of study, but we investigate the time of the first minima in the echo as a proxy to the equilibration timescale of the system.For all of the calculations we have performed, we have observed an approach to an equilibrium value in the onebody flavor content with a timescale which appears insensitive to the total number of spins.The equilibration we observe in flavor content is subsequent to a development of one-body entanglement on a timescale which is similarly insensitive to N .Our numerical observations are consistent with a time-to-equilibrium which scales simply as O(µ −1 ).In Fig. 6 we show the timescales we extracted for the equilibration of these quantities for a range of system sizes.Because our computational method retains all amplitudes, we are computationally limited in the system sizes N we can investigate.
For states that are close to polarized product states in the ⃗ e 3 direction, we can argue that the time-scale to equilibrium cannot scale to zero as N ≫ 1.In the Appendix, we calculate to quadratic order in the Taylor series expansion of the time evolution of ⟨σ 3,i (t)⟩ about t = 0. We find that, for the typical states we consider here, the linear term in the expansion vanishes, and the quadratic term scales as µ 2 /N .We expect such a Taylor series to have a radius of convergence scaling as µ −1 .Noting ⟨σ 3,i (t = 0)⟩ ≈ ±1, then truncating the series at quadratic order and estimating when ⟨σ 3,i (t)⟩ ≈ 2⟨ Ĵ3 ⟩/N ≪ 1, we would conclude t ∼ √ N /µ.
However, such a time-scale is outside the expected region of convergence for our series, making such a conclusion not self-consistent.We can, though, conclude that t → 0 as N ≫ 1 is impossible: the smallness of the quadratic term, and the suppression of higher orders as t → 0, would not allow a self-consistent solution to the equilibrium condition in arbitrarily small neighborhoods of the origin.We intend to follow up these observations with investigations of these timescales both analytically and with computational methods which may allow substantially larger system sizes, such as tensor network methods.

C. Path integral description of time evolution and equilibrium distributions
In the results presented in the previous section, it is striking that at late times each individual neutrino flavor expectation value is approximately given by ⟨σ 3,i ⟩ ≈ 2⟨ Ĵ3 ⟩/N .We observe this approximate flavor isotropy even when there is substantial correlation between the initial flavor content and momenta of the neutrinos (as in the initial split configuration of Fig. 4).This is behavior we have observed across a variety of system sizes, and initial correlations between flavor and momenta.The total ⟨ Ĵ3 ⟩ is of course conserved by the Hamiltonian, but there is no a priori reason to expect all the spins to reach the same equilibrium value.To clarify this point and to define the equilibrium values for general systems of neutrinos (and both neutrinos and anti-neutrinos), we consider the time evolution in a path integral approach.
The time-evolved expectation value of an operator expressed as an expansion of overlaps with the initial state for a system described by the GOE produces random phases in the off-diagonal elements that would normally be present in Eq. (6).Rewriting the expectation value as a sum over product states of ν e and ν τ spins |n⟩ the amplitude of the state in state |n⟩ as a function of time is: where we have taken the initial state |ψ 0 ⟩ as an arbitrary state in the |n⟩ basis.Equivalently, writing the state |ψ(t)⟩ as a sum over the states |n⟩ where the A n (t) are the magnitudes of the overlaps (real and positive) and the phases are ϕ n (t) .The magnitudes and phases are given by: The expectation value of σ3,i averaged over a time window of size ∆ is: and we have used the facts that σ3,i is diagonal in the |n⟩ basis and that the average over non-diagonal energy eigenstates over time goes to zero.This is a sum over diagonal matrix elements in the |n⟩ basis, each with a positive coefficient, suggesting the time-evolved state can be written as an incoherent sum of the states |n⟩.
For an initial state with overlaps with many eigenstates, the random phases (e −iEαt ) for a GOE would translate to random phases in the ν e − ν τ product state basis ϕ n (t).This would not be true near the ground state where the phases cannot be random, but initial product states we wish to describe are not near the ground state of this Hamiltonian.Random phases also keep ⟨σ ±,i ⟩ = ⟨σ ∓,i ⟩ = 0 for all times for an initial product state in |n⟩.
The Hamiltonian can always be divided into diagonal ( Ĥd νν ) and off-diagonal ( Ĥod νν ) pieces, where, in the |n⟩ basis, the off-diagonal pieces are proportional to (1 which is simply a weighted permutation operator exchanging anti-parallel flavor spins i and j.The full path integral describing the propagation with exp[−i Ĥνν t] is then a sum over all paths with all-to-all two-body spin exchanges while the diagonal piece just reduces to pure phases, as do two-body exchanges between parallel spins. We can write the path integral as a sum over all powers of Ĥod νν operators.Summing over the number of nondiagonal operators at random times, with diagonal evolution (pure phases) between the non-diagonal terms, we can rewrite the time evolution with a path integral as Here, we can see that Ĥd νν introduces pure phases into each term of the path integral, while the off-diagonal terms Ĥod νν induce transitions from one product state to another.In this equation, the sum of all times t i has a resultant of t, and we have separated the path integral into terms with a specific number of insertions (indexed by m) of the off-diagonal operators.
In principle one could sample these paths by their absolute magnitudes as is done in many quantum Monte Carlo (QMC) approaches, in particular diagrammatic Monte Carlo [76].If the initial product state connects to many basis states and the Ĥd νν introduces random phases, the final state is an incoherent sum of product states.The off-diagonal terms Ĥod νν induce transitions as a series of spin swaps between one product state and another.Such a sampling of swaps conserves the expectation value of Ĵ± , Ĵ3 and Ĵ2 of an initial product state since exchanges of the original spins have the same expectation values of total and projected spin as the original state.
All permutations of the original spins can be reached by a series of two-body swaps.For an initial product state of many neutrinos, we can easily calculate the expectation value of ⟨ Ĥk νν ⟩ for small k.These quantities must be conserved by the time evolution.If the angular distributions are similar, as they are deep inside a proto-neutron star, all product states of permutations of the initial spins should, when averaged over time, have equal magnitudes of overlaps A. Each permutation of the initial product state will have approximately the same ⟨ Ĥd νν ⟩, with a variance inversely proportional to the number of neutrinos since there are N 2 terms in the Hamiltonian.Unitarity requires that the average of the squared absolute magnitudes is inversely proportional to the number of spins.
More generally the time-averaged magnitudes of the amplitudes are expected to depend on the ⟨ Ĥd νν ⟩ of the individual permutations.These can arise, for example, from interference between different insertions of off diagonal operators leading to the same final state.If we assume time-averaged overlaps vary smoothly with ⟨ Ĥd νν ⟩, we can calculate this dependence by requiring that the lowest N moments of the Hamiltonian are conserved.The linear dependence could be parametrized as a temperature as in ETH, because we also know that the absolute magnitudes of the time-averaged amplitudes are real and positive.Higher order moments put further constraints on the evolution.Example calculations for small numbers of neutrinos are discussed below.Crucially the low order moments of Ĥνν can be calculated exactly from the initial product state, or from an incoherent sum of physically reasonable product states.
The path integral of Eq. ( 14) can be described as a random walk in the basis states.For large times we expect the phases of individual ν e − ν τ product states to become random if the initial product state is in the middle of the spectrum.The short time evolution of an amplitude of a particular state (a n (t) ≡ ⟨n|ψ(t)⟩) up to order δt is governed by There are O(N 2 ) terms in the sum over states m ∈ P (n) where P (n) is the set all of states which can be produced by one pair permutation in Ĥod νν , each matrix element with a substantial random component.
Assuming phases of individual ν e − ν τ product states at large enough time separations are random, the expression for the path integral can be described as a random walk in the basis states.The complex amplitude for a specific state at a given large time will be distributed as a two-dimensional Gaussian centered at zero describing its real and imaginary parts.The average magnitude (squared) is governed by unitarity.Integrating over times after equilibration should produce a constant absolute magnitude of the overlap, with a variance decreasing approximately inversely proportional to the square root of the time integrated over.Expectation values are obtained as the incoherent sum of expectation values in individual purely ν e -ν τ product states.
We investigated the behavior of the time average of |A n (t)| at late times exactly for small N .The above discussion suggests that we should expect where the t subscript indicates a time average, and N m is the total number of states in the Hilbert space with quantum number m = ⟨n| Ĵ3 |n⟩.In Fig. 7  √ N m .Both of the considered cases show some structure in the magnitudes versus ⟨n| Ĥνν |n⟩, as they must in order to conserve all of the moments of the Hamiltonian.For a more complex Hamiltonian with terms violating the conservation of J 2 or components of J or for time-varying Hamiltonians the symmetries are further reduced [77], making the approximations considered here even more accurate.Even for the case of a specific product initial state and a static Hamiltonian with these symmetries the behavior of the absolute magnitudes A is smooth in energy.The equilibrium distribution in energy and angle reached by a pure neutrino system is simply the weighted average of the initial distributions in energy and angle.For the simple two-body interaction discussed here, the total number of neutrinos of each flavor is conserved.This equilibration would produce the horizontal dashed lines in Fig. 4 and is reached quite quickly, as we discussed previously.
For an incoherent sum over orthogonal states, the path integral can be modeled as a classical process including exchanges of pairs of spins.A classical swap network can be implemented for a large number of neutrinos.The potential dependence of the amplitudes on ⟨ Ĥνν ⟩ could be implemented through a Metropolis Monte Carlo with the parametrization of the energy dependence determining the accept/reject probability.For cases with minimal dependence, such as equal initial angular distributions, the accept/reject probability could be fixed to 1.
Because of the lack of coherence, it could be that nonforward scattering would also be important.The magnitude of these amplitudes would be similar to the forward scattering.However the phases from different magnitudes of the neutrinos' momenta would oscillate rapidly on the MeV −1 scale, rendering their contributions very small.In principle all amplitudes should be consistently included and summed but these highly oscillatory pieces are unlikely to significantly interfere with the forward scattering amplitudes.The same picture of random diagonal phases and rapid flavor exchanges would occur and not change this picture, particularly since the angular dependence is not measurable.At larger distances where the neutrino flux is reduced, off-diagonal vacuum oscillation terms will be important, leading to further flavor evolution.Potentially this could be modeled as a system of single neutrino spins evolving according to their onebody Hamiltonian with random swaps imposed.However the equilibrium we discuss should dominate at smaller radii where much of the dynamics and chemical evolution occurs.At larger radii the evolution is governed by simple one-body dynamics due to the expected geometric reduction in the strength of the two-body interaction.
The case of a mixture of neutrinos and antineutrinos of different flavors is more interesting.The total number density of a particular flavor of neutrino (or antineutrino) is denoted by n(ν x ), and the net lepton flavor is given by n(ν x ) − n(ν x ).The net lepton number is preserved by the path integral picture with random phases and two-particle exchanges.The conservation of these lepton number conservation conditions (totaling three, one for each flavor) plus total number conservation yield only four conditions, though, when six must be satisfied to define equilibrium, neutrinos and antineutrinos of three flavors each.Assuming random phases renders the quantum problem essentially classical, then it is possible to determine the equilibrium distributions very easily.
In environments where both neutrinos and antineutrinos are present, flavor off-diagonal evolution from processes like ν e νe ⇆ ν µ νµ yields a rapid equilibrium.In such an equilibrium state, the time averaged flux into a particular state must equal the time averaged flux out of a state.Since the magnitude of the Hamiltonian matrix elements are symmetric under time reversal this implies the probability of the product of the densities of neutrinos and anti-neutrinos in each flavor must satisfy n(ν e )n(ν e ) = n(ν µ )n(ν µ ) = n(ν τ )n(ν τ ), where n represents the density of neutrinos of a given flavor.These two additional conditions on the product of neutrino and antineutrino densities, along with the four previous, completely determine the equilibrium condition.
We consider a simply implemented three-neutrinoflavor swap network of 24000 total neutrinos and antineutrinos, each labeled only by their energy, as we expect that the flavor evolution will rapidly tend to isotropy in momentum given the above arguments.We distribute them randomly in energy with a distribution that mimics an energy weighted Boltzmann distribution of the form P(E) ∝ Ee −βE which we show in the top panels for all neutrino and antineutrino species in Fig. 8.We start with initial distributions with fractional n(ν e ) = 1/2 and an average energy of 11.9 MeV, n(ν e ) = 1/6 with an average energy of 15 MeV, and all other species with n(ν x ) = 1/12 and an average energy of 18.2 MeV.
With these initial populations, we perform swaps among the neutrinos and antineutrinos of different energies by randomly selecting two neutrinos from the distribution and swapping their energies.In this process, if a neutrino and antineutrino of the same type are selected, we also allow them to change to a different flavor, with probability 1  3 for each flavor.We perform this swapping process on average 250 times per neutrino in the distribution, and we show the final flavor configuration after the swap network in the two lower panels of Fig. 8.
For the swapped distribution, we find that the average energy of each species is approximately equal and is Ē ≈ 14.5 MeV.The number of electron flavor neutrinos and antineutrinos are reduced to 0.38 and 0.05, respectively by transformation to other flavors with equal products of number of neutrinos times antineutrinos in all flavors.The other flavor populations are increased slightly to 0.14.While the individual flavor populations have been adjusted, the conserved differences have been respected.This rapid approach to equilibrium can significantly impact the dynamics and chemical evolution of supernovae and neutron star mergers.

IV. CONCLUSION
In this work we have provided evidence the generic ν − ν coherent forward scattering Hamiltonian is nonintegrable, and the resultant level spacing statistics behave equivalently to those of a random matrix drawn from a Gaussian orthogonal ensemble.This behavior is generically expected in non-integrable Hamiltonians as there is not an extensive set of conserved quantities which would admit degeneracy in the spectrum.Such Hamiltonians display chaotic behavior, and a large category of initial conditions are expected to thermalize in the sense that one-body expectation values will obtain some equilibrium value from which they stray only transiently, and which can be predicted from an energy diagonal grandcanonical thermal distribution.Artificial impositions of symmetry can categorically change the nature of the coherent forward scattering Hamiltonian.We should only make such simplifying assumptions in the pursuit of understanding if doing so does not qualitatively change the behavior of the system we are studying.This work strongly encourages more careful consideration on how to appropriately take a thermodynamic limit in these systems.
The evidence we present here strongly suggests that late-time one-body expectation values can be obtained a priori from a thermal partition function.While fully diagonalizing the Hamiltonian for large N is prohibitively expensive, even in the lower dimensional block-decimated invariant subspaces, Monte Carlo methods may be feasibly implemented to evaluate the partition function at larger values of N .Such a scheme may provide a method for feasibly determining the late time one-body flavor expectation values in the fast oscillation regime without explicitly solving the many-body Schrödinger equation or numerically diagonalizing the Hamiltonian.
While the time-independent Hamiltonian and initial conditions we consider here are quite simple, we do not expect the neglected effects, including non-forward scattering, spatially resolved initial states or time dependence in the Hamiltonian, to fundamentally increase the coherence in the system.In at least some cases a simple classical picture of equilibrium, as determined by assuming an evolution to an incoherent sum of product states, should provide an accurate approximation to one-body observables.In particular, this path integral picture can be used to define an intriguing equilibrium distribution for systems of neutrinos and antineutrinos of multiple flavors.
Furthermore in the near future quantum computers may provide an avenue for performing the coherent time evolution of the quantum many-body system (see e.g.[44,[79][80][81][82][83] for recent attempts on current generation devices), and may be capable of evaluating expectation values using explicit time evolution for similarly large N , which can be compared to statistical partition functions obtained on either classical or quantum computers.As such, quantum computers may facilitate comparisons between these two predictions.
Finally, as the number of spins included increases, the block-decimated subspaces of the Hamiltonian grow combinatorially large which we should expect to drive the system even closer to the predictions made utilizing the partition function.While we have provided evidence which suggests that the equilibration we have observed here will be obtained when Ĥνν dominates the evolution, a remaining open question is the precise rate at which this thermal equilibrium is achieved.We generally observe equilibration on a timescale which is on the order of ∼ 10µ −1 , however we have not proven that this will occur on such short time scales generically.The examples we have studied do not display sufficient sensitivity to the system size in the equilibration timescale to make a clear determination of the relationship.We leave a more thorough investigation of the approach to equilibrium to future work.

FIG. 1 .
FIG.1.Extracted statistics (black histogram) for the ratio of sequential level spacings of the Hamiltonian for the non-integrable generic parameterization ( Ĥνν , left panel) and for the integrable axially symmetric simplification ( Ĥax νν , right panel) for N = 16 in the j = 2, m = 2 subspace.Overlaid as a red dashed line is the predicted probability distribution for the ratio of sequential level spacings for the non-integrable (GOE, left panel) distribution and the integrable (Poisson, right panel) distribution.

FIG. 3 .
FIG. 3. Magnitude of the off-diagonal matrix elements ⟨Eα| Ô|E β ⟩ with Ô = σ3 on the first neutrino in the energy eigenbasis as a function of the absolute value of the energy difference |Eα − E β | (in units of µ).The left panel shows results for the non-integrable Hamiltonian Ĥνν (blue data) while the right panel shows results for the integrable case Ĥax νν (red data), in both cases we only look at the j = 2 and m = 2 subspace.Also shown is the running average of the matrix elements' magnitude obtained taking the median over a window of 50 matrix elements (yellow curve on the left panel and black curve on the right panel).

FIG. 4 .
FIG.4.In the top panel we show the initial condition (black squares) and the average over the late time oscillations of ⟨σ3,i⟩ obtained from the solution to the many-body Schrödinger equation for each neutrino.The thin black dotted line indicates the conserved value 2⟨ Ĵ3⟩/N .Blue diamonds are the expectation value predicted from the grand-canonical partition function fit from the initial condition, while green circles indicate the energy mixed state distribution (EMSD) average.The lower panel shows the difference (∆⟨σ3,i⟩) between the exact solution and the partition function fit (blue diamond) and the numerical solution to the MB Schrödinger equation and the microcanonical approximation (green circles).The filled red area represents the RMS deviation of the oscillations of the exact solution about its mean.

FIG. 5 .
FIG. 5. Evolution of ⟨σ3⟩ for neutrino i = 1 compared to the prediction utilizing the grand-canonical partition function (black dotted line) at early time.The line and arrow indicate when the time-averaging begins for the evaluation of the data plotted in Fig. 4. Inset: ⟨σ3⟩ as a function of time for every neutrino on the entire considered time domain.

3 FIG. 6 .
FIG.6.Times for which: the average time (TS, indicated by black diamonds) the one-body entropies reach 95% of their maximum; the Loschmidt echo (TL, black squares) reaches its first local minimum; and the average time (TP 3 , indicated by black circles) the the one-body ⟨σ3⟩ reaches its first turning point.
we show two cases.The top panel represents √ N m × |A n | t for the quantum state whose one-body expectation values are shown in Figs. 4 and 5 for all states |n⟩ in the m = 1 subspace, plotted against the corresponding diagonal elements of Ĥνν .The lower panel is similar but for an N = 15 state which initially had purely 10 ν e and 5 ν τ resulting in unit overlap with the m = 2.5 subspace.For the lower panel, spins were chosen ν e and ν τ randomly, the v x/y components of velocity were chosen with random azimuthal angles on the unit sphere, and the v z components were chosen randomly from the uniform interval [0, 1].This figure shows that for both cases, in the respective m subspaces, the time-averaged magnitudes |A n | t ≈ O(1)/

15 FIG. 7 .
FIG.7.Time-averaged magnitudes of state amplitudes (|An|t) for the N = 16 (top) and N = 15 (bottom) cases in the m subspaces (top: m = 1, and bottom: m = 2.5) with the largest total overlap with the total state plotted against the respective diagonal elements of the Hamiltonian in the product state (|n⟩) basis.In each case, the state projected into the specified m subspace has been normalized to unity, and the time-averaged amplitude magnitudes have been scaled by √ Nm where Nm is the total number of states in the specified m subspace.

FIG. 8 .
FIG.8.Initial (top) and Final (bottom) neutrino energy spectra (nν i (E)) for each flavor (denoted i) in a simple calculation assuming random phases from the diagonal parts of the Hamiltonian.The product of the number of electron neutrinos and antineutrinos equilibrates with the other flavors, while the difference is maintained for each flavor.
FIG.2.Diagonal matrix elements ⟨Eα| Ô|Eα⟩ with Ô = σ3 on the first neutrino in the energy eigenbasis as a function of the energy Eα (in units of µ).The left panel shows results for the non-integrable Hamiltonian Ĥνν while the right panel shows results for the integrable case Ĥax νν , in both cases we only look at the j = 2 and m = 2 subspace.