Entanglement and Many-Body effects in Collective Neutrino Oscillations

Collective neutrino oscillations play a crucial role in transporting lepton flavor in astrophysical settings, such as supernovae, where the neutrino density is large. In this regime, neutrino-neutrino interactions are important and simulations in the mean-field approximation show evidence for collective oscillations occurring at time scales much larger than those associated with vacuum oscillations. In this work, we study the out-of-equilibrium dynamics of a corresponding spin model using Matrix Product States and show how collective bipolar oscillations can be triggered by many-body correlations if appropriate initial conditions are present. We find entanglement entropies scaling at most logarithmically in the system size suggesting that classical tensor network methods could be efficient in describing collective neutrino dynamics more generally. These observation provide a clear path forward, not only to increase the accuracy of current simulations, but also to elucidate the mechanism behind collective flavor oscillations without resorting to the mean field approximation.

Neutrinos play a pivotal role in extreme astrophysical events like core-collapse supernovae and neutron star mergers, where they are responsible for both reinvigorating a stalled shock-wave and controlling the conditions for nucleosynthesis in the ejected material [1][2][3][4]. In these environments with a large neutrino density ρ ν , neutrino flavor evolution is substantially modified by neutrino-neutrino scattering processes which can lead to self-sustained collective flavor oscillations [5][6][7][8][9][10][11][12][13][14]. Since neutrinos in supernovae are emitted with fluxes and spectra that are strongly flavor dependent [2], the presence of collective flavor oscillations could then lead to important effects [15][16][17]. In some situations however (see eg. [18,19]) these could appear only at too large radii to modify the explosion mechanism in supernovae.
For large neutrino densities near the protoneutronstar, flavor oscillations can occur at frequencies of the order of the neutrino-neutrino interaction energy µ = √ 2G F ρ ν instead of the much smaller vacuum oscillation frequency ω = ∆m 2 /(2E), a phenomenon known as fast flavor conversion [20][21][22][23] (see [24] for a review). In these expressions G F is Fermi's constant, E the neutrino energy and ∆m 2 = m 2 2 − m 2 1 the squared mass gap in the two flavor approximation used here.
The standard approach to study these collective phenomena is to consider only the leading order forwardscattering neutrino-neutrino interactions where only flavor can be exchanged among neutrinos (for the full description see eg. [25,26]). In this scheme, a system with N neutrinos is represented as a collection of SU (2) flavor isospins governed by the Hamiltonian [27] [28] with σ i = (σ x i , σ y i , σ z i ) the vector of Pauli matrices acting on the i-th neutrino. In the expression above, E k are the neutrino energies and, in the flavor basis {|↑ , |↓ } = {|ν e , |ν x }, the orientation of the (normalized) vector B = (sin(2θ), 0, − cos(2θ)) is related to the mixing angle θ. The geometry of the problem is encoded in the two-body coupling matrix J ij = (1 −p i ·p k ) witĥ p k = p k /| p k | and p k the momentum of the k-th neutrino. Finally, note that this Hamiltonian is obtained under the assumptions of infinite interaction times and that neutrinos have infinite plane-wave solutions, conditions that might not always be appropriate (see eg. [29]). This Hamiltonian model is the base of all the past studies of (two-flavor) collective neutrino oscillations in anisotropic but homogeneous settings (see [27] for an explicit derivation of the traditional evolution equations in the mean field approximations from Eq. (1)) . The inclusion of spatial inhomogeneity in the neutrino density distribution, important to describe more realistically the environment encountered in a supernova explosion, will be the subject of future work. The focus here is improving the understanding of the simpler homogeneous case as a fundamental stepping stone for more realistic simulations of flavor dynamics in astrophysical environments.
Collective neutrino oscillations are understood as being caused by unstable modes in the mean field dynamics generated by the Hamiltonian in Eq. (1) which can amplify initially small flavor perturbations exponentially fast [23,30,31]. The role of entanglement and quantum effects in the out-of-equilibrium dynamics [32] of neutrinos has received renewed interest recently [33,34] but is still not well understood, with seemingly conflicting results in the past predicting either a vanishingly small contribution in the large system size limit [35,36] or substantial flavor evolution over time scales τ F = µ −1 log(N ), which can remain relevant for large systems [30,37]. We will refer to oscillations on the time scale τ F to be "fast" in contrast to "slow" oscillations occurring for τ L = µ −1 √ N . This convention is different from the more common definition in the literature where "fast" and "slow" oscillations are associated with time scales ≈ µ −1 and ≈ √ µω respectively, for any system size.
In this work, we study the model Eq. (1) for large systems with up to N = 128 neutrino amplitudes in a simplifying limit and show how coherent flavor oscillations arXiv:2102.10188v2 [hep-ph] 18 Nov 2021 at the fast scale τ F can occur even without vacuum mixing. These first-principle simulations, obtained without resorting to the mean-field approximation or exploiting the symmetries present in the model, are made possible by using a Matrix Product State (MPS) [46] representation of the many-body neutrino state. This technique (described in detail in Sec. II below) allows for efficient simulation of the real-time dynamics in situations where the entanglement entropy (which provides a measure of quantum correlations) is relatively small and growing only as the logarithm of the system size. In the context of collective neutrino oscillations, this is potentially very beneficial for two main reasons: first, MPS calculations are a promising route to perform full many-body simulations of neutrino dynamics in medium to large system sizes in order to understand the impact of commonly used simplifications like the mean-field approximation, second, and possibly more importantly, the direct connection between the simulation efficiency of MPS calculations and the amount of quantum correlations generated by the neutrino dynamics can be used as a direct probe of the role of entanglement and its time evolution in the instabilities that lead to the presence of collective modes in neutrino oscillations. An example of the possible insights that can be gained from the entanglement properties of a neutrino system is the connection between dynamical phase transitions and the presence of collective oscillations discussed in Ref. [57]. Those result were made possible by using the MPS techniques introduced in the present work and have the potential to complement criteria like linear stability analysis [23] in the search of conditions for collective modes to be active in complex astrophysical scenarios like supernova explosions.
The rest of this paper is organized as follows. In Sec. I we introduce the details of the model and the initial conditions used in this work. Sec. II provides an introduction to MPS simulations and the techniques used to approximate the time evolution under the neutrino Hamiltonian in Eq. (1). The results of the simulations and described in Sec. III and in Sec. IV we provide a summary and future perspectives. The Appendices contain a direct comparison with the mean field approach and a discussion of the numerical accuracy of our implementation.

I. NEUTRINO SPIN MODEL
We consider a simple model, in the mass basis, in which the neutrinos are divided into two equally sized "beams" (A and B) and initialized as the product state |↑ . This corresponds to choosing all the neutrinos in the A beam to be |ν 2 and all the neutrinos in the B beam to be |ν 1 . We will consider the simple case where the energies of the neutrinos in beam A are all equal to E A and similarly for the neutrinos in beam B. Finally, we will neglect the geometric information encoded in the coupling matrix J ij and take the single angle approximation by setting J ij = 1 for all neutrinos [27,31] (see also [38]). Note that these are not necessarily realistic conditions, as neutrinos are produced in flavor states, but nevertheless contain many of the features relevant for the astrophysical problem. As we will see, the choice of using the mass basis is also useful to better expose differences with the mean-field treatment.
The resulting neutrino Hamiltonian, similar to the one considered in Ref. [33], can be written as follows where A and B are sets of indices for the neutrinos in the A and B beam respectively while we set the one body energies as ω k = ∆m 2 /(2E k ). Since our initial state is completely polarized along the z-axis, the mean field approximation for the dynamics generated by H predicts no flavor evolution at all: any time evolution of the flavor polarization is completely driven by many-body correlations. The Hamiltonian in Eq. (2) is integrable and a complete solution could be obtained using the Bethe ansatz [27]. This was used in Ref. [33] to study a timedependent generalization of the model for small systems. At this point it is convenient to introduce the total spin operator J = i σ i /2 and correspondingly J A and J B for the neutrinos in the two beams. The full Hamiltonian in Eq. (2) commutes with both the z component of the total spin, J z , and its magnitude, J 2 z . As the initial state |Ψ 0 is an eigenstate of J z , we can rewrite the relevant Hamiltonian for our setup as follows where we dropped irrelevant constant factors and introduced a single parameter δ ω = (ω A − ω B )/2 for the energy difference between the two beams. This model is a variant of the Lipkin-Meshov-Glick (LMG) model [39] explored extensively in the past [40][41][42][43][44]. The model with no vacuum term merits special attention as it is diagonal in the angular momentum basis and exact results for the flavor dynamics can be obtained as shown in Refs. [35,36]. In this model, collective flavor oscillation are shown to occur at time scales much longer than τ F with the scaling τ L = µ −1 √ N expected from incoherent scattering from the background neutrinos [35,45]. Due to this strong scaling with system size, when N 1 the flavor dynamics is essentially frozen and the mean field prediction of no evolution eventually becomes exact. This has been taken as an indication that the mean field description is appropriate for large systems and that entanglement does not play any major role [35,45].
As a way to track flavor oscillations, we compute the survival probability p(t) = (1 − σ z 1 (t) )/2 for a neutrino (we take the first) to be found in it initial flavor state (in this case |↓ ). By using the exchange symmetry for neutrinos within each beam of both the Hamiltonian and the initial state |Ψ 0 , together with energy conservation, we can write the expectation value of the total spin as When δ ω /µ < 0 the total spin, starting already at a relatively small value, can only decrease further during time evolution. Since J 2 is positive semi-definite, this introduces a lower bound p(t) ≥ 1 − µ/(2N δ ω ) and the neutrinos remain frozen in the initial state in the thermodynamic limit. In the opposite limit, for δ ω /µ ≥ 1/4 flavor oscillations become parametrically small as p(t) ≥ 1 − µ/(4δ ω ) but remain finite in the N 1 limit. In the intermediate regime 0 < δ ω /µ 1/4 we instead expect to find large amplitude flavor oscillations.
We solve for the real-time evolution of the system using a Matrix Product State(MPS) ansatz for the many-body wave function [46][47][48]. The approach has a controllable error and it's computationally efficient for systems with low levels of bipartite entanglement [46,49] and is described in detail in the next Section.

II. MATRIX PRODUCT STATE SIMULATION
In this section we review the Matrix Product State (MPS) technique used to generate the results in this work. Additional information can be found in the original paper [46] and in more recent reviews [48,50]. In order to implement long-range interactions efficiently with a MPS, we use the same strategy introduced recently for circuit based simulations in [51], a brief explanation of the construction is presented in the last subsection. In the last section we also compare directly the meanfield approximation to an MPS-based simulation with unit bond dimension.

A. MPS Representation
Consider a general quantum state |Ψ of a system with N spin-1/2 particle, it's wavefunction can be expressed as a N -component tensor as with |σ 1 · · · σ N one of the 2 N basis states spanning the total Hilbert space of the system. A Matrix Product State is an exact representation of the N -component tensor Ψ σ1···σ N as a product of N , smaller, 3-component tensors A σ k α,β . A constructive procedure to find this representation, starting from the original tensor Ψ σ1···σ N , can be obtained by applying N times a Singular Value Decomposition (SVD) of appropriately reshaped tensors. This decomposition allows to represent an arbitrary rectangular matrix M of dimension D A ×D B as the following product of three matrices where U is a D A × min(D A , D B ) matrix with orthonormal columns, Σ is a diagonal matrix of dimension min(D A , D B ) × min(D A , D B ) with either zero or positive entries and V † is a min(D A , D B ) × D B matrix with orthonormal rows. An important property of the SVD, which we will use repeatedly in the derivation below, is that the approximation of the original matrix M , with rank r ≤ min(D A , D B ), using a second matrix M with rank r < r can be done optimally, in the Frobenius norm, by truncating the sum in Eq. (6) so that we keep only the r largest singular values Σ kk . We will assume that the SVD is performed with the convention that singular values in Σ are ordered as Σ 11 ≥ · · · ≥ Σ N N . With this convention, the optimal approximation corresponds to truncating the sum as Starting from the first spin on the left, we first reshape Ψ σ1···σ N into the rectangular matrix Ψ σ1(σ2···σ N ) of dimension 2 × 2 N −1 . An explicit SVD of this matrix using Eq. (6) will then produce the following decomposition with r 1 ≤ 2 the rank of the original rectangular matrix and in the second line we have reshaped the product of Σ and V † to a vector W . In the last line, we have further reshaped the matrix U into two row vectors, with components A [1]σ1 α1 = U σ1α1 , and the vector W to a 2r 1 × 2 N −2 matrix Ψ [1] (α1σ2)(σ3···σ N ) . At this point, we can repeat the procedure by applying a SVD to this last matrix and proceed from left to right until we have reached the last spin of the chain. The resulting decomposition becomes where in the last line we suppressed the lower indices and replaced them by matrix multiplications. Using this construction, a general state |Ψ can be written as a Matrix Product State. Note that the rank of the SVD in the center of the chain could be as large as 2 N/2−1 (assuming N even) and storage requirements for representing a general state will then scale exponentially with the system size. As we will see soon, the maximum rank r max = max k r k needed to describe accurately a given state is related to the maximum bipartite entanglement entropy in the system. In order to see this, it is first useful to introduce a different construction for the MPS. Consider a similar iterative decomposition procedure for the tensor Ψ σ1···σ N as the one described above, but now starting from the right-most site. The first step leads to where we have introduced the two column vectors B β N σ N and the residual matrix has maximum dimensions given by 2 N −2 × 2r N . Proceeding with the same procedure until the left edge, we find another matrix product state with amplitudes The two representations are related by a gauge transformation [48,50], and they are characterized by different normalization properties. In particular, we have For this reason, the state resulting from the decomposition in Eq. (9) in terms of A tensors is called a leftcanonical MPS, while the construction in terms of B tensors in Eq. (11) is called right-canonical [48,50].
In order to see more clearly the connection between entanglement and the maximum rank in the decomposition, also called bond dimension, it is useful to introduce yet another canonical form. This is obtained by starting with a decomposition from the left up to a chosen site k, at this point we start with a decomposition from the right and we stop when we reach site k + 1. We can write the result of this decomposition as with the diagonal matrix Σ [k] containing the singular values at the bond between the k and the k + 1 sites. The site k is called the orthogonality center, since tensors to the left are left-normalized while tensors on the right are right-normalized. This property allows us to write the state in this representation as where we have defined two sets of states and (16) Thanks to the normalization properties of the A and B tensors from Eq. (12), we see that these states form orthonormal sets of states for the left and right Hilbert spaces respectively. We can then interpret Eq. (14) as the Schmidt decomposition of the state |Ψ on the leftright bipartition. The number of states in the sum can then be identified with the Schmidt rank, which is a natural measure of entanglement between states on the left partition and states on the right one.
In particular, the Von-Neumann entanglement entropy for the left-right bipartition at site k of the state |Ψ can be expressed as a function of the singular values as The Schmidt rank r k gives the number of non-zero singular values so that the entanglement is bounded by We can now see that MPS with small bond dimensions correspond to quantum states with little entanglement [46]. Before moving on to the time-dependent simulation, it's important to mention that the mixedcanonical form in Eq. (13) is also particularly convenient to find more efficient approximations to an initial MPS.
Given an MPS |Ψ with bond dimension r k at the bond between site k and k + 1, we can find another MPS |Ψ with lower bond dimension r k < r k so that the error in 2-norm is bounded from above by [52] |Ψ − |Ψ Thanks to the bi-orthogonality of the Schmidt decomposition Eq. (14), we can obtain |Ψ from |Ψ by keeping only the r k largest singular values in the sum. If the singular values decrease sufficiently fast, a good approximation to the original state could be obtained with an MPS with a much smaller bond dimension [46].

B. Time evolution of a MPS
We now turn to explain how one can simulate the timeevolution for the general neutrino Hamiltonian in Eq.(1) of the main text. The basic operation we will use repeatedly in the following is the application of an SU (4) operation on a pair of nearest neighbor indices (k, k + 1). A generic 4 × 4 unitary U of this type can be represented as a 4-component tensor G σ1σ2 with each index having dimension equal to 2. When we contract this tensor with the MPS, only the tensors around the bond (k, k + 1) get transformed. More specifically, we have that the matrices get updated into We can now restore the MPS form by performing one final SVD of C σ k σ k+1 properly reshaped as a matrix The SVD requires O((r k−1 + r k+1 )r 2 k+1 ) operations and is efficient only if we are able to keep a small bond dimension by truncating the MPS as shown above.
With nearest-neighbor two-site operators at our disposal, we can now present the scheme used to simulate time evolution under Hamiltonians of the form with arbitrary long range interactions. For convenience we first rewrite the Hamiltonian as a sum of, not necessarily nearest-neighbor, pair Hamiltonians h ij as The time evolution operator under this Hamiltonian can be approximated for short times using where the second sequence of products are performed in reverse order. This is the standard second-order Trotter-Suzuki approximation for the time evolution operator [53]. Each one of the operators appearing in Eq. (26) is a SU (4) unitary acting on all pairs of sites for a total of N (N − 1) gates. In order to replace long range gates into nearest-neighbor ones, it is convenient to introduce the SWAP operation. This is defined as When applied to a pair of sites, it exchanges the state as SWAP |φ ⊗ |ψ = |ψ ⊗ |φ .
Using this operation we can transport two sites next to each other, apply the unitary T ij = exp(−iδh ij ) as a nearest-neighbor operation and finally swap the two sites back to their original position [54]. Given the all-to-all nature of the pair interaction in Eq. (25), a naive application of this strategy will result in O(N 3 ) nearest-neighbor operations to implement one step of time evolution using the approximation in Eq. (26) for the propagator.
In this work we employ instead the swap-network scheme introduced in Ref. [51] for studying neutrino dynamics on a digital quantum computer with linear connectivity. The construction is inspired by earlier work on fermionic swap-networks [55] and relies on the observation that, if we perform N layers of nearest-neighbor swap operations by alternating the even bonds with the odd bonds, we can bring every site next to every other site at least once using a total of N (N − 1)/2 operations. At the end of these N steps, the sites will be left in reverse order. It is now clear that, if we define a modified propagator S ij as follows (with h ij from Eq. (25)) and apply nearest-neighbor S ij operations instead of swaps as before, at the end of the sequence we will have implemented the first product of pair operators in the approximation Eq. (26) for the propagator. If we now apply the sequence in reverse order we can complete a full time step and obtain a state where sites have the original order. One can express this more directly by explicitly denoting the sites where the unitary operation is applied as superscript indices S a,b i,j . With this notation, and assuming N even for simplicity, we can express the total unitary Q e k on the k-th even layer and Q o k on the k-th odd layer as follows where i(a, k, e/o) and j(a, k, e/o) denote the amplitude residing on site a at the k-th even or odd step. The mapping between sites and amplitudes evolves as the algorithm proceeds but it is straightforward to keep track of this and assign the correct amplitude index at every point in the step. The full second-order step from Eq. (26) can then be expressed as A schematic depiction of a complete step for a system with N = 4 neutrinos is displayed in Fig. 1. In order to avoid cumbersome notation, in the figure only the lower indices (ij) in the propagators S ij are used and these identify the neutrino amplitude associated with the specific sites the propagator is acting upon at any given step.
In this work we perform a SVD truncation to MPS form after each application of the nearest neighbor operation S ij using a threshold of 10 −10 for the discarded eigenvalues. Better truncation schemes are possible (see eg. [48]) and their potential advantages will be explored in future work.
Finally, we decompose the full time evolution for a total time t into steps of size δ t as We then approximate each short-time propagator for time δ using the swap-network scheme described above.
For all the results presented in the main text, we used a fixed time-step δ = 0.05µ −1 which we found provided a good compromise between precision and efficiency for the systems studied here (see also Appendix A for a more complete discussion). As already mentioned above, the overall scheme remains efficient provided the singular values on each bond decay sufficiently fast to zero that a good approximation to the MPS can be obtained using a low bond dimension. Thanks to the (at most) logarithmic scaling of the entanglement entropy observed in the systems studied in main text, the maximum bond dimension required to keep the truncation error below 10 −10 was ≈ N/2 allowing us to easily explore systems with N 100. The software developed to carry out the simulations reported in this work used the C++ version of the ITensor library which provides efficient representations for MPS and more general tensors [56].

III. RESULTS
The qualitative difference on the collective flavor oscillations in the different dynamical phases can be observed directly from the time evolution of the flavor polarization. Fig. 2 shows the average flavor polarization of the neutrinos in the A beam J A z (t) /(N/4) = 1 − 2p(t), for a system of N = 96 neutrinos and different values for the asymmetry energy δ ω . Due to the conservation of J z , the polarization of beam B is simply the inverse. We can see strong flavor oscillations for 0 < δ ω µ which for δ ω < µ/2 can even lead to a net flavor inversion in the beam (red areas in Fig. 2). In the stable phases for large values of δ ω /µ (bottom panel) or negative asymmetries (top panel), the amplitude of these oscillations is greatly reduced. More detail on these flavor oscillations is given  in panel (a) of Fig. 3 where we show the results for the survival probability p(t) and the same system of N = 96 neutrinos. For negative values of δ ω the polarization experiences small oscillations around the initial state, as expected from our the discussion on the evolution of the total angular momentum based on Eq. (4).
In the phase with 0 < δ ω µ the survival probability experiences instead fast oscillations on time scales much shorter than in the high density limit δ ω = 0. If, in the latter case, the first minimum of the survival probability is reached for a time t P that scales with τ L ∝ √ N , in the former case we have instead evolution at the fast scale τ F ∝ log(N ) found in Refs. [30,37], but now for a physical model with SU (2) invariance (see also [57]). The difference in scaling of t P with system size is evident from the results presented in panel (b) of Fig. 3. For all values of δ ω = 0 a good fit to data is obtained with the ansatz µt P (N ) = a t log(N ) + c t . The best fit values for the parameters extracted from our simulations are reported in Tab. I, together with the extrapolated value of the minimum of the survival probability p min in the N 1 limit. Even though it doesn't directly quantify the final flavor content in the long time limit, this quantity is still remarkably useful as it provides an upperbound on the maximum amplitude of flavor oscillations. In particular we see that, for δ ω outside the critical region (0.0, 1.0), the latter converges to p min ≈ 1 in the thermodynamic limit and no flavor evolution is present. Interestingly, the fluctuation timescales in the unstable region seem to scale approximately as a t ∝ µ/δ ω similarly to the more familiar bipolar oscillations [58][59][60]. The instability observed in our simulations is also similar to bipolar oscillations in that, if we take the beam B to be composed of anti-neutrinos (with negative ω B [31] in normal hierarchy) and beam A of neutrinos, the appearance of the instability leading to oscillations in our setting depends on the neutrino hierarchy: for inverted hierarchy, the single particle energies change sign and the instability cannot occur because δ ω < 0 irrespective of the energies, for normal hierarchy we have always positive energy differences and oscillations can occur when ω A +|ω B | µ. The situation is reversed exchanging neutrinos with anti-neutrinos. Finally, fast oscillations can also occur if both beams are composed by neutrinos (or anti-neutrinos) provided that with + for ν and − forν in the normal hierarchy [61].
In panel (c) of Fig. 3 we show the, normalized, expectation value of the total angular momentum measured using the relation in Eq. (4) derived above. As we can see from the figure, in the unstable region the system can acquire a considerable fraction of the maximum possible total angular momentum, especially for δ ω = µ/2 (blue line) where it reaches a maximum of ≈ 80%.
As an alternative way to quantify quantum correlations in the evolved state, we compute the half-chain entanglement entropy (see eg. [62]) defined explicitly as Here the density matrix ρ B = T r A [ρ(t)] is obtained by tracing the full density matrix ρ(t) =|Ψ(t) Ψ(t)| of the system at time t over the first N/2 neutrinos. The presence of different dynamical phases can be seen directly from the results on the half-chain entanglement entropy S N/2 (t) presented in Fig. 4. Panel (b) shows the maximum value S max that is reached by the entanglement entropy for different values of δ ω and system sizes while in panel (a) we show examples of the time evolution of S N/2 (t) for a system of N = 96 neutrinos. In the phase for δ ω < 0, the behavior is similar to a gapped system with S max independent of system size while we find S N/2 (t) to oscillate in time (red curve) with an approximately constant frequency. We recover a similar result for δ ω µ when, as expected from the discussion above, the system becomes again approximately frozen in the initial configuration. As shown in panel (a) for the case δ ω = µ (orange line), the oscillations in the entanglement entropy are much less regular than in the δ ω < 0 phase and do not return close to zero. At the transition point δ ω = 0 we find that, after an initial rapid growth phase, the entropy S N/2 (t) reaches a peak at S max ≈ log 2 (N/2) and then plateaus at a slightly smaller value with oscillations around the average. This maximum value for the entanglement entropy is reminiscent to the one in ground states of spin systems at a quantum critical point [63,64] and reflects the absence of a gap in the Hamiltonian. This connection allows to understand why it is possible to carry out an efficient simulation even in the infinite interaction range considered here [49]. We find a similar behavior for the maximum entropy reached in the dynamics for the intermediate regime 0 < δ ω /µ 1 as shown in panel(b). For non-zero values of the energy asymmetry however, the half-chain entanglement entropy S N/2 (t) shows strong oscillations analogous to those found in similar spin models undergoing a dynamical phase transition [57]. We note that an increased dynamics of the entanglement entropy has been already associated in the past to the presence of a dynamical phase transition, but was usually associated with a increase of the time-derivative of S N/2 at the critical times [65,66].
The time t ent needed to reach the first maximum in the entanglement entropy is shown in panel (c) of Fig. 4 and has three distinct scaling regimes: for δ ω < 0 is approximately independent from system size, at the critical point we recover the familiar t P ≈ √ N and for positive δ ω > 0 we find a logarithmic evolution t P ≈ log(N ) which seems to approach a constant again for large value of δ ω . In this regime the system seems to enter a second "gapped" phase, similar to the one for δ ω < 0, but simulations with larger systems might be required to confirm this. We note at this point that the slow increase of entanglement entropy with system size (and with time) is precisely what allows us to perform a computationally tractable approximation to the time evolution using Matrix Product States [46,49].

IV. SUMMARY AND PERSPECTIVE
Entanglement is a fundamental feature of interacting quantum many body systems [62,63] which, somewhat surprisingly, can be shown to play no important role in determining ground-state properties of systems with allto-all interactions like for the neutrino forward-scattering Hamiltonian of Eq. (1) (see eg. [67]). This is at the heart of the expectation that a mean-field treatment of the dynamics would also be appropriate. However, the role played by entanglement in out-of-equilibrium settings, like those relevant for flavor evolution in a supernova environment, is not as well understood [33,34].
In this work we have solved for the real time dynamics of a simpler system described by the Hamiltonian in Eq. (2) sharing many similarities to the model used to describe bipolar oscillations [31,59]. In order to overcome the exponential complexity of a direct solution of the problem we use a simple technique, based on a MPS representation [46], whose computational complexity scales with the amount of entanglement generated by the dynamics. The results for this simple model suggest that the maximum bipartite entropy reached in collective neutrino oscillations scales only as log(N ) with the logarithm of the system size. This favorable scaling allows for the efficient calculation of the dynamics for relatively long times and systems exceeding 100 neutrinos. The computational efficiency could also be improved by using more complex tensor networks like MERA [68] which naturally describe states with logarithmically increasing entropy.
As shown in an accompanying paper [57], the behavior observed here is also found in similar (non-integrable) spin systems whose out-of-equilibrium dynamics crosses a quantum critical point, suggesting a strong link between dynamical phase transitions and fast collective oscillations. This link provides us directly with conditions, like Eq. (33), to determine when instabilities can occur. These conditions provide a generalization of those commonly obtained with linear stability analysis to the long time, non-linear, regime.
The artificial setting considered here, with a strongly asymmetric initial state and no effect from the vacuum mixing angle, was chosen to highlight differences with a typical mean-field treatment of the problem where the dynamics is simply absent. We expect however the qualitative conclusions reached here, fast oscillations caused by a dynamical phase transition and a general slow growth of the entanglement entropy with system size, to be valid also in more general situations where multi-angle effects, and more complicated energy distributions are considered [57]. In this regime, MPS simulations could be used reliably to track flavor dynamics up to relatively large, and more realistic, neutrino systems. Important extensions of this work will be to include effects due to finite-mixing angles, the interaction with external matter, the generalization to the 3 flavor case and the more realistic time-dependent setting.
Based on the results presented here, we still cannot completely exclude the existence of neutrino configurations able to generate substantially more entanglement during the dynamics. Interestingly, the presence of these classically-hard configurations would be signalled clearly by the failure of the MPS method to scale efficiently to large system sizes. If found, these would be ideal models to explore using near-term quantum simulators like trapped-ion systems [69] or future circuit based quantum computers [51], which do not necessarily suffer from an exponentially increasing computational cost when large entanglement is present. In future studies, the effects discussed here will be integrated into more realistic simulations to better ascertain the role of entanglement in the dynamical evolution of dense neutrino environments of astrophysical importance. In this appendix we provide more information about the choice of time step in the MPS simulations used in the main text. We first recall the expression for the ap- proximate evolution operator used in our simulations and the very conservative error bound (see e.g. [70] for tighter ones) presented in Eq. (26) of the main text where · are operator norms. This approximation is then useful for short-time evolutions only and in order to perform accurate simulations for a long time t it is necessary to break the interval [0, t] in small steps of size δ t as described in Eq. (32) of the main text. Using the union bound, one can straightforwardly generalize Eq. (A3) to arbitrary long times as follows This results implies that, for a fixed evolution time t, the error scales quadratically with the time step and that to guarantee a maximum error it is sufficient to choose where t/δ is the total number of steps in the simulation. In order to estimate the contribution of this error to the observables displayed in the main text we performed additional simulations for the N = 96 system at both δ ω = 0 and δ ω = µ/2 using a time step δ = 0.2µ −1 four times larger than the one used for the simulations described in the main text. We show the comparison between these two choices of time step for a system with N = 96 neutrinos and two values of the energy asymmetry δ ω for the entanglement entropy in Fig. 5 and for the survival probability in Fig. 6. The top panels show with continuous lines the results for δ = 0.05µ −1 (the choice used for the results in the main text) and the dots indicate the results with δ = 0.2µ −1 which are expected to have ≈ 16 times larger time-step errors. The observed deviations are extremely small for all times and, as shown in the bottom panels of both figures, increase slightly with the total evolution time but are always below ≈ 1% for the entanglement entropy and ≈ 0.2% for the survival probability. Importantly, for times close to the maximum of the entropy the errors are well below the percent level indicating that the time step error is not a significant contribution to the uncertainty of the results presented in the main text.

Appendix B: Direct comparison with mean-field
In this last section we show results obtained from a direct comparison between the MPS representation used in the work presented in the main text and a more standard mean-field solution of the evolution equations. This is done in order to better elucidate the role of the approximations used in the time-evolution algorithm described in the previous section and of possible finite-size effects which are absent in the mean-field.
The model considered in the main text is in the mass basis and in the mean-field approximation there is no flavor evolution. By rotating the initial state toward the flavor basis one can however produce collective oscillations also within the approximation. In this section, we will compare the flavor evolution in the mean-field approximation and the thermodynamical limit N → ∞ with a MPS simulation at finite N where at each time step we truncate the bond dimension to r = 1 on all links.
For this discussion we will use the same Hamiltonian in Eq.(3) of the main text, but rotated to the flavor basis with J A , J B the total spin operators in the two beams and B = (sin(2θ), 0, − cos(2θ)) is related to the mixing angle θ. In the following, we will also keep the same initial condition |Ψ 0 as in the text, but now in the flavor basis. The equation of motion for the spin operators can then be expressed as follows One can eliminate the explicit dependence on the number of amplitudes N by expressing these equations in terms of polarization vectors P A = 4 N J A and P B = 4 N J B . In the mean field approximation we have then The solution of these coupled equations can be simplified with the change of variables proposed in Refs. [59,71]. In Fig. 7 we compare the result for the time evolution of the survival probability as obtained by either solving the equations in Eq. (B3) (black solid line) or using the MPS time-evolution scheme described above with a fixed bond dimension of r = 1 (this figure is equivalent to the solid curve in Fig. 2(b) of Ref. [71], upon rescaling of the time variable). The different colored curves correspond to different choices for the time-step δ used to simulate time-evolution with Eq. (32). The results show a discrepancy with the meanfield results that grows as a function of time and depends on the chosen time-step. This effect is due to the approximate way the evolution is maintained at the mean field level: the local SVD-based truncation of a link tensor to r = 1 does not guarantee to select the globally optimal tensor. This deficiency of the time evolution scheme described in Sec. II B is one of the main motivation behind alternatives like the Time Dependent Variational Principle(TDVP) approach [72] which, instead, approximates globally optimal tensors remaining directly on the manifold of matrix product states with fixed bond dimension. As the operators G σ1σ2 σ 1 σ 2 used in the updates from Eq. (21) become closer to the identity, this truncation error becomes smaller. We can clearly see this effect in Fig. 7 as we lower the size of the time step. These results show the importance of checking convergence of results with respect to the choice of δ (as done here in Appendix A) and also the potential benefits of incorporating ideas like the TDVP scheme in future studies of neutrino oscillations with MPS.