Dynamical Phase Transitions in models of Collective Neutrino Oscillations

Collective neutrino oscillations can potentially play an important role in transporting lepton flavor in astrophysical scenarios where the neutrino density is large, typical examples are the early universe and supernova explosions. It has been argued in the past that simple models of the neutrino Hamiltonian designed to describe forward scattering can support substantial flavor evolution on very short time scales $t\approx\log(N)/(G_F\rho_\nu)$, with $N$ the number of neutrinos, $G_F$ the Fermi constant and $\rho_\nu$ the neutrino density. This finding is in tension with results for similar but exactly solvable models for which $t\approx\sqrt{N}/(G_F\rho_\nu)$ instead. In this work we provide a coherent explanation of this tension in terms of Dynamical Phase Transitions (DPT) and study the possible impact that a DPT could have in more realistic models of neutrino oscillations and their mean-field approximation.

When considering astrophysical settings with large neutrino densities, neutrino-neutrino scattering processes can play an important role in shaping the flavor evolution and can lead to collective oscillations in a neutrino cloud [1,2]. This mechanism has been found to play an important role in extreme environments like the early-universe [3][4][5] or core-collapse supernovae and binary neutron-star mergers [6][7][8][9][10][11][12]. In the latter situations for example, fast neutrino flavor oscillations can lead to important consequences for the revival of the shock wave and nucleosynthesis in the ejected material [13][14][15].
In this work we study simple models of neutrinoneutrino interactions in the forward-scattering limit, when only flavor can be exchanged among neutrinos. For simplicity we also assume that only two flavor of neutrinos mix: ν e corresponding to the electron flavor and ν x , a combination of µ and τ flavors [16]. In this model, neutrinos are mapped into SU (2) flavor isospins evolving at low densities under the vacuum Hamiltonian [17] with σ i = (σ x i , σ y i , σ z i ) the vector of Pauli matrices acting on spin i. The one-body coefficients ω k are connected to the squared mass gap ∆ m = m 2 2 −m 2 1 by ω k = ∆ m /(2E k ), with E k the neutrino energy. The neutrino mass hierarchy is reflected in the sign of the gap: for normal hierarchy we consider ∆ m > 0, while for inverted hierarchy we take ∆ m < 0 [17,18]. The orientation of the "magnetic field" vector B k = (sin(2θ), 0, − cos(2θ)) is related to the mixing angle θ. Importantly, the collective oscillations discussed in this work are not related to the presence of off-diagonal components in the H vac Hamiltonian and in order to avoid confusion we will use a global SU (2) rotation to move to the mass basis |↓ = |ν x and |↑ = |ν e with a diagonal vacuum Hamiltonian.
With the addition of the forward-scattering weak interaction among neutrinos, the full Hamiltonian reads [17] where the interaction strength is given by µ = √ 2G F ρ ν , with G F the Fermi constant and ρ ν the neutrino number density. The geometry of the problem is encoded in the coefficients of the two-body coupling matrix J ij as with p k the momentum associated with the k-th neutrino.
In the low density limit µ ω k , the neutrinos oscillate independently with their own frequency ω k . The presence of the forward-scattering interaction can allow collective effects to develop when µ ω k giving rise to interesting phenomena like synchronization [4,8,19,20], bipolar oscillations [21][22][23] and spectral splits/swaps [24][25][26][27][28]. Due to the computational complexity of solving directly for the dynamics generated by the Hamiltonian H F S for large systems, much of the current understanding of collective oscillation phenomenology is derived within meanfield approaches (see [18] for a review) which, owing to the infinite range of the interaction in Eq. (2), are expected to become increasingly correct as we approach the thermodynamic limit N 1 (this is true in general at least for the ground-state energy, see eg. [29]).
In this work we are interested in understanding the out-of-equilibrium dynamics of the spin model for large but finite systems in order to understand the rate of convergence to the mean field result. Early work by Friedland and Lunardini [30] studied the Hamiltonian in Eq. (2) in the limit where the vacuum term is negligible (high density) and assuming the geometry is isotropic. In this limit, the Hamiltonian is proportional to the total angular momentum operator and therefore easily diagonalizable. The exact solution shows that substantial flavor evolution occurs only for the time scales τ L ≈ µ −1 √ N associated with incoherent scattering. The result is fully consistent with Ref. [31] which argued, using a shorttime approximation, that no entanglement is generated in the many-body evolution of the system and that the mean-field picture of incoherent scattering is correct.
The original study in Ref. [30] was motivated by earlier work by Bell, Rawlinson and Sawyer [32] which presented numerical evidence from a similar model, where however SU (2) invariance was explicitly broken, supporting a very different result: neutrino flavor evolution occurring on much shorter time scales τ S = µ −1 , independently of system size. Despite the infinite range of the pair interaction in H F S , one can expect the time for information to propagate throughout the whole system to be lower bounded by the information signaling time scaling as τ si ∝ log(N ) instead (see eg. [33]). Later work by Sawyer [34] provided additional numerical evidence, with larger system sizes, suggesting indeed the presence of collective flavor oscillations on a fast time scale τ F ≈ µ −1 log(N ).
In the present work, we propose an explanation for the emergence of these different time scales, in apparently very similar models for the neutrino forward scattering problem, as a consequence of the presence of a Dynamic Phase Transition (DPT) [35,36] in the spin system. The models considered in [32,34], and described in more detail in Sec. I below, give rise to fast oscillations with times scaling as τ F by introducing however an unphysical perturbation that breaks the SU (2) invariance of the neutrino Hamiltonian in Eq. (2). As shown recently in a companion paper [37], the presence of of the vacuum Hamiltonian H vac can also produce fast oscillations with times scaling as τ F . In Sec. II we provide additional details about these results and establish a stronger connection with the underlying DPT. Finally, we provide a summary and conclude in Sec. III.

I. HIGH DENSITY LIMIT
It is reasonable to expect that collective effect would be enhanced in the high density limit where µ 1 and the neutrino-neutrino coupling is strong. In the next two sections we will study neutrino systems in the limit where µ |ω k | and neglect the vacuum one-body part from the full Hamiltonian. This contribution will be reintroduced and shown to play an important role in Sec. II below.

A. Single angle approximation
We start our discussion with the model obtained using a very common simplification: the single angle approximation. This amounts to neglect the spatial information encoded in the coupling matrix J ij from Eq. (3) and replace it with it's average value. Here and in the following we will take, without loss of generality, the coupling to be J ij = 1. The final Hamiltonian, after neglecting the one-body vacuum term, can then be written as where we introduced the total flavor spin J = 1 2 i σ i . This model is similar to the Lipkin-Meshov-Glick (LMG) model [38] which, together with it's variants, has been explored extensively in the past [39][40][41][42][43]. The Hamiltonian H sa is diagonal in the angular momentum basis |j, m with j ∈ 0, . . . , N/2 and eigenvalues given by The ground-state is the singlet |0, 0 and the gap to excited states with total spin less than N η/2 vanishes in the thermodynamic limit for any η < 1. Owing to the high degree of symmetry of this model, analytical solutions can be found for the evolution of any observable quantity as a function of time. In particular, a useful observable considered also in Refs. [30,32,37] is the flavor persistence p(t), defined as the probability of measuring one of the neutrinos in the same flavor state it had at the beginning of time evolution. Throughout this work we will consider an initial product state defined as In this case the flavor persistence can be expressed explicitly as the following expectation value where |Ψ(t) = exp(−itH) |Ψ 0 is the time evolved state and, without loss of generality, we have considered the first neutrino which started in the heavy flavor state |↓ at time t = 0. Here and in the following, we will denote the two sets of spins initialized with opposite polarizations in |Ψ 0 as A and B, with corresponding total spin operators J A = (X A , Y A , Z A ) and J B respectively. In order to expose the role of Dynamical Phase Transitions in the collective oscillation phenomenon, we want to describe the full time evolution of the initial state |Ψ 0 under the Hamiltonian in Eq. (4) as a quantum quench [44]. In this setup one starts with an initial Hamiltonian H 0 sa , of which |Ψ 0 is a ground state of, and suddenly changes to the final Hamiltonian H sa given above. With our choice of initial state |Ψ 0 , the initial Hamiltonian we consider in this case can be chosen as where we have indicated with A and B the set of indices for the spins of the A and B group. The Hamiltonian H 0 sa has two degenerate ground-states corresponding to |Ψ 0 and to it's spin-reversed partner obtained by applying the Pauli X operator to each spin: |Ψ 1 = i σ x i |Ψ 0 . The full Hamiltonian used for our quantum quench can then be express compactly as follows with (µ(0), ν(0)) = (0, 1) and (µ(t), ν(t)) = (1, 0) ∀t > 0. The system described by the full Hamiltonian H(t) undergoes a quantum phase transition between a gapped phase for ν(t) µ(t) to a gapless phase for ν(t) µ(t). Contrary to the gapless Hamiltonian H sa , the full Hamiltonian H(t) in Eq. (9) is not diagonal in the coupled angular momentum basis |J, M . Using a mean-field calculation, which is exact in the thermodynamic limit, we find for µ > 0 a critical point at ν = 0 in the thermodynamic limit (see Appendix A 1 for more details). The quench dynamics under consideration here will therefore terminate at the quantum critical point.
In order to define and characterize in general a Dynamical Phase Transition (see [36] for a review) one usually starts by introducing the Loschmidt echo as with |Φ the initial (pure) state at t = 0 and H f the final Hamiltonian of the quench. The quantity L(t) is a fidelity measure [45] that quantifies the probability for the system to return to it's initial state. A DPT is then characterized by non-analiticities in the rate function where N is the total number of particles in the system and λ(t) an intensive "free energy" [35,46]. The rate λ(t) plays here the role of a non-equilibrium equivalent of the thermodynamic free-energy. Notably, other definitions of DPT are possible, for instance using time averaged order parameters [47][48][49] and there are known cases where the two definitions of criticality are incompatible [50]. In the rest of this work we consider only DPT characterized using the Loschmidt echo and leave for future work a more detailed connection to dynamical order parameters. Due to the degeneracy in the ground-space of the initial Hamiltonian H 0 sa , the Loschmidt echo in Eq. (10) needs to be generalized. As shown in Refs. [49,51] a consistent generalization can be found by considering the total probability P (t) of returning to the ground-space where we introduced the two Loschmidt echoes associated with both ground-states. In the thermodynamic limit N 1 only one of the two contribution will dominate resulting in the asymptotic scaling [51] (14) up to exponentially small corrections. The rate functions λ 0 (t) and λ 1 (t) correspond to the definition in Eq. (11) but applied to L 0 (t) and L 1 (t) separately. A DPT can then occur whenever L 0 (t) and L 1 (t) intersect at some finite value t * for the evolution time [49,51]. According to the phase diagram described above, our initial state is quenched up to the critical point and this could lead to a finite value of the crossing time t * for any finite N .
In order to test this scenario, we performed numerical simulations using the Time Evolving Block Decimation (TEBD) algorithm with Matrix Product States (MPS) [52] implemented using the iTensor library [53]. The appealing property of this class of algorithms is that their computational cost scales with the amount of entanglement generated by the real-time dynamics and can then be used efficiently when quantum correlations are sufficiently weak. The implementation of the time evolution operator U (t) = exp(−itH) follows the swap network scheme employed also in past quantum simulations [54]. Additional details on this computational scheme can be found in the companion paper Ref. [37].
The results in the main panel of Fig. 1 show the two Loschmidt echoes L 0 (t) and L 1 (t) for system of different size. In marked difference with the nearest neighbour case studied in Ref. [51], the crossing time t * shows a rapid evolution with system size on time scales proportional to τ L = µ −1 √ N . From the results of our simulations we extract a value of t * /τ L = 1.34(2) for the crossing time. The divergence of τ L with the system size N indicates that this is not technically a DPT, in the sense that the crossing of Loschmidt echoes is a finitesize effect that will vanish in the thermodynamic limit.
The results of our simulation for the the flavor persistence p(t), defined explicitly in Eq. (7), are shown in Fig. 2. We recover the result reported in Ref. [30]: the minimum of the persistence is achieved at times t P ∝ τ L . This is clearly indicated by the inset (b) of Fig. 2 which shows the persistence as a function of the rescaled time t = t/ √ N . The dependence on system size is minimal. The right hand panels show more in detail the system size dependence of t P , in panel (c), and of the value p min (N ) of the persistence at it's minimum, in panel (d). The latter is plotted as a function of 1/N to emphasize the power law scaling of p min (N ) = p min −c/N (the solid green curve in panel (d)). The result of these fit for the minimum time is t P /τ L = 2.10(5) while p min = 0.357(2) in the infinite system size limit. The first two data points in panel (d) of Fig. 2 correspond to N = 8 and N = 10 and we see that one needs to reach N = 16 before deviations from the 1/N behavior are apparent.
All of these time scales quickly diverge for large system sizes N 1 and the mean-field solution, which predicts for this models no time evolution at all, becomes eventually exact in the thermodynamic limit.
In order to quantify quantum correlations in the evolved state, we compute the half-chain entanglement entropy (see eg. [55]) defined as the reduced density matrix obtained by tracing the full density matrix of the neutrino system at time t, denoted as ρ(t), over the first N/2 spins belonging to the A group defined above. We see from the results in Fig. 3 that, after an initial growth, the entropy S N/2 (t) reaches a peak and then plateaus at a value S max ≈ log 2 (N/2) with oscillations around the average. The maximum value S max for the entanglement entropy is reminiscent to the one in ground states of one dimensional spin systems at a quantum critical point [56,57] and reflects the absence of a gap in the Hamiltonian H sa in Eq. (4). The qualitative behavior of S N/2 (t) is remarkably close to the one observed with a similar model (but different initial conditions) in Ref. [58]  where the entanglement entropy was observed to peak and then plateau when the system was quenched at the critical point of a DPT. The observed time scale to reach the peak, also connected to the Eherenfest time t Ehr [58], was found there to scale as t Ehr ≈ log(N ) similarly to the fast scale τ F while away from the quantum critical point t Ehr ≈ √ N like τ L . From our simulation we find that in our case, despite being at the critical point, the entropy grows more slowly and reaches the peak on the slow time scale t ent ≈ √ N . In order to account for finite size effects, we perform a fit to the data shown in the inset of Fig. 3 using The optimal parameter for the leading order term is found to be aτ L = 1.16(4) while the finite size corrections b and c are O(10). This time scale is very similar, and always strictly smaller, to the crossing time t * when the DQPT occurs (see red dashed line in inset of Fig. 3). A separate test of whether the entanglement time t ent scales algebraically (case α) or logarithmically (case β) in system size can be obtained by estimating the time to reach S N/2 = S max using two limiting cases For the single angle setup considered in this section we found a good fit to data only for the model from "case α" (shown as purple dotted line in Fig. 3) with optimal parameters A = 2.14(4) and B = 2 respectively.
We note that this slow increase of the entanglement entropy with system size and with time is at the hearth of the classical simulatability of the neutrino model in the single-angle approximation with Matrix Product States: the maximum bond dimension needed only scales linearly with N to obtain converged results. The TEBD scheme employed here, and in the accompanying paper [37], is however not optimal for long range interactions and further progress could be made using more sophisticated simulation techniques like the Time Dependent Variational Principle [59] as well as different tensor network [60,61] or neural network states [62].
The first calculations showing a many-body "coherent speedup" of flavor oscillations at the shorter time-scale τ S = µ −1 were obtained in Refs. [32,34] using a neutrino Hamiltonian that explicitly breaks the global SU (2) flavor invariance of the Hamiltonian H F S in Eq. (2). The symmetry-breaking term used in both cases is The control parameter here is ∆ and for ∆ = 1 the original SU (2) invariant interaction is recovered. As we will see below the geometry of the problem encoded in the angular factors J ij will play now an important role. In this section we will consider a very simple situation: two neutrino beams, one with N/2 neutrinos starting in the |↓ state and one with N/2 neutrinos starting in |↑ . These correspond to the sets A and B defined above. Neutrinos belonging to the same group interact with the same strength J AA = J BB while neutrinos belonging to different beams interact with a coupling J AB . Using the total flavor spin operators J A and J B introduced above, we can write the full Hamiltonian used in this quench as plus an inconsequential constant factor that we ignore. The limit in which the beams are very collimated corresponds to the choice J AA = 0, and the weak interactions are relevant only across beams. For the rest of this section we will measure energies in units of (µJ AB ) and use directly the dimensionless parameter Γ = J AA /J AB . The quench dynamics we will consider in this section starts in the limit ∆ → ∞ with Γ → 0 which corresponds to the starting Hamiltonian H 0 sa considered above, with |Ψ 0 and |Ψ 1 as it's two degenerate ground-states.
The equilibrium phase diagram of the two beam Hamiltonian H tb is now much richer than with the single angle approximation (see Fig. 4). Using a mean-field approach (see App. A 2 for details) we can identify 4 distinct phases depending on the value of the SU (2) breaking parameter ∆ and on the ratio Γ of the two body couplings which specifies the relative orientations of the beams.
For collimated beams with Γ < 1 we find two gapped phases, one with anti-ferromagnetic order in the z direction at large positive values of ∆ (denoted by AF M ) and one with ferromagnetic order along the z direction for sufficiently negative values of ∆ (denoted by F M ). These two phase are separated by a gapless phase, indicated by XY in Fig. 4, where anti-ferromagnetic order is preserved in the xy−plane but is lost in the z direction. The only ordered phase present for Γ ≥ 1 is the F M phase for ∆ < 0 while a disordered gapless phase emerges for positive values of ∆ (denoted by DIS in Fig. 4).
The results within the single angle approximation described in the previous section (and in Ref. [30]) correspond to the trajectory indicated by the F L arrow in Fig. 4 and ending at the full dot (which indicates the single angle point). The dashed lines emanating from that point indicate parameter values for which, due to conservation laws, the dynamics is indistinguishable from the one obtained with the F L quench. Note that this holds also for quenches that are apparently crossing a phase boundary. This is in agreement with previous studies showing that a DPT can fail to appear even in quenches that crossed a phase boundary (see eg. [49,63,64]).
Here we study in some detail the quench used in the original paper by Bell et al. in Ref. [32] (denoted by the BRS arrow in Fig. 4) and comment on the qualitative differences with the single angle case explored in the previous section. Further exploration of the interplay between the equilibrium phase boundaries displayed in Fig. 9 and the presence of a DPT would be very interesting. However, since in order to describe neutrino interactions we are not allowed to break the SU (2) invariance explicitly, we cover here only the simpler case needed to explain the findings of Refs. [32,34] and proceed in the next section to consider instead the SU (2)-invariant problem considered already in Ref. [37] which shows similar features.
We start by looking at both the time evolution of the flavor persistence p(t) and the crossing time of the two Loschmidt echoes from Eq. (13). The main panel of Fig. 5 shows the flavor persistence p(t) for various system sizes (solid lines) together with the equivalent result in the single angle approximation from the previous section (dotted lines). It is clear that flavor evolution happens much faster in the BRS quench, with sustained oscillations for long times. The frequency of these oscillations, as measured by the time t P to reach the first minimum, follows the fast time scale τ F = µ −1 log(N ) as µt P = 2.04(5) log(N ) + 1.6(1) (see panel (b) of Fig. 5). This is in agreement with the expectations from results presented in Refs. [32,34] and much faster than in the single angle approximation (shown as the green dashed line in Fig. 5(b)) we studied above and in Ref. [30].
The Loschmidt echoes L 0 (t) and L 1 (t) are also found to cross at shorter time scales than those found in Sec. I A.
The results for the crossing time t * as a function of the system size N are presented in panel (c) of Fig. 5 and again follow the fast time scale with t * /τ F = 1.56 (4).
The stark difference with the single angle case can also be observed in the evolution of the half-chain entanglement entropy S N/2 defined in Eq. (15). The main panel  of Fig. 6 shows the entanglement entropy for different system sizes N = 8, 16, 24, 32, 48, 64, 96, 128 (solid lines in the main panel). The behavior in this case is qualitatively different from the results shown in Fig. 3 for the single angle approximation: the entanglement entropy itself oscillates in time, reaching values as high as S max (dashed lines in Fig. 6) multiple times. In the results shown in Fig. 6 we see two distinct peaks whose times scale with the fast time scale τ F as t 1 ent /τ F = 1.3(1) and t 2 ent /τ F = 3.9(1) respectively (these fits are shown in the inset of Fig. 6 as continuous lines). To corroborate these findings we also shown in the main panel is the "case β" fit from Eq. (17) which very accurately matches the evolution of the entropy maximum.
The results shown in this section were obtained using Γ = 0 as in the original model from Ref. [32] but we confirmed the presence of the same logarithmic time scale also for larger values up to Γ ≈ 0.7 as observed also in previous work as reported in Ref. [34]. The original model from Ref. [32] also used a more complex angular distribution than the two beam geometry employed here and in Ref. [34], unfortunately, due to the explicit N dependence of the angular distribution used there, it was not possible to obtain a smooth extrapolation in system size as we have done with the other models in this work. We have found in a few selected cases at fixed N that, with our initial state |Ψ 0 , more complex angular dependence actually slows down the dynamics as compared to the two beam geometry. This effect is likely due to frustration of some of the interaction terms and in future work we plan to assess more quantitatively the role of multi-angle effects by using model geometries that have a well-defined scaling with system size.

II. INTERMEDIATE DENSITY REGIME
The fast flavor oscillations observed in the models of the previous section are unfortunately not directly relevant to neutrino physics since the correct Hamiltonian is SU (2) flavor invariant also in the general case. The previous result, however, points to the fact that oscillations at the time scale τ F can appear when one crosses a quantum critical point and we have a DPT in the quantum quench. By tuning appropriately the one body part of the forward-scattering Hamiltonian in Eq. (2) we can orchestrate this to happen also in a physically relevant scenario closely related to the model used in describing bipolar collective oscillations (see eg. [22,65]).
In this section we will consider the same model we introduced in the companion paper [37] where the system is still decomposed in the two beams A and B but now with two different energies where we have also used the single angle approximation for the coupling matrix J ij in the interaction. The Hamiltonian commutes with the z component of the total flavor spin J z = Z A +Z B and, given our initial state |Ψ 0 , it's expectation vale remains zero at all times. Using spin operators for the neutrinos in the two beams and denoting the spin difference along the z axis as D z = Z B − Z A , we can write the full Hamiltonian as (cf. [37]) where we introduced δ ω = (ω A − ω B )/2 for the energy difference between the two beams and dropped an irrelevant constant. The equilibrium phase diagram depends on the sign of the two body interaction µ: • for a ferromagnetic coupling µ < 0, there is a second order transition at δ ω = ±|µ| between polarized phases with D z = ∓N/2 and a broken phase with ferromagnetic order in the xy plane [41].
• for an anti-ferromagnetic coupling µ > 0, the transition between gapped polarized phases is of first order and at δ ω = 0 instead [39].
On the other hand, the Loschmidt echo Eq. (10) characterizing a DPT is invariant upon inversion of the full Hamiltonian H ID → −H ID and we can therefore expect the dynamical phase diagram to display features of both cases above and depend instead only on the relative sign of the two couplings constants µ and δ ω . This is indeed the case as shown in the results presented in Ref. [37] which we briefly summarize here. Using energy conservation together with the known initial state |Ψ 0 whose energy expactatin value reads we can express the instantaneous value of the total angular momentum as a function of the staggered spin polarization D z as follows with initial conditions J 2 (0) = D z (0) = N/2. As was show in the accompanying paper [37], this relation between the total angular momentum and the flavor asymmetry in the two beams is sufficient to characterize qualitatively the entire out-of-equilibrium dynamics. For completeness we provide a more complete derivation of those results with more details in the following.
In the case where the energy asymmetry δ ω /µ < 0 is negative, the total spin, which starts already at a relatively small value, can only decrease further during time evolution. Since the operator J 2 is positive semi-definite this introduces a constraint on the fluctuations that D z can experience, in particular and the change in polarization per spin vanishes in the thermodynamic limit. This suggests that for δ ω /µ < 0 the system experiences negligible flavor evolution and is always close to the initial state, this was called the frozen phase in Ref. [37]. In the opposite limit δ ω /µ > 0 instead, the fluctuations become parametrically small at low densities (corresponding to δ ω /µ 1 ) but remain finite also in the N 1 limit This inequality provides a nontrivial bound on the spin, or flavor, fluctuations only for large δ ω > µ/4. Based on the discussion of the equilibrium phase diagram of this model, we expect to find the system in the gapped polarized phase, with D z large and positive, for sufficiently large δ ω /µ values. An estimate for the transition can be obtained by considering the minimum value of δ ω /µ for which the first order fluctuations preserve the sign of the order parameter. This can be obtained by ensuring with Using the fact that J z = J 2 z = 0 for our initial state, we can find the following upperbound on the variance We therefore expect the system to be in the polarized phase and experience little flavor evolution when and possibly at somewhat smaller values due to the bound Eq. (27) being not tight. Finally, in the regime 0 < δ ω /µ ≤ 1/4 the total spin J 2 , and correspondingly the flavor difference D z , can experience strong fluctuations bounded by As expected from this qualitative discussion, the dynamical phase diagram delineated above corresponds to a combination of the equilibrium phase diagrams of both the ferromagnetic and anti-ferromagnetic cases, with the exception that the transition at large δ ω /µ appears shifted to larger values than δ ω /µ = 1.
As shown also in Ref. [37], the presence of these different dynamical phases is directly visible in the time evolution of the half-chain entanglement entropy for different values of the one body energy asymmetry δ ω /µ. In the frozen phases for either δ ω /µ < 0 or δ ω /µ 1 the entanglement entropy remains small with a maximum value independent of system size. For negative energy asymmetry δ ω /µ the entropy experiences fast oscillations which bring S N/2 back to zero periodically. This is shown in the top panel of Fig. 7 showing the evolution of the half-chain entropy for a system of N = 96 neutrino amplitudes across the different dynamical phases. In the gapless region 0 < δ ω /µ 1 the entanglement entropy shows strong fluctuations as a function of time, with maximum values close to S max = log 2 (N/2) and monotonically decreasing for increasing one-body energy asymmetry (see also Fig.3 of Ref. [37]). The special case δ ω /µ = 0 matches the behavior presented in Fig. 3   late times. The scaling of time scales in the half-chain entropy for the unstable region 0 < δ ω /µ 1 shows a logarithmic behavior as expected from the presence of a DPT into a gapless phase, similarly to what we have found for the SU (2)-broken model in Sec.I B.
In order to establish a closer connection to dynamical phase transitions as defined in the previous sections, we now consider the evolution of the Loschmidt echo L(t) for different values of the asymmetry parameter δ ω in all three dynamical phases. Contrary to the situation in Sec. I A and Sec. I B, the initial Hamiltonian we consider here (namely H ID with δ ω < 0 and µ = 0) has a unique ground-state. For all quenches with δ ω = 0 considered in this section, we have always found L 1 (t) ≈ 0 for large system sizes and a DPT will not appear as a crossing of echoes as before, but instead as sharp peaks in the Loschmidt rate λ(t) defined in Eq. (11) above. This is illustrated in of Fig. 8 where the Loschmidt rate λ(t) is shown for different values of δ ω in a system of N = 96 neutrino amplitudes. The purely two-body case at δ ω = 0 has a DPT generated by crossing Loschmidt echoes at t = t * (shown as a dot in Fig. 8), followed by additional sharp features at later times. For negative values of δ ω , in the frozen phase, the rate λ(t) remains smooth at all times, while for positive δ ω sharp features start to appear at even shorter times than the t * crossing time and a DPT can occur in the system. Obtaining an estimate for the critical time where a DPT might occur in this case is complicated by it's expected evolution with system size, in parallel to the case δ ω = 0 considered in Sec. I A above. This has prevented a reliable extraction of a unique critical time t * in the unstable region 0 < δ ω /µ 1 using results up to N = 128 and a single value for the time-step of the evolution (here we used 0.05µ −1 as in Ref. [37]). This observation highlights the usefulness of entanglement measures such as the half-chain entropy as a more robust indicator of the presence of qualitative changes in the dynamical phase of a many-body quantum system. Future explorations employing either semi-classical approaches, like those used for instance in [50], or specialized simulations exploiting more directly symmetries of the system, are expected to be able to clarify the role of fidelity measures as the Loschmidt echo in characterizing the different dynamical phases found in models of neutrino flavor evolution.

III. SUMMARY AND CONCLUSIONS
The presence of collective oscillations in the dynamical evolution leading to neutrino flavor transport has long been recognized as an important effect in describing the dynamics of astrophysical environments like supernovae and the early universe [1,2,4,6]. Early explorations by Sawyer and coworkers [32,34,66] suggested that quantum correlations, in the many-body spin system corresponding to a neutrino cloud, could lead to a coherent speed-up of collective oscillations, with possibly important consequences for the dynamics of these environments. This idea, which invites caution on the interpretation of results for the neutrino flavor evolution obtained using mean-field approximations (which neglect quantum entanglement), has been challenged in the past by presenting counter-examples in solvable models where the qualitative prediction of the mean-field are matched by the exact solution [30,67]. The absence of entanglement in the neutrino dynamics more generally has also being argued as a justification for the mean-field approach to the problem [31]. This debate has recently re-emerged thanks to works like Ref. [68] and Ref. [69] which showed that entanglement is indeed produced when solving exactly the many-body neutrino problem encoded in the forward scattering Hamiltonian of Eq. (2) and it's timedependent generalizations. The explored systems were however too small (N = O(10)) to draw general conclusions applicable to the large collections of neutrino amplitudes needed for realistic simulations.
Exploiting the expectation that the entanglement entropy is unlikely to grow too large in these many-body systems, due to the infinite range of interactions in the spin model of Eq. (2), the present work extends the idea presented in the companion paper Ref. [37] to use a Matrix Product State (MPS) representation in order to efficiently describe the neutrino wave-function as it evolves from an initial product state. As explained in more detail in Ref. [37], this approach is ideal for low levels of bipartite entanglement in the system and allows to easily simulate systems with ≈ 100 neutrino amplitudes with modest computational resources. This simulation strategy is used here with two main goals, the first one was to validate the early small scale simulations by Sawyer et al. [32,34] which, correctly, predicted flavor evolution to occur (in their model) at a fast time-scale τ F ≈ µ −1 log(N ). This shows that indeed many-particle neutrino interactions cause a novel coherent effect not captured by the mean-field approximation. A similar effect is also found in the more familiar bipolar oscillations described in detail in Ref. [37] and Sec. II of the present work. The second goal was to explain the presence of this fast time scale as being generated by an underlying Dynamical Phase Transition. This observation explains the absence of the effect in the exactly solvable models discussed in Refs. [30,67] and provides a more direct link between the presence of coherently-enhanced flavor oscillations and non-negligible levels of entanglement in the many-body state generated by the dynamics.
The work presented here and in the accompanying paper Ref. [37] opens the way to accurate many-body simulation of the full quantum dynamics of neutrino flavor transport with controllable errors. The use of entanglement-efficient methods, like the MPS representation used here, will allow for the first time a more direct comparison with popular approximation methods working in the mean field for large system sizes. This will be critical to allow for the inclusion of rich energy/angle distributions and avoid the limitations of special symmetric points like the model studied in Ref. [30] and covered in Sec. I A of the present work. Possible failures of this program would be associated to situations where the entanglement entropy grows substantially with system size. The identification of the parameter regimes where this happens would shed light on potentially interesting candidates to study using quantum computing devices as recently explored in Ref. [54]. Finally, a better understanding of the dynamical phase diagram of neutrino models, as the one described in Eq. (2) and it's generalization to the full 3 flavor case, would help identify the conditions (beyond linear stability analysis) required for collective oscillations to appear in complex environments like supernovae explosions by an appropriate analysis of simulation results. Work is ongoing to extend the results presented in this work to more realistic conditions in order to better asses the impact of entanglement in astrophysical settings with large neutrino densities.
with µ, ν positive constants. As discussed in the main text, in the limit ν = 0 the system is gapless and the groundstate has zero total angular momentum and zero energy. In the limit µ = 0 instead, the system has two degenerate ground-states which, in the angular momentum basis |s A , m A ⊗|s B , m B of the two set of spins with total angular momenta S A and S B , we can write as In these configurations the system has an antiferromagnetic order across beams characterized by Z A Z B = −N 2 /16. In the gapless phase the order parameter is zero. The expectation value of the full Hamiltonian in either of the anti-ferromagnetic states reads and becomes negative for a sufficiently large antiferromagnetic coupling ν > 8µ/N . In the thermodynamic limit we expect the critical point to be at ν = 0 for any µ > 0. As we will see in a more general case below, if we allow ν to become negative other phases emerge.

Phase diagram of the two-beam model
In this section we provide more details on the calculation of the mean-field phase diagram presented in Fig. 4 of the main text. This corresponds to the ground-state phase diagram of the following Hamiltonian (cf. Eq. (19)) with a positive coupling constant Γ = J AA /J AB . The order parameters of interest here are the average staggered magnetizations of the two beams We start the discussion of the equilibrium phase diagram by considering first some special cases: • at the SU (2) symmetric point, corresponding to ∆ = 1, we have the following Hamiltonian For Γ < 1 the system is in an anti-ferromagnetic gapless phase characterized by M V AB = − N 16 and undefined values for M XY AB and M Z AB . For Γ > 1 we have instead a disordered gapless phase characterized by a vanishing order parameters M V AB = M XY AB = M Z AB = 0. At the single angle point Γ = 1, the three order parameters are undefined. Note that, when the initial state is |Ψ 0 from Eq. (6), the resulting evolution is the same for any value of Γ since S 2 A and S 2 B are conserved quantities.
• at the single angle point Γ = 1 we have instead with Z tot = Z A +Z B the total spin in the z direction (and similarly for X tot and Y tot ). The groundstate of this model for ∆ < 0 is a (gapped) ferromagnet with M Z AB = N 16 , for ∆ ≥ 0 the groundstates are the singlet states with zero total spin and with undefined order parameters. Given our initial state |Ψ 0 , and the fact that [Z tot , H Γ=1 ] = 0, the time evolution is exactly equivalent to the the single angle case studied above for any value of ∆.
• for collimated beams with Γ = 0 we have simply For ∆ > 1 the ground states are |Ψ 0 and the spin reversed partner |Ψ 1 introduced in Sec. I A of the main text. The system has anti-ferromagnetic order with M Z AB = − N 16 and there is a finite energy gap to excited states. For −1 < ∆ < 1 the system is gapless with M Z AB = 0, in fact we have a continuum of zero-energy modes polarized in the XY plane with M XY AB = − N 16 . Finally, for ∆ < −1 the system is a ferromagnet along the Z direction with M Z AB = N 16 and M XY AB = 0.
In order to get the rest of the phase diagram we will compare energies of the different phases in the mean field limit. Let's first rewrite the Hamiltonian as The mean field states we will consider here are: together with the disordered state |Φ DIS with zero total spin in beam A and B. In the expression above we use the notation |± to indicate the eigenstates of the Pauli X operator with positive and negative eigenvalue respectively. The corresponding expectation values for the energy in the full Hamiltonian Eq. (A4) are The resulting phase diagram is depicted in Fig. 9. Along the critical lines separating the different phases we have the following • boundaries between AF M and DIS and between XY and DIS (dashed black curve in Fig. 9): all the order parameters are undefined due to the degeneracy of the spectrum for states with different values of the total spin in the two beams but zero total angular momentum.