Exact solution of multi-angle quantum many-body collective neutrino flavor oscillations

I study the flavor evolution of a dense neutrino gas by considering vacuum contributions, matter effects and neutrino self-interactions. Assuming a system of two flavors in a uniform matter background, the time evolution of the many-body system in discretized momentum space is computed. The multi-angle neutrino-neutrino interactions are treated exactly and compared to both the single-angle approximation and mean field calculations. The time unit chosen is $\mu_0^{-1}=(\frac{G_F}{2\sqrt{2}V})^{-1}$. The mono-energetic two neutrino beam scenario is solved analytically. I proceed to solve flavor oscillations for mono-energetic cubic lattices and quadratic lattices of two energy levels. In addition I study various configurations of twelve, sixteen, and twenty neutrinos. I find that when all neutrinos are initially of the same flavor, all methods agree. When both flavors are present, I find fast collective oscillations and rapid flavor equilibration develop both for mono-energetic scenarios and when neutrinos have different momentum magnitudes. However, this is not seen in the mean field treatment. The difference can be ascribed to non-negligible flavor polarization correlations being present. Even in dense matter environments I find rapid flavor equilibration from the many-body treatment while the mean field approximation shows no flavor oscillation at all. Entanglement entropy is significant in all such cases. The relevance for supernovae or neutron stars mergers is contingent upon the value of the normalization volume $V$ and the large $N$ dependence of the timescale associated with oscillations. In future work, I intend study this dependence using larger lattices and include anti-neutrinos.


I. INTRODUCTION
Neutrinos are some of the most abundant particles found in nature, produced during the early universe [1][2][3][4][5], from stars like our sun during their lifetime [6,7], and in copious amounts during core collapse supernovae [8][9][10][11][12][13][14][15]. From solar, atmospheric, and reactor experiments we know that they are massive, and the mass eigenstates are different from the weak interaction eigenstates [16][17][18]. As a result, neutrinos oscillate between their three flavors during their propagation in vacuum [19]. For instance, the solar neutrino problem, the discrepancy between the electron neutrino flux emitted from the sun and observed on earth, was a long standing problem [20,21]. The resolution came from the modified dispersion relation in matter which induces resonant flavor conversions commonly referred to as the Mikheyev-Smirnov-Wolfenstein (MSW) effect [22][23][24]. Neutrino-matter interactions induce a self-energy which is the origin of this effect.
In this work I study flavor oscillations under the influence of the self-energy induced by neutrino-neutrino interactions, called collective oscillations. This type of oscillations has been studied extensively in the literature [25][26][27][28][29][30][31][32][33][34][35]. In core-collapse supernovae and merging of neutron stars a very large number of neutrinos are produced [36]. In such situations even the average energies of different flavors are different since only electron neutrinos participate in charged-current weak interactions. * ermalrrapaj@gmail.com Neutrino-neutrino scattering, being between particles of the same type, is of a different nature than neutrinomatter scattering, and gives rise to forward scattering terms in the many-body Hamiltonian which contribute to oscillations [37,38]. These terms are dependent only on the angle between neutrinos and couple neutrinos of different energies making flavor evolution a rather intricate many-body problem. Neutrino transport and flavor oscillations are active areas of research, and no exact treatment has yet been provided. Due to the amount of energy neutrinos carry during a core collapse supernova, they can have a significant impact on the explosion. For instance, the kinetic energy of the material ejected is only about 1% of the neutrino energy [39,40]. If the majority of neutrinos are of electron flavor, they can be re-absorbed by matter and provide the energy required for the explosion. Knowing the flavor evolution is crucial in understanding the role neutrinos play in this environment.
The mean field method was first proposed in [37,38] and has been widely applied subsequently. The time evolution equations in this method become nonlinear and analytical solutions are available only in special cases. However, the method needs to be compared to exact solutions to test its region of validity.
Seminal work on the role of neutrino-neutrino interactions in flavor oscillations has been conducted in [41][42][43]. However, at the time, only a part of the full Hamiltonian was considered, that is a simplified version of neutrino-neutrino interactions, which here I call singleangle approximation, without the presence of the mass term which causes flavor mixing.
Later an algebraic approach was developed which re-tains the many-body nature of the Hamiltonian and where two body operators are not replaced by effective one body counterparts [44]. In this approach, the manybody system is formulated as an SU(2) coherent state path integral of two flavors (SU (3) for three flavors). This approach, as demonstrated in this article, is particularly useful for finding exact solutions and constants of motion. An alternative approximation to the mean field is to replace the angular dependence in neutrinoneutrino interactions by the average angle but keep the two body term; this is commonly referred to as the singleangle approximation. This simplified Hamiltonian can be studied within the Richardson-Gaudin framework [45,46] and one can solve for the eigenvalues [47,48]. Using the single-angle approximation, one can calculate the spectral split for a large number of neutrinos [49]. However, there is no a priori physical motivation for replacing the angular dependence of neutrino-neutrino interactions by an average value. This averaging replaces the complicated angular structure in the two-body neutrino interactions with the squared total angular momentum (in flavor space) of the system, which is not present in the original Hamiltonian. This is the Casimir operator [50] of the SU(2) algebra and commutes with the SU(2) generators, thus it introduces symmetries not present in the original Hamiltonian. From this perspective, the singleangle approximation also needs to be tested for validity. In this work, I solve the many-body multi-angle Hamiltonian and compare to the single-angle approximation and the multi-angle mean field method. As a first step, I assume a simplified scenario of electron neutrinos ν e , and an additional flavor which can be considered as a superposition of tau and muon neutrinos denoted by ν x , in matter at constant density. In addition, I discretize the momentum space, which allows for the treatment of the many-body problem on the lattice. This appears to be one of the first attempts at an exact treatment of neutrino-neutrino interactions on the lattice. The role of anti-neutrinos is postponed for future works. In addition, time evolution for large lattices is computationally expensive for classical computers as the degrees of freedom increase exponentially with the number of lattice points. Hence, I focus on smaller lattices. The Hamiltonian for two neutrino flavors is an ideal candidate for quantum computing given its resemblance to the Ising model [51]. This technology allows one to overcome the limitations of classical computers and explore much larger lattices than currently available. I plan to explore this research direction in future work.
Section II describes the operator algebra required for a treatment of neutrino oscillations in a discretized momentum lattice. Then, section III studies the equations of motion in all the three methods outlined here. Throughout this work, the time unit chosen for plotting results In section IV a system of two neutrino beams is analyzed, and in sections V and VI cubic and two energy level quadratic systems are studied respectively. In section VII a system of twelve neutrinos is studied in vacuum and dense matter. In section VIII two different configurations of sixteen neutrinos in dense matter are studied. I conclude the analysis with a system of twenty neutrinos in section IX. Section X studies the entropy of entanglement of each neutrino in the manybody wave-function. The findings are summarized and conclusions are drawn in section XII.

II. TWO FLAVOR OSCILLATIONS AS A LATTICE OF SU(2) ALGEBRAS
By introducing creation and annihilation operators for one neutrino with three momentum p, the generators of an SU(2) algebra can be written [44,52,53]: Indeed, there is an SU(2) algebra associated with each momentum value, all commuting with each other. To be mathematically rigorous, I discretize the momenta using a box quantization so that I get an SU(2) algebra instead of the usual current algebra. The sum of these operators over all possible values of momenta also generates a global SU(2) algebra. The Hamiltonian is comprised of three contributions: the mass (vacuum) term, matter interactions, and neutrino self-interactions. I start with the first two terms, for which one needs to pick either the mass or the flavor basis to express them. I opted to work in the flavor basis, where B = sin (2θ)x − cos (2θ)ẑ. As Eq. (2) shows, matter effects are equivalent to a modified neutrino oscillation frequency and mixing angle, In Eq. (2), the first integral represents neutrino mixing and the second integral represents neutrino forward scattering off the background matter. The SU(2) algebra operators are represented by Pauli matrices. The electron density N e in the second term of Eq. (2) is inside the momentum integral since electron densities encountered by neutrinos traveling in different directions can be different, but for simplicity here it is assumed to be constant. While this assumption is obviously unacceptable over large distances, here I focus on oscillations that happen over very short ranges for which the approximation should hold. Neutrino-neutrino forward-scattering contributions, which are the basis of collective oscillations, are described by the Hamiltonian where ϑ pq is the angle between neutrino momenta p and q and V is the normalization volume. The (1 − cos ϑ pq ) term in the integral above means that neutrinos traveling in the same direction do not forward scatter off each other. Contributions from collisions are proportional to G 2 F and are sub-leading with respect to forwardscattering.
The total Hamiltonian H ν + H νν describes collective neutrino oscillations and the neutrino masses is the first term in Eq. (2). In writing H νν neutrino masses were set to zero.

III. MANY-BODY HAMILTONIAN IN MOMENTUM LATTICE
The Hamiltonian described in the previous section can be written as a sum of operators at different lattice points,Ĥ In the Heisenberg picture, operators evolve in time due to their commutation relation with the Hamiltonian, Only the mass term is responsible for flavor oscillations in either basis, as the neutrino self interaction term commutes with J k . This leads to the following operator being a constant of motion, This result is expected from previous work in symmetries and invariants of neutrino oscillations. For a complete analysis the reader is referred to [52,53]. One can also explore what happens to a particular lattice point (direction), In this case B · J p is not an invariant of motion due to neutrino-neutrino interactions. As I explore in this work, they lead to fast collective oscillations which can be observed by looking at oscillations of individual lattice points. Equations (5) and (6) are derived in appendix A. As this is one of the first lattice treatments of neutrino oscillations it is important to compare to the approximations commonly used in literature.

Mean field approximation
Neutrino flavor oscillations are commonly studied in the mean field approximation, in which products of two one body operators are approximated by the product of a one body operator and the expectation of the other. While this approximation greatly simplifies calculations, its validity is questionable. The underlying assumption upon which it is based is that two or higher many particle correlations are too small to have a significant impact. In this work I find that when this assumption fails, also the mean field method fails to provide the correct flavor evolution in time. For the moment, let us describe the method. I will return upon this issue later in this article. For the Hamiltonian under consideration, the mean field treatment is given as follows, where P p = J p . The time evolution of the SU(2) algebra generators is, As the expectation value of the polarization P q also evolves in time, an additional equation of motion is needed for completeness and consistency, By solving this equation one studies neutrino oscillations in the mean field approximation. As I show in the next sections, the underlying assumption for this approximation is valid when only one flavor is initially present. I find that when two flavors are initially present, fast collective oscillations are obtained in the exact treatment but do not appear in the mean field method. In the recent literature, linear instability analysis has been employed in collective flavor oscillations and fast conversions have been found [54][55][56][57][58][59]. For these fast conversions to develop, among many other criteria, anti-neutrinos must be present. An additional criterion for the two beam scenario is that the angle between the two neutrino modes must be acute [60]. In my analysis fast oscillations develop without the need for anti-neutrinos (for all lattices considered here) and regardless of the angle in the two beam case, which is quite an interesting difference. I also compare to the single-angle approximation on the lattice. For the mean field method I perform a multi-angle treatment.

IV. ANALYTICAL SOLUTION FOR TWO NEUTRINO BEAMS
If two neutrino beams with the same energy but different directions are considered, one can solve the system exactly by making use of the Dirac basis [61] Any four dimensional matrix M can be decomposed in this orthogonal basis as follows, The Hamiltonian for this system is Since both beams have the same energy, ω = ω 1 = ω 2 . By working in the Heisenberg picture I calculate any operator as a function of time, O(t) = e iHt Oe −iHt . In appendix B, the decomposition of the time evolution operator e −iHt in the Dirac basis is provided. Once the evolution operator is known, the generators (J k = J k 1 I 2 + I 1 J k 2 ) of the SU(2) algebra can be calculated at any moment in time (appendix B), From the time dependence of these operators one can verify that B · J(t) is a constant of motion, As collective oscillations appear by looking at individual points (directions) in the momenta lattice, I focus on point 1 and study the time evolution of the polarization for different initial conditions, and, From these equations I find that the frequency of oscillations is proportional to ω when only one flavor is present and proportional to µ when both flavors are present. Collective oscillations do not depend on the energy scale or the mixing angle, provided α = 0. In addition, the angle between the two neutrinos does not prevent oscillations from developing as long as they are not parallel to each other. To understand the difference between these two initial conditions I study the flavor polarization correlations. Here I define the correlation between two operators as For the two beam scenario I find From these correlations I expect the mean field approximation to match the exact solution for initial conditions where only one flavor is present, and maximal disagreement when the net polarization is 0 (equal amounts of each flavor).

A. Single-angle approximation
The average angle for the two beam scenario is cos ϑ = (1 + cos ϑ 12 )/2. The respective time evolution of the expectation value of J Z 1 I 2 is identical to the multi-angle results for initial conditions of only electron or x flavor. When both flavors are present the results are different, Correlations vanish for a single flavor and differ for mixed flavor initial conditions, The attentive reader might be puzzled as to why I am comparing multi and single-angle calculations for a system of two particles where indeed there is only one angle present between them. As shown in append B, the single-angle approximation allows for neutrinos to interact with themselves, differently from the original Hamiltonian. Even for a small system of two particles this causes a difference with exact results. This difference seems mainly quantitative rather than qualitative at this stage. However, as the number of particles increases, so does the number of angles between them, and one needs to compare with multi-angle results.

B. Mean field approximation
The total polarization and difference of polarizations are defined as This simplifies calculations as equations of motion partially decouple, The first equation has the analytical solution, where P 0 = P tot (t = 0). More details can be found in appendix C. The second equation is solved numerically for generic initial conditions. However, here I am interested in two specific cases: two neutrino beams of one flavor (|ν e ν e ) and one neutrino of each flavor (|ν e ν x ).

Case: |νeνe
Initially both neutrino beams are of electron flavor, Since the time evolution equations are of first order, once one knows the initial conditions, the solutions can be found, The expectation values for individual polarizations are, By trigonometric identities one can prove that Eqs. (14) and (9) are identical. In this case all methods agree, there are no correlations present, and the frequency of oscillations depends on ω but not on µ.

Case: |νeνx
If the initial wave-function has equal amounts of electron neutrinos ν e and x neutrinos ν x , The total polarization is a constant of motion and the difference of polarizations depends only on the mass term, Equation (15) is identical to Eq. (13) in the previous section and, thus, has the same solution. The expectation values for individual polarizations are, The time evolution for individual modes depends only on the mass term, in sharp contrast with the result from the exact solution. For initial conditions without net polarization the mean field method shows no collective oscillations, which is quite interesting, as from the previous section, this is the scenario where correlations are present.

C. Numerical results for |νeνx
To simplify the numerical implementation of the equations in the previous section, µ 0 = √ 2 G F V is set as the energy (inverse time) unit. I focus on mode 1, and assume the initial wave-function is |Ψ 0 = |ν e1 ν x2 . and x, sin θ vacuum =0.297 [62]. By varying the two couplings and the angle between the two modes one can see that fast collective flavor oscillations develop in vacuum. As Eq. (10) shows, the mass term ω p is responsible for the small oscillation frequency, and the high frequency is linearly proportional to µ. In other words, these are frequency modulated oscillations where ω p plays the role of the original frequency and µ is the modulating frequency. Oscillations develop regardless of the sign of cos (α). However, the mean field calculation fails to find any collective oscillatory behavior and only displays the usual vacuum oscillation with frequency ω p .
Next I turn my attention to flavor oscillations in the presence of matter. While a detailed study is beyond the scope of this work, I focus on the adiabatic approximation (constant matter density). In other words, I study length scales for which matter density does not vary drastically. As Eq. (3) shows, the effective mixing angle in dense media is much smaller than in vacuum and the effective mass coupling increases. In Fig. 2 I used α = 0.01 as an illustration. Collective oscillations develop for both multi and single-angle calculations while the mean field approximation shows no oscillations at all. It seems that as long as there is some mixing, no matter how small the mixing angle, fast oscillations develop and the frequency of oscillations in this case depends only on µ. The explanation for the mean field result can be found by analyzing Eq. (15). If the mixing angle is very small and there is no net polarization initially, P tot (0) = 0 and P diff (0) =ẑ.
Then, the rate of change of the polarization difference is This is a rather interesting finding as collective oscillations in dense media could have quite dramatic implications for core collapse supernovae and neutron star binary mergers. One can perform an order of magnitude estimate of the ratio between the neutrino mean free path from charged reactions and for collective oscillations. To estimate the ratio of these two length scales, I assume the normalizing volume in the two body term to be V ∼ n −1 ν , This result seems to suggest that once collective oscillation develop they will effectively equilibrate neutrino flavors between collisions with matter. However, the result depends on the numerical value of V , which I discuss in section XI. I also consider the case for which the initial wave-function is only electron flavor. From Eqs. (9) and (12) the single-angle approximation gives the same result as the multi-angle treatment and they both show ordinary vacuum oscillations. In addition, no correlations are present, and the mean field approximation result in Eq. (14) agrees with the exact solution. In this section I consider a three dimensional cube as depicted in Fig. 3. The basis used in this case is a direct product of eight Pauli matrices which allows for the exact computation of the time evolution operator. I study os- cillations in matter with an equal number of electron and x flavor in the initial wave-function. Figure 4 shows the correlation of all points and the polarizations for two lattice points with initial flavors ν e and ν x respectively for a mixing angle α = 0.01. The mean field shows no flavor evolution, while the many-body methods depict rapid flavor equilibration. In the single-angle approximation, as correlations change between maximal and minimal values, so does the flavor mixing. The multi-angle calculation shows strong correlations, and the flavor polarization oscillates close to zero as a result. The presence of correlations explains the difference between exact calculations and mean field similarly to the two beam case. As before, correlations disappear when all neutrinos start with the same flavor and all methods agree. I verified numerically that B·J is a constant of motion. Flavor evolution in the vacuum with no net initial polarization and alternatively with only one flavor initially is provided in appendix D. in medium with ω/2 = µ = µ0 and α = 0.01.

Initial conditions with partial polarization
In the previous sections I have studied case when all neutrinos start with one flavor, or both flavors are present in the same amount. However, an intermediate scenario is needed where there is some flavor polarization in the initial wave-function, but not all neutrinos are of the same type. As an illustration, I performed the flavor evolution for the cubic lattice with an initial wavefunction of six electron neutrinos and two x neutrinos, |Ψ 0 = |ν (6) e ν (2) x . Figure 5 displays the flavor evolution for two lattice points which start oscillations with ν e and ν x respectively. Collective oscillations still develop but are not as pronounced as in the previous cases. Similar to previous results, the mean field method does not show any collective behavior. In addition, both the single-angle and mean field approximations display oscillations with a contrary phase with respect to the exact solution when it comes to the flavor of lower concentration. In the multi-angle lattice treatment, the x flavor quickly synchronizes in phase with the electron flavor, but this does not happen for the two other methods.

VI. TWO ENERGY LEVELS: QUADRATIC LATTICE
So far mono-energetic flavor oscillations have been considered. Due to their interactions with matter during a core collapse supernova, electron neutrinos decouple in lower density regions with respect to other flavors and the resulting energy spectrum is lower [40,63,64]. To account for this difference in spectrum, in this section I work with two quadratic lattices of different momenta values as shown in Fig. 6. I pick a momentum magnitude ratio of three as a conservative representation of the mean energy ratio found in numerical simulations [40,63]. I denote by ω p the frequency caused by the lower momentum value. B. Initial Conditions with partial polarization I verified that if the initial wave-function contains only one neutrino flavor, no collective oscillations occur and the mean field agrees with the lattice calculation. In addition, correlations vanish as well.
Next I consider a mixture of 75% electron flavor, |Ψ 0 = |ν (6) e ν (2) x like in the case of the cubic lattice. From Fig. 8 one can see that collective oscillations develop and from fig. 9 that correlations are present. Lattice points 4 and 8 show flavor oscillations for the larger lattice with ν x and ν e flavors initially; lattice point 1 is a representative from the smaller lattice. The amplitude of collective oscillations is less pronounced and it oscillates close to zero. The mean field still differs from the exact solution.
As a final check, I study oscilla- tions in a dense medium in Fig. 10. I doubled the range of the time axis to observe the pattern of oscillations which is rather peculiar. For the smaller lattice which initially is electron flavor, the mean field shows no oscillations, while the multi and single-angle show collective oscillations develop, but the polarization never becomes negative. For the larger lattice, the mean field displays suppressed oscillations with a low frequency. The manybody methods shows collective oscillations with a high modulation frequency. Figure 11 confirms the presence of correlations among neutrinos of different energies even in a dense medium.   The results depicted in these figures agree with previous section: there is rapid flavor equilibration in the many-body methods, particularly in the multi-angle case, while the mean field shows no such behavior. Not surprisingly, correlations are significant when this happens. Especially in a dense matter medium, the qualitative difference is quite visible as the mean field shows no oscillations at all. The neutrinos remain in the initial flavor. However, the multi-angle shows rapid flavor equilibration, more so for the flavor initially present in smaller amount. The single-angle shows a result in between the multi-angle and mean field.

VIII. SYSTEMS OF SIXTEEN NEUTRINOS
In this section I continue the many-body study of the neutrino flavor oscillations by considering two different configurations of sixteen neutrinos.
At first I consider two cubic lattices of electron and x flavor respectively, with the x flavor having momentum magnitude three times higher than electron neutrinos. In Fig. 16 I plot the flavor polarization and correlation in a dense medium with a reduced mixing angle  The second configuration has three times more electron neutrinos with a momentum magnitude three times lower.
The flavor polarization and correlations depicted in Fig. 18 show the same trend; with oscillations not present in the mean field approximation and flavor equilibration in the many body calculation.

IX. SYSTEM OF TWENTY NEUTRINOS
I conclude the many-body study with a system of twenty neutrinos. The configuration is similar to fig. VII with eight neutrinos in the smaller cubic lattice and three quadratic lattices of x flavor neutrinos. Similarly to previous sections, the ratio of the momenta magnitude is three. Figure 19 shows the flavor polarization for two neutrinos in the lattice and overall correlation as with α = 0.01. The trend found in previous sections is also present here. One would expect the mean field to approach the many-body results as the particle number is increasing, but this does not happen here. However, given the relatively small number of particles considered in this work, no statements can be made regarding large particle number behavior. Larger lattices will be studied in future work.

X. ENTANGLEMENT ENTROPY
From the configurations studied in this work a common feature emerges; whenever polarization correlations are significantly present, the mean field approximation fails to capture the time evolution described by the manybody method. To have a better understanding, in this section I study the entanglement entropy of each neutrino in the many-body wave-function versus the rest. Since the initial wave-functions considered in this work are eigenstates of flavor, they are pure states. And, as they time evolve due to the Hamiltonian, they remain pure, although not eigenstates of flavor anymore. For such wave-functions, the Von-Newman entanglement entropy serves as a measure of entanglement between the neutrinos. For a given density matrix it is defined as, S = − Tr (ρ ln ρ) , ρ = |Ψ Ψ| (16) Here I focus on a bipartition of the wave-function with one neutrino in one subspace, and the rest of the neu- trinos in the other. For an in depth description of entanglement and its various measures used in quantum computing, the reader is referred to [65][66][67][68]. The bipartition chosen allows one to study whether one neutrino can be factored out from the rest of the wave-function. If this is possible, the single particle mean field should describe the many-particle system. The wave-function can be decomposed based on this partition, The label i denotes some basis of orthonormal wavefunctions of N − 1 neutrinos. The reduced density matrix is obtained by tracing out the N − 1 neutrinos, ρ ν = Tr Ψ i ν (N −1) ρ Ψ ν (N ) and the respective entanglement entropy is S ν = −Tr (ρ ν ln ρ ν ). This quantity is zero when there is no entanglement and ln(2) when the neutrino is maximally entangled with the other neutrinos in the wave-function. Thus, if S ν = 0 one expects a single particle mean field description to be appropriate. As the entropy increases, the mean field result should get further away from the exact calculation. In appendix D in fig. 24 I find both correlations and entanglement entropy vanish when the mean field approximation agrees with the many-body result. This happens when only one flavor is initially present.
When both flavors are initially present, correlations and entanglement develop. For two neutrino beams, one For all lattices studied in this work, when both flavors are initially present, I find that neutrinos get rapidly entangled. Figure 21 displays the entanglement entropy for systems of eight, twelve, sixteen and twenty particles as function of time for initial electron and x flavor respectively. In [41] the authors find that entanglement does not develop when only the single-angle two body term in the Hamiltonian is considered. Later, in [42], entanglement was found for that Hamiltonian when systems of fourteen particles were studied. In the case at hand, I am studying the full multi-angle Hamiltonian and I find entanglement when both flavors are initially present. Due to the fact that the Hamiltonian dimensions increase exponentially with the particle number, it is very computationally expensive to study much larger systems. While the main focus of this work is the multi-angle Hamiltonian, it is worth pondering what happens for larger and larger systems. A simple assumption would be that correlations would decrease as the number of particles increases, thus entanglement would decrease as well. Unfortunately, if this question is posed for the transverse Ising spin chain the answer contradicts this assumption. In [69] the authors find that the entanglement entropy increases with time until it reaches a maximal value. This value of the entanglement depends on the initial state and does not decrease and the particle number increases. A conclusive statement about the multi-angle Hamiltonian can be made once larger systems are studied, which is to pursued in the future. In the meantime, I will focus on the effects of the two-body operator only, and for simplification, I will consider the single-angle approximation and compare it to the mean field result in the next section.

A. Many-body treatment
In this section, in order to study infinitely large systems, I will ignore the mass term in the Hamiltonian, and study the system in the single-angle approximation. This is the Hamiltonian studied in [43], where, Here I follow the same steps as in [43], and proceed to calculate the polarization and entanglement entropy as functions of time. The initial wave-function chosen is Based on [43] one can express the wave-function at a given time as a direct product of a single neutrino state and the rest, where η(N, J) is the Clebsch-Gordan coefficient [70], The reduced density matrix is obtained by tracing over the subspace of 2N − 1 neutrinos, The function A(N, µ, t) is given as follows, where 2 F 1 (a, b; c; z) is the Gauss hypergeometric function. The expression in front of |ν e ν e | is the so called survival probability for the neutrino to remain in the initial flavor [43], while the expression in front of |ν x ν x | is the probability for it to change flavor. In the full Hamiltonian B · J is a constant of motion while for this Hamiltonian J z is. This means that the coefficients for |ν e ν x | and |ν x ν e | are zero in this case by symmetry, but not in the full Hamiltonian. The entanglement entropy is, In addition, one can also study the polarization of individual neutrinos from the expressions above, The complete derivation of |Ψ(t) can be found in appendix E. The derivation of the survival probability has already been performed in [43]. The authors did not study entanglement entropy or compare the polarization from the exact solution with the mean field in that article. Their focus was on the time scale associated with flavor equilibration, and found it to be τ −1 ∼ µ 0 √ N . As µ 0 depends on the normalization volume V , the authors took V ∼ cm 3 , and concluded that the time scale was too large to be significant for supernovae (τ ∼ 10 22 s). However, one could argue that this quantity should be related to the neutrino density, V −1 ∼ n ν . If nν n0 ∼ (O)(10 −10 ), where n 0 = 0.16 fm −3 is the nuclear saturation density which is reached at the core of proto-neutron stars, the time scale becomes relevant (τ ∼ 10 −6 s). As the goal of this work is the mathematical solution for the time evolution due to the multi-angle Hamiltonian, I refrain from making any conclusive statements as to which choice is appropriate. I proceed in the next subsection to compare to the mean field result for the same Hamiltonian and initial wave-function.

B. Mean Field method
The equations of motion are derived from eq. 8 by setting ω p = 0, In fig. 22 I plot the polarization and the entanglement entropy as functions of time for the system of 1000 neutrinos. The many body result is provided by the analytical expression if eq. 22 and the mean field is solved numerically. The qualitative difference between the many-body result and the mean field is quite striking. The mean field shows no oscillations at all, while the exact result shows rapid flavor equilibration. This is agreement with all the configurations studied so far. The lack of oscillations in the vacuum predicted by the mean field can be explained by comparing this system with the two particle system in section IV B. All the 500 initial electron neutrinos are equivalent for the given Hamiltonian, and the same can be said for the x flavor ones. The system of equations in 22, then, can be written in terms of the sum of the polarizations, and the difference between electron and x flavor. From section IV B, the equations of motion for the 1000 neutrino system can be written as, The mass term, responsible for the time evolution of these two quantities is zero, so they are both constants of motion. No matter how large the system, if the initial wavefunction has equal numbers of both flavors, the mean field will show no oscillation, On the other hand, lim t→∞ F ( √ N µt) = 0, so lim t→∞ P z νe (t) = lim t→∞ P z νx (t) = 0.
Actually, the asymptotic value is reached very quickly as fig. 22 shows. Not surprisingly, neutrinos get rapidly entangled with each other and this quantity reaches its maximum rather quickly as shown in the same plot. This result agrees qualitatively with what I have found in the previous sections.

XII. CONCLUSION
Neutrino flavor oscillations are important for the thermodynamic evolution of core collapse supernovae and neutron star mergers, and also impact the nucleosynthesis of heavy elements in these environments. Due to the complicated nature of this many-body system, the mean field approximation is widely used in numerical simulations. However, it needs to be compared to an exact many-body treatment in order to assess its region of validity. By considering uniform matter and discretized momentum space, I am able to solve exactly the time evolution of flavor oscillations. I provide analytical results for the two neutrino beam scenario and numerical results for cubic and quadratic lattices of eight, twelve, sixteen and twenty particles. I find that when both flavors are present, fast collective oscillations develop and there is flavor equilibration but these are not observed in the mean field treatment. This qualitative difference can be ascribed to correlations developing among neutrinos for the mixed flavor scenario as neutrinos rapidly entangle in time units When only one flavor is present, correlations vanish, there is no entanglement, oscillations vanish, and the mean field agrees with the exact solution. While the single-angle approximation on the lattice does show fast collective oscillations, there are differences from the multi-angle result when the two neutrino flavors are initially present but not in equal amounts. To confirm my findings I calculate the flavor polarization and entanglement entropy in the large N limit by studying only the two-body contribution in the single-angle approximation based on [43].
The implications for supernovae might be quite interesting as the neutrino mean free path due to weak reactions could be larger than the mean free path for collective oscillations depending on the choice of the normalization volume V . To give a more definitive answer, anti-neutrinos must be included and larger lattices need to be considered. I plan to follow up on this in future work.
Similarly, for the individual lattice point I find

(B5)
Given this decomposition, one can compute any operator as function of time by working in the Heisenberg picture, O(t) = e iHt Oe −iHt . In the main article I have provided J x (t) and J z (t). Here I also provide J y (t), I also provide the flavor polarization for each lattice point, been provided in the main article, here I provide the rest, ν e1 , ν e2 |I 1 J Z 2 |ν e1 , ν e2 (t) = − ν x1 , ν x2 |I 1 J Z 2 |ν x1 , ν x2 (t) = 1 2 cos(2α) sin 2 (ω/2 t) + cos 2 (ω/2 t) , A similar procedure is employed for the single-angle approximation, but for every pair interaction I use the same angle, η = (1 + cos ϑ 12 )/2. The qualitative difference between the single and multi-angle results lies in the fact that in the multi-angle case, a particle can not interact with itself, while the single-angle approximation allows for this to happen: For the two particle system this means an overall constant difference between the two Hamiltonians as there is only one angle between the two neutrino beams. However, for a greater number of particles, there are many angles present, and the difference between the two Hamiltonians is not a simple constant. For initial conditions with partial polarization, as shown in sections V and VI, there are qualitative differences in the respective flavor evolution.

Electron flavor initial conditions
In the two beam case I found complete agreement among the exact solution and various approximations if only one flavor is present initially. As Fig. 24 shows, the same happens for a cubic lattice and only the usual vacuum oscillations are present. There are no correlations or entanglement entropy in agreement with the result for the two neutrino beams.

Appendix E: Infinitely large systems
The derivation here is based on the work in [43] and the reader is referred to that article for an in depth set up of the system at hand. The initial wave-function can be decomposed by means of Clebsch-Gordan coefficients in the basis of total angular momentum with zero z-projection, Since the goal is to understand what happens to a particular neutrino, here I focus on the first one, which happens to be initially of electron flavor.
This is a bipartition of a system of 2N into a subsystem of 1 and a subsystem of 2N − 1 neutrinos. However, exactly the same derivation can be performed for any other neutrino, for instance the last one, which is initially of x flavor. In addition, each of the wave-functions in eq. E4 can be decomposed as direct products,