Entanglement and Many-Body effects in Collective Neutrino Oscillations
IIQuS@UW-21-02
Entanglement and Many-Body effects in Collective Neutrino Oscillations
Alessandro Roggero InQubator for Quantum Simulation (IQuS), Department of Physics,University of Washington, Seattle, WA 98195, USA (Dated: February 23, 2021)Collective neutrino oscillations play a crucial role in transporting lepton flavor in astrophysicalsettings, such as supernovae, where the neutrino density is large. In this regime, neutrino-neutrinointeractions are important and simulations in mean-field approximations show evidence for collectiveoscillations occurring at time scales much larger than those associated with vacuum oscillations. Inthis work, we study the out-of-equilibrium dynamics of a corresponding spin model using MatrixProduct States and show how collective bipolar oscillations can be triggered by quantum fluctuationsif appropriate initial conditions are present. The origin of these flavor oscillations, absent in themean-field, can be traced to the presence of a dynamical phase transition, which drastically modifiesthe real-time evolution of the entanglement entropy. We find entanglement entropies scaling at mostlogarithmically in the system size suggesting that classical tensor network methods could be efficientin describing collective neutrino dynamics more generally.
Neutrinos play a pivotal role in extreme astrophysi-cal events like core-collapse supernovae and neutron starmergers, where they are responsible for both reinvigo-rating a stalled shock-wave and controlling the condi-tions for nucleosynthesis in the ejected material [1, 2].In these environments with a large neutrino density ρ ν ,neutrino flavor evolution is substantially modified byneutrino-neutrino scattering processes which can lead toself-sustained collective flavor oscillations [3–12]. Sinceneutrinos in supernovae are emitted with fluxes and spec-tra that are strongly flavor dependent [2], the presenceof collective flavor oscillations could then lead to impor-tant effects [13–15]. For large neutrino densities nearthe proto-neutronstar, flavor oscillations can occur at fre-quencies of the order of the neutrino-neutrino interactionenergy µ = √ G F ρ ν instead of the much smaller vacuumoscillation frequency ω = ∆ m / (2 E ), a phenomenonknown as fast flavor conversion [16–19] (see [20, 21] forreviews). In these expressions G F is Fermi’s constant, E the neutrino energy and ∆ m = m − m the mass gapin the two flavor approximation used in this work.The standard approach to study these collective phe-nomena is to neglect incoherent scattering and onlyconsider neutrino-neutrino interactions in the forward-scattering limit where only flavor can be exchangedamong neutrinos. In this limit, a finite system with N neutrinos is represented as a collection of SU (2) flavorisospins governed by the Hamiltonian [22] [23] H fs = N (cid:88) k =1 ∆ m E k (cid:126)B · (cid:126)σ k + µ N N (cid:88) i We consider a simple model, in the mass basis, in whichthe neutrinos are divided into two equally sized ”beams”( A and B ) and initialized as the following product state | Ψ (cid:105) = N/ (cid:79) n =1 |↓(cid:105) ⊗ N/ (cid:79) m =1 |↑(cid:105) . (2)This corresponds to choosing all the neutrinos in the A beam to be | ν (cid:105) and all the neutrinos in the B beam tobe | ν (cid:105) . We will consider the simple case where the ener-gies of the neutrinos in beam A are all equal to E A and a r X i v : . [ h e p - ph ] F e b similarly for the neutrinos in beam B . Finally, we will ne-glect the geometric information encoded in the couplingmatrix J ij and take the single angle approximation bysetting J ij = 1 for all neutrinos [20, 22] (see also [33]).The resulting neutrino Hamiltonian, similar to the oneconsidered in Ref. [25], can be written as follows H = − ω A (cid:88) i ∈A σ zi − ω B (cid:88) i ∈B σ zi + µ N (cid:88) i The qualitative difference on the collective flavor oscil-lations in the different dynamical phases can be observeddirectly from the time evolution of the flavor polariza-tion. In Fig. 1 we show the average flavor polarization ofthe neutrinos in the A beam (cid:104) J Az ( t ) (cid:105) / ( N/ 4) = 1 − p ( t ),for a system of N = 96 neutrinos and different values forthe asymmetry energy δ ω . Due to the conservation of J z , = / = = / = / = / = Time after quench t [ ] F l a v o r p o l a r i z a t i o n o f b e a m A FIG. 1. (Color online) Flavor polarization per particle (cid:104) J Az ( t ) (cid:105) / ( N/ 4) of neutrinos in the A beam as a function oftime for six values of the energy asymmetry parameter δ ω /µ (from top to bottom): − . , . , . , . , . , . Time after quench t [ µ −1 ] F l a vo r p e r s i s t e n ce p ( t ) δ ω = - µ/2δ ω = 0δ ω = µ/4δ ω = µ/2δ ω = µ Time after quench t [ µ −1 ] < J ( t ) > / J m a x System size N M i n i m u m ti m e t P [ µ − ] t~ √ N (a) (b)(c) t~log(N)t ~ const FIG. 2. (Color online) Panel (a) shows the flavor persis-tence p ( t ) in a system of N = 96 neutrinos initialized in | Ψ (cid:105) and quenched at different values of the one body asymme-try parameter: the red line at the top is for a negative value δ ω = − µ/ 2, the black line contains only the two body poten-tial and the green, blue and orange lines correspond to pos-itive asymmetries δ ω = (0 . , . , . µ respectively. Panel(b) shows the evolution of the minimum time t P as a func-tion of system size and panel (c) shows the (normalized) timedependent expectation value of the total angular momentumfor different values of δ ω in a system of N = 96 neutrinos. the polarization of beam B is simply the inverse. We cansee strong flavor oscillations for 0 < δ ω (cid:46) µ which for δ ω < µ/ δ ω /µ (bottom panel) or negative asymmetries(top panel), the amplitude of these oscillations is greatlyreduced. More detail on these flavor fluctuations is givenin panel (a) of Fig. 2 where we show the results for thepersistence p ( t ) and the same system of N = 96 neutri-nos. We see that for negative values of δ ω the polarizationexperiences small oscillations around the initial state, asexpected from our the discussion on the evolution of thetotal angular momentum based on Eq. (5).In the broken phase for 0 < δ ω (cid:46) µ the flavor persis-tence experiences instead fast oscillations on time scalesmuch shorter than in the high density limit δ ω = 0.If, in the latter case, the minimum of the persistence isreached for time scales proportional to τ L ∝ √ N , in theformer case we have instead evolution at the fast scale τ F ∝ log( N ) found in Refs. [24, 29], but now for a physi-cal model with SU (2) invariance (see also [45]). The log-arithmic dependence on system size of the time to reachthe minimum t P is evident from the results presented inpanel (b) of Fig. 2. For all values of δ ω (cid:54) = 0 a good fit todata is obtained with the ansatz µt P ( N ) = a t log( N )+ c t .The best fit values for the parameters extracted from oursimulations are reported in Tab. I, together with the ex-trapolated value of the minimum of the persistence p min δ ω /µ a t b t c t p min -0.5 0 0 1.8(2) 1.00(1)0.0 0 2.10(5) 0 0.36(1)0.01 5.4(2) 0 -8(1) 0.35(2)0.05 2.58(6) 0 0 0.31(4)0.125 1.65(6) 0 1.4(2) 0.32(3)0.25 1.22(4) 0 1.6(2) 0.41(3)0.5 1.06(4) 0 1.4(2) 0.61(4)1.0 0.96(5) 0 0 0.99(2)TABLE I. Optimal parameters of the model discussed in themain text µt P ( N ) = a t log( N ) + b t √ N + c t for time to reachthe minimum of the persistence t P . The minimum of thepersistence in the thermodynamic limit is extracted using theansatz p min ( N ) = p min + p c /N γ with exponents γ = 1 for δ ω = 0 and γ = 0 . in the N (cid:29) δ ω out-side the critical region (0 . , . p min ≈ a t ∝ (cid:112) µ/δ ω similarly to the more familiar bipolar oscil-lations [46–48]. The instability observed in our simula-tions is also similar to bipolar oscillations in that, if wetake the beam B to be composed of anti-neutrinos (withnegative ω B [20] in normal hierarchy) and beam A ofneutrinos, the appearance of the dynamical phase tran-sition leading to oscillations in our setting depends onthe neutrino hierarchy: for inverted hierarchy, the singleparticle energies change sign (accounting for ∆ m < δ ω < ω A + | ω B | (cid:46) µ . The situationis reversed exchanging neutrinos with anti-neutrinos. Fi-nally, fast oscillations can also occur if both beams arecomposed by neutrinos (or anti-neutrinos) provided that0 ≤ ± E B − E A E A E B (cid:46) √ G F ρ ν ∆ m , (6)with + for ν and − for ¯ ν in the normal hierarchy [49].In panel (c) of Fig. 2 we show the, normalized, expec-tation value of the total angular momentum measuredusing the relation in Eq. (5) derived above. We use J max = N ( N + 2) / δ ω = µ/ ≈ (cid:104) (cid:126)J ( t ) (cid:105) = (cid:104) J z ( t ) (cid:105) = 0 validat any point in the time evolution, we see that any value -2 -1 0 1 2 One body energy difference δ ω [ µ] T i m e t o m a x i m u m t e n t [ µ − ] Time after quench t [ µ −1 ] E n t a ng l e m e n t e n t r opy S N / δ ω = µδ ω = µ/2δ ω = 0 δ ω = - µ/2 M a x i m u m e n t r opy S m a x N = 16N = 32N = 64N = 128 N = 96 t~ √ N t ~ l og ( N ) t ~ c o n s t t ~ c on s t (a) (b)(c) FIG. 3. (Color online) Half chain entanglement entropy S N/ for N = 96 neutrinos and different values of the one-bodyenergy difference: the orange line is δ ω = µ , the blue line isfor δ ω = µ/ 2, the red line for δ ω = − µ/ δ ω = 0. Panel (b) shows the maximum entropy S max forthree system sizes N = 16 , , , 128 and different values of δ ω . Panel (c) shows the time t ent to reach the maximum. of (cid:104) J (cid:105) that deviates from the initial value of N/ J , explored in the dynamics with δ ω /µ < 0, isobtained for singlet states with zero total spin and vari-ance [52, 53]. Higher values of total angular momentum,as obtained in the unstable region, are related instead tosymmetric Dicke states with m excitations [52, 54].As an alternative way to quantify quantum correlationsin the evolved state, we compute the half-chain entangle-ment entropy (see eg. [55]) defined explicitly as S N/ ( t ) = − T r [ ρ B ( t ) log ( ρ B ( t ))] . (7)Here the density matrix ρ B = T r A [ ρ ( t )] is obtained bytracing the full density matrix ρ ( t ) = | Ψ( t ) (cid:105)(cid:104) Ψ( t ) | of thesystem at time t over the first N/ S N/ ( t ) presented in Fig. 3. Panel (b) shows themaximum value S max that is reached by the entangle-ment entropy for different values of δ ω and system sizeswhile in panel (a) we show examples of the time evolu-tion of S N/ ( t ) for a system of N = 96 neutrinos. Inthe frozen phase, for δ ω < 0, the behavior is similar toa gapped system with S max independent of system sizewhile we find S N/ ( t ) to oscillate in time (red curve) withan approximately constant frequency. We recover a sim-ilar result for δ ω (cid:38) µ when, as expected from the dis-cussion above, the system becomes again approximatelyfrozen in the initial configuration. As shown in panel (a)for the case δ ω = µ (orange line), the oscillations in theentanglement entropy are much less regular than in thefrozen phase and do not return close to zero. At the tran-sition point δ ω = 0 we find that, after an initial rapid growth phase, the entropy S N/ ( t ) reaches a peak at S max ≈ log ( N/ 2) and then plateaus at a slightly smallervalue with oscillations around the average. This maxi-mum value for the entanglement entropy is reminiscentto the one in ground states of spin systems at a quantumcritical point [56, 57] and reflects the absence of a gapin the Hamiltonian. We find a similar behavior for themaximum entropy reached in the dynamics for the inter-mediate regime 0 < δ ω /µ (cid:46) S N/ ( t ) shows strong os-cillations analogous to those found in similar spin modelsundergoing a dynamical phase transition [45]. We notethat an increased dynamics of the entanglement entropyhas been already associated in the past to the presenceof a dynamical phase transition, but was usually associ-ated with a increase of the time-derivative of S N/ at thecritical times [58, 59].The time t ent needed to reach the first maximum inthe entanglement entropy is shown in panel (c) of Fig. 3and has three distinct scaling regimes: for δ ω < t P ≈ √ N and for positive δ ω > t P ≈ log( N ) whichseems to approach a constant again for large value of δ ω .The system seems to enter a second gapped phase in thisregime but simulations with larger systems might be re-quired to confirm this.We note at this point that the slow increase of en-tanglement entropy with system size (and with time) isprecisely what allows us to perform a computationallytractable approximation to the time evolution using Ma-trix Product States [41, 44]. SUMMARY AND PERSPECTIVE Entanglement is a fundamental feature of interactingquantum many body systems [55, 56] which, somewhatsurprisingly, can be shown to play no important role indetermining ground-state properties of systems with all-to-all interactions like for the neutrino forward-scatteringHamiltonian of Eq. (1) (see eg. [60]). This is at the heartof the expectation that a mean-field treatment of the dy-namics would also be appropriate. However, the roleplayed by entanglement in out-of-equilibrium settings,like those relevant for flavor evolution in a supernova en-vironment, is not as well understood [25, 26].In this work we have solved for the real time dynam-ics of a simpler system described by the Hamiltonian inEq. (3) sharing many similarities to the model used to de-scribe bipolar oscillations [20, 47]. In order to overcomethe exponential complexity of a direct solution of theproblem we use a simple technique, based on a MPS rep-resentation [41], whose computational complexity scaleswith the amount of entanglement generated by the dy-namics. The results for the simple model described bythe integrable Hamiltonian in Eq. (3) suggest that themaximum bipartite entropy reached in fast neutrino os-cillations scales only as log( N ) with the logarithm of thesystem size. This favorable scaling allows for the efficientcalculation of the dynamics for relatively long times andsystems exceeding 100 neutrinos. The computational ef-ficiency could also be improved by using more complextensor network states like MERA [61] which naturallydescribe states with logarithmically increasing entropy.As shown in the accompanying paper [45], the behav-ior observed here is also found in similar (non-integrable)spin systems whose out-of-equilibrium dynamics crossesa quantum critical point, suggesting a strong link be-tween dynamical phase transitions and fast collective os-cillations. This link provides us directly with conditions,like Eq. (6), to determine when instabilities can occur.These conditions provide a generalization of those com-monly obtained with linear stability analysis to the longtime, non-linear, regime.The artificial setting considered here, with a stronglyasymmetric initial state and no vacuum mixing effect,was chosen to highlight differences with a typical mean-field treatment of the problem where the dynamics issimply absent. We expect however the qualitative con-clusions reached here, fast oscillations caused by a dy-namical phase transition and a general slow growth of theentanglement entropy with system size, to be valid alsoin more general situations where both multi-angle effectsand more complicated energy distributions are consid-ered [45]. In this regime, MPS simulations could be usedreliably to track flavor dynamics up to relatively large,and more realistic, neutrino systems.Based on the results presented here, we still cannotcompletely exclude the existence of neutrino configura-tions able to generate substantially more entanglementduring the dynamics. Interestingly, the presence of theseclassically-hard configurations would be signalled clearlyby the failure of the MPS method to scale efficiently tolarge system sizes. If found, these would be ideal mod-els to explore using near-term quantum simulators liketrapped-ion systems [62] or future circuit based quantumcomputers [63], which do not necessarily suffer from anexponentially increasing computational cost when largeentanglement is present. In future studies, the effectsdiscussed here will be integrated into more realistic sim-ulations to better ascertain the role of entanglement inthe dynamical evolution of dense neutrino environmentsof astrophysical importance.I want to thank Joseph Carlson, Vincenzo Cirigliano,Huaiyu Duan, Joshua Martin, Amol Patwardhan, ErmalRrapaj and Martin Savage for the many useful discus-sions about the subject of this work. This work was sup-ported by the InQubator for Quantum Simulation underU.S. DOE grant No. DE-SC0020970 and by the Quan-tum Science Center (QSC), a National Quantum Infor- mation Science Research Center of the U.S. Departmentof Energy (DOE). [1] R. D. Hoffman, S. E. Woosley, and Y.-Z. Qian, “Nu-cleosynthesis in neutrino-driven winds. II. implicationsfor heavy element synthesis,” The Astrophysical Journal , 951–962 (1997).[2] Hans-Thomas Janka, “Explosion mechanisms of core-collapse supernovae,” Annual Review of Nuclear and Par-ticle Science , 407–451 (2012).[3] Martin J. Savage, Robert A. Malaney, and George M.Fuller, “Neutrino Oscillations and the Leptonic Chargeof the Universe,” Astrophys. J. , 1 (1991).[4] James Pantaleone, “Dirac neutrinos in dense matter,”Phys. Rev. D , 510–523 (1992).[5] James Pantaleone, “Neutrino oscillations at high densi-ties,” Physics Letters B , 128 – 132 (1992).[6] Sergio Pastor and Georg Raffelt, “Flavor oscillationsin the supernova hot bubble region: Nonlinear effectsof neutrino background,” Phys. Rev. Lett. , 191101(2002).[7] A B Balantekin and H Y¨uksel, “Neutrino mixing and nu-cleosynthesis in core-collapse supernovae,” New Journalof Physics , 51–51 (2005).[8] George M. Fuller and Yong-Zhong Qian, “Simultane-ous flavor transformation of neutrinos and antineutrinoswith dominant potentials from neutrino-neutrino forwardscattering,” Phys. Rev. D , 023004 (2006).[9] Huaiyu Duan, George M. Fuller, J. Carlson, and Yong-Zhong Qian, “Coherent development of neutrino flavor inthe supernova environment,” Phys. Rev. Lett. , 241101(2006).[10] Alexander Friedland, “Self-refraction of supernova neu-trinos: Mixed spectra and three-flavor instabilities,”Phys. Rev. Lett. , 191102 (2010).[11] Joshua D. Martin, Changhao Yi, and Huaiyu Duan,“Dynamic fast flavor oscillation waves in dense neutrinogases,” Physics Letters B , 135088 (2020).[12] Meng-Ru Wu and Irene Tamborra, “Fast neutrino con-versions: Ubiquitous in compact binary merger rem-nants,” Phys. Rev. D , 103007 (2017).[13] Gianluigi Fogli, Eligio Lisi, Antonio Marrone, andAlessandro Mirizzi, “Collective neutrino flavor transi-tions in supernovae and the role of trajectory averaging,”Journal of Cosmology and Astroparticle Physics ,010–010 (2007).[14] Yong-Zhong Qian, George M. Fuller, Grant J. Mathews,Ron W. Mayle, James R. Wilson, and S. E. Woosley,“Connection between flavor-mixing of cosmologically sig-nificant neutrinos and heavy element nucleosynthesis insupernovae,” Phys. Rev. Lett. , 1965–1968 (1993).[15] Yong-Zhong Qian and George M. Fuller, “Neutrino-neutrino scattering and matter-enhanced neutrino flavortransformation in supernovae,” Phys. Rev. D , 1479–1494 (1995).[16] R. F. Sawyer, “Speed-up of neutrino transformations ina supernova environment,” Phys. Rev. D , 045003(2005).[17] R. F. Sawyer, “Neutrino cloud instabilities just above theneutrino sphere of a supernova,” Phys. Rev. Lett. , , 042–042 (2016).[19] Ignacio Izaguirre, Georg Raffelt, and Irene Tamborra,“Fast pairwise conversion of supernova neutrinos: A dis-persion relation approach,” Phys. Rev. Lett. , 021101(2017).[20] Huaiyu Duan, George M. Fuller, and Yong-ZhongQian, “Collective neutrino oscillations,” Annual Reviewof Nuclear and Particle Science , 569–594 (2010),https://doi.org/10.1146/annurev.nucl.012809.104524.[21] Sovan Chakraborty, Rasmus Hansen, Ignacio Izaguirre,and Georg Raffelt, “Collective neutrino flavor conver-sion: Recent developments,” Nuclear Physics B , 366– 381 (2016), neutrino Oscillations: Celebrating the No-bel Prize in Physics 2015.[22] Y. Pehlivan, A. B. Balantekin, Toshitaka Kajino, andTakashi Yoshida, “Invariants of collective neutrino oscil-lations,” Phys. Rev. D , 065008 (2011).[23] Note the difference in the definition of µ there.[24] R.F. Sawyer, “’Classical’ instabilities and ’quantum’speed-up in the evolution of neutrino clouds,” (2004),arXiv:hep-ph/0408265.[25] Michael J. Cervia, Amol V. Patwardhan, A. B. Bal-antekin, S. N. Coppersmith, and Calvin W. John-son, “Entanglement and collective flavor oscillations ina dense neutrino gas,” Phys. Rev. D , 083001 (2019).[26] Ermal Rrapaj, “Exact solution of multiangle quantummany-body collective neutrino-flavor oscillations,” Phys.Rev. C , 065805 (2020).[27] Alexander Friedland and Cecilia Lunardini, “Do many-particle neutrino interactions cause a novel coherent ef-fect?” Journal of High Energy Physics , 043–043(2003).[28] Alexander Friedland, Bruce H. J. McKellar, and IvonaOkuniewicz, “Construction and analysis of a simplifiedmany-body neutrino model,” Phys. Rev. D , 093002(2006).[29] Nicole F. Bell, Andrew A. Rawlinson, and R.F. Sawyer,“Speed-up through entanglement—many-body effects inneutrino processes,” Physics Letters B , 86 – 93(2003).[30] Martin Eckstein, Marcus Kollar, and Philipp Werner,“Thermalization after an interaction quench in the hub-bard model,” Phys. Rev. Lett. , 056403 (2009).[31] Bruno Sciolla and Giulio Biroli, “Dynamical transitionsand quantum quenches in mean-field models,” Journalof Statistical Mechanics: Theory and Experiment ,P11003 (2011).[32] Bruno Sciolla and Giulio Biroli, “Quantum quenches, dy-namical transitions, and off-equilibrium quantum criti-cality,” Phys. Rev. B , 201110 (2013).[33] Since [ H, J A ] = [ H, J B ] = 0 we are free to change theintra-beam couplings at will thus changing the geome-try. In the general case we can substitute µ → µ (1 − cos( θ AB )), with θ AB the angle between the momenta inthe A and B beams, for all the results in this work.[34] H.J. Lipkin, N. Meshkov, and A.J. Glick, “Validity ofmany-body approximation methods for a solvable model:(i). exact solutions and perturbation theory,” NuclearPhysics , 188 – 198 (1965).[35] Julien Vidal, R´emy Mosseri, and Jorge Dukelsky, “En- tanglement in a first-order quantum phase transition,”Phys. Rev. A , 054101 (2004).[36] Julien Vidal, Guillaume Palacios, and Claude Aslan-gul, “Entanglement dynamics in the lipkin-meshkov-glickmodel,” Phys. Rev. A , 062304 (2004).[37] Julien Vidal, Guillaume Palacios, and R´emy Mosseri,“Entanglement in a second-order quantum phase transi-tion,” Phys. Rev. A , 022107 (2004).[38] Jos´e I. Latorre, Rom´an Or´us, Enrique Rico, and JulienVidal, “Entanglement entropy in the lipkin-meshkov-glick model,” Phys. Rev. A , 064101 (2005).[39] Pedro Ribeiro, Julien Vidal, and R´emy Mosseri, “Exactspectrum of the lipkin-meshkov-glick model in the ther-modynamic limit and finite-size corrections,” Phys. Rev.E , 021106 (2008).[40] Alexander Friedland and Cecilia Lunardini, “Neutrinoflavor conversion in a neutrino background: Single- ver-sus multi-particle description,” Phys. Rev. D , 013007(2003).[41] Guifr´e Vidal, “Efficient classical simulation of slightly en-tangled quantum computations,” Phys. Rev. Lett. ,147902 (2003).[42] Guifr´e Vidal, “Efficient simulation of one-dimensionalquantum many-body systems,” Phys. Rev. Lett. ,040502 (2004).[43] Sebastian Paeckel, Thomas K¨ohler, Andreas Swoboda,Salvatore R. Manmana, Ulrich Schollw¨ock, and ClaudiusHubig, “Time-evolution methods for matrix-productstates,” Annals of Physics , 167998 (2019).[44] Norbert Schuch, Michael M. Wolf, Frank Verstraete, andJ. Ignacio Cirac, “Entropy scaling and simulability bymatrix product states,” Phys. Rev. Lett. , 030504(2008).[45] A. Roggero, “Dynamical phase transitions in models ofcollective neutrino oscillations,” (in preparation (2021)).[46] V. Alan Kosteleck´y and Stuart Samuel, “Self-maintainedcoherent oscillations in dense neutrino gases,” Phys. Rev.D , 621–627 (1995).[47] Steen Hannestad, Georg G. Raffelt, G¨unter Sigl, andYvonne Y. Y. Wong, “Self-induced conversion in denseneutrino gases: Pendulum in flavor space,” Phys. Rev. D , 105010 (2006).[48] Huaiyu Duan, George M. Fuller, J. Carlson, and Yong-Zhong Qian, “Analysis of collective neutrino flavor trans-formation in supernovae,” Phys. Rev. D , 125005(2007).[49] Note that this condition can be modified by the presenceof an MSW-like matter term which couple differently tothe two flavors.[50] A. Sørensen, L.-M. Duan, J. I. Cirac, and P. Zoller, “En-tanglement in a second-order quantum phase transition,”Nature , 63–66 (2001).[51] G´eza T´oth, Christian Knapp, Otfried G¨uhne, andHans J. Briegel, “Optimal spin squeezing inequalities de-tect bound entanglement in spin models,” Phys. Rev.Lett. , 250405 (2007).[52] G´eza T´oth, Christian Knapp, Otfried G¨uhne, andHans J. Briegel, “Spin squeezing and entanglement,”Phys. Rev. A , 042334 (2009).[53] G´eza T´oth, “Entanglement detection in optical latticesof bosonic atoms with collective measurements,” Phys.Rev. A , 052327 (2004).[54] John K. Stockton, J. M. Geremia, Andrew C. Doherty,and Hideo Mabuchi, “Characterizing the entanglement of symmetric many-particle spin- systems,” Phys. Rev.A , 022112 (2003).[55] J. Eisert, M. Cramer, and M. B. Plenio, “Colloquium:Area laws for the entanglement entropy,” Rev. Mod.Phys. , 277–306 (2010).[56] G. Vidal, J. I. Latorre, E. Rico, and A. Kitaev, “En-tanglement in quantum critical phenomena,” Phys. Rev.Lett. , 227902 (2003).[57] G. Refael and J. E. Moore, “Entanglement entropy ofrandom quantum critical points in one dimension,” Phys.Rev. Lett. , 260602 (2004).[58] Markus Heyl, “Dynamical quantum phase transitions:a review,” Reports on Progress in Physics , 054001(2018).[59] Asmi Haldar, Krishnanand Mallayya, Markus Heyl,Frank Pollmann, Marcos Rigol, and Arnab Das, “Sig-natures of quantum phase transitions after quenchesin quantum chaotic one-dimensional systems,” arXive-prints , arXiv:2004.02905 (2020), arXiv:2004.02905[cond-mat.stat-mech].[60] Fernando G. S. L. Brand˜ao and Aram W. Har-row, “Product-state approximations to quantum states,”Communications in Mathematical Physics , 47–80(2016).[61] G. Vidal, “Class of quantum many-body states that canbe efficiently simulated,” Phys. Rev. Lett. , 110501(2008).[62] J. Zhang, G. Pagano, P. W. Hess, A. Kyprianidis,P. Becker, H. Kaplan, A. V. Gorshkov, Z.-X. Gong,and C. Monroe, “Observation of a many-body dynami-cal phase transition with a 53-qubit quantum simulator,”Nature , 601–604 (2017).[63] B. Hall, A. Roggero, A. Baroni, and J. Carlson, “Quan-tum simulation of collective neutrino oscillations,” (inpreparation (2021)).[64] Ulrich Schollw¨ock, “The density-matrix renormalizationgroup in the age of matrix product states,” Annals ofPhysics , 96 – 192 (2011), january 2011 Special Issue.[65] F. Verstraete and J. I. Cirac, “Matrix product states rep-resent ground states faithfully,” Phys. Rev. B , 094423(2006).[66] Masuo Suzuki, “General theory of fractal path integralswith applications to many-body theories and statisticalphysics,” Journal of Mathematical Physics , 400–407(1991).[67] E M Stoudenmire and Steven R White, “Minimally en-tangled typical thermal state algorithms,” New Journalof Physics , 055026 (2010).[68] Ian D. Kivlichan, Jarrod McClean, Nathan Wiebe, CraigGidney, Al´an Aspuru-Guzik, Garnet Kin-Lic Chan, andRyan Babbush, “Quantum simulation of electronic struc-ture with linear depth and connectivity,” Phys. Rev. Lett. , 110501 (2018).[69] Matthew Fishman, Steven R. White, and E. MilesStoudenmire, “The ITensor software library for tensornetwork calculations,” (2020), arXiv:2007.14822. Matrix Product State Simulation In this supplemental material we provide a self-contained explanation of the Matrix Product State(MPS) technique used to generate the results in this work. Additional information can be found in the origi-nal paper [41] and in more recent reviews [43, 64]. In or-der to implement long-range interactions efficiently witha MPS, we use the same strategy introduced recently forcircuit based simulations in [63], a brief explanation ofthe construction is presented in the last subsection. MPS Representation Consider a general quantum state | Ψ (cid:105) of a system with N spin-1 / N tensor as | Ψ (cid:105) = (cid:88) σ ,...,σ N Ψ σ ··· σ N | σ · · · σ N (cid:105) , (8)with | σ · · · σ N (cid:105) one of the 2 N basis states spanningthe total Hilbert space of the system. A Matrix Prod-uct State is an exact representation of the rank- N ten-sor Ψ σ ··· σ N as a product of N , smaller, rank-3 tensors M σ k α,β . A constructive procedure to find this representa-tion, starting from the original tensor Ψ σ ··· σ N , can beobtained by applying N times a Singular Value Decom-position (SVD) of appropriately reshaped tensors. Thisdecomposition allows to represent an arbitrary rectangu-lar matrix M of dimension D A × D B as the followingproduct of three matrices M ij = (cid:88) k U ik Σ kk V † kj , (9)where U is a D A × min( D A , D B ) matrix with orthonor-mal columns, Σ is a diagonal matrix of dimensionmin( D A , D B ) × min( D A , D B ) with either zero or posi-tive entries and V † is a min( D A , D B ) × D B matrix withorthonormal rows. An important property of the SVD,which we will use repeatedly in the derivation below, isthat the approximation of the original matrix M , withrank r ≤ min( D A , D B ), using a second matrix M (cid:48) withrank r (cid:48) < r can be done optimally, in the Frobenius norm,by truncating the sum in Eq. (9) so that we keep onlythe r (cid:48) largest singular values Σ kk . We will assume thatthe SVD is performed with the convention that singu-lar values in Σ are ordered as Σ ≥ · · · ≥ Σ NN . Withthis convention, the optimal approximation correspondsto truncating the sum as M (cid:48) ij = r (cid:48) (cid:88) k U ik Σ kk V † kj . (10)Starting from the first spin on the left, we first re-shape Ψ σ ··· σ N into the rectangular matrix Ψ σ ( σ ··· σ N ) of dimension 2 × N − . A SVD of this matrix will thenproduce the following decompositionΨ σ ( σ ··· σ N ) = r (cid:88) α U σ α Σ α α V † α ( σ ··· σ N ) = r (cid:88) α U σ α W α σ ··· σ N = r (cid:88) α A [1] σ α Ψ [1]( α σ )( σ ··· σ N ) , (11)with r ≤ V † to a vector (cid:126)W . In the last line, we have furtherreshaped the matrix U into two row vectors, with compo-nents A [1] σ α = U σ α , and the vector (cid:126)W to a 2 r × N − matrix Ψ [1]( α σ )( σ ··· σ N ) . At this point, we can repeat theprocedure by applying a SVD to this last matrix and pro-ceed from left to right until we have reached the last spinof the chain. The resulting decomposition becomesΨ σ ··· σ N = (cid:88) α ,...,α N − A [1] σ α A [2] σ α α · · · A [ N − σ N − α N − α N − A [ N ] σ N α N − = A [1] σ A [2] σ · · · A [ N − σ N − A [ N ] σ N , (12)where in the last line we suppressed the lower indices andreplaced them by matrix multiplications. Using this con-struction, a general state | Ψ (cid:105) can be written as a MatrixProduct State. Note that the rank of the SVD in the cen-ter of the chain could be as large as 2 N/ − (assuming N even) and storage requirements for representing a generalstate 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 isrelated to the maximum bipartite entanglement entropyin the system. In order to see this, it is first useful tointroduce a different construction for the MPS. Considera similar iterative decomposition procedure for the ten-sor Ψ σ ··· σ N as the one described above, but now startingfrom the right-most site. The first step leads toΨ ( σ ··· σ N − ) σ N = r N (cid:88) β N U ( σ ··· σ N − ) β N Σ β N β N V † β N σ N = r N (cid:88) β N W (cid:48) σ ··· σ N − β N V † β N σ N = r N (cid:88) β N Ψ [ N ]( σ ··· σ N − )( σ N − β N ) B [ N ] σ N β N , (13)where we have introduced the two column vectors B [ N ] σ N β N = V † β N σ N and the residual matrix has maximumdimensions given by 2 N − × r N . Proceeding with thesame procedure until the left edge, we find another ma-trix product state with amplitudesΨ σ ··· σ N = B [1] σ B [2] σ · · · B [ N − σ N − B [ N ] σ N . (14) The two representations are related by a gauge transfor-mation [43, 64], and they are characterized by differentnormalization properties. In particular, we have (cid:88) σ k A [ k ] σ k † A [ k ] σ k = (cid:88) σ k B [ k ] σ k B [ k ] σ k † = . (15)For this reason, the state resulting from the decompo-sition in Eq. (12) in terms of A tensors is called a left-canonical MPS, while the construction in terms of B ten-sors in Eq. (14) is called right-canonical [43, 64].In order to see more clearly the connection betweenentanglement and the maximum rank in the decomposi-tion, also called bond dimension , it is useful to introduceyet another canonical form. This is obtained by startingwith a decomposition from the left up to a chosen site k , at this point we start with a decomposition from theright and we stop when we reach site k + 1. We can writethe result of this decomposition asΨ σ ··· σ N = A [1] σ · · · A [ k ] σ k Σ [ k ] B [ k +1] σ k +1 · · · B [ N ] σ N , (16)with the diagonal matrix Σ [ k ] containing the singular val-ues at the bond between the k and the k + 1 sites. Thesite k is called the orthogonality center , since tensors tothe left are left-normalized while tensors on the right areright-normalized. This property allows us to write thestate in this representation as | Ψ (cid:105) = r k (cid:88) α k Σ [ k ] α k α k (cid:12)(cid:12) Φ Lα k (cid:11) ⊗ (cid:12)(cid:12) Φ Rα k (cid:11) , (17)where we have defined two sets of states (cid:12)(cid:12) Φ Lα k (cid:11) = (cid:88) σ ,...,σ k (cid:104) A [1] σ · · · A [ k ] σ k (cid:105) α k | σ · · · σ k (cid:105) , (18)and (cid:12)(cid:12) Φ Rα k (cid:11) = (cid:88) σ k +1 ,...,σ N (cid:104) B [ k +1] σ k +1 · · · B [ N ] σ N (cid:105) α k | σ k +1 · · · σ N (cid:105) . (19)Thanks to the normalization properties of the A and B tensors from Eq. (15), we see that these states form or-thonormal sets of states for the left and right Hilbertspaces respectively. We can then interpret Eq. (17) asthe Schmidt decomposition of the state | Ψ (cid:105) on the left-right bipartition. The number of states in the sum canthen be identified with the Schmidt rank , which is a nat-ural measure of entanglement between states on the leftpartition and states on the right one.In particular, the Von-Neumann entanglement entropyfor the left-right bipartition at site k of the state | Ψ (cid:105) canbe expressed as a function of the singular values as S k = − r k (cid:88) α k =1 Σ [ k ]2 α k α k log (cid:16) Σ [ k ]2 α k α k (cid:17) . (20)The Schmidt rank r k gives the number of non-zero sin-gular values so that the entanglement is bounded by S k ≤ log ( r k ) . (21)We can now see that MPS with small bond dimen-sions correspond to quantum states with little entan-glement [41]. Before moving on to the time-dependentsimulation, it’s important to mention that the mixed-canonical form in Eq. (16) is also particularly convenientto find more efficient approximations to an initial MPS.Given an MPS | Ψ (cid:105) with bond dimension r k at the bondbetween site k and k + 1, we can find another MPS | Ψ (cid:48) (cid:105) with lower bond dimension r (cid:48) k < r k so that the error in2-norm is bounded from above by [65] (cid:107)| Ψ (cid:105) − | Ψ (cid:48) (cid:105)(cid:107) ≤ r k (cid:88) α k = r (cid:48) k +1 Σ [ k ]2 α k α k . (22)Thanks to the bi-orthogonality of the Schmidt decompo-sition Eq. (17), we can obtain | Ψ (cid:48) (cid:105) from | Ψ (cid:105) by keepingonly the r (cid:48) k largest singular values in the sum. If thesingular values decrease sufficiently fast, a good approx-imation to the original state could be obtained with anMPS with a much smaller bond dimension [41]. Time evolution of a MPS We now turn to explain how one can simulate the time-evolution for the general neutrino Hamiltonian in Eq. (1)of the main text. The basic operation we will use repeat-edly in the following is the application of an SU (4) oper-ation on a pair of nearest neighbor indices ( k, k + 1). Ageneric 4 × U of this type can be represented asa rank 4 tensor G σ σ σ (cid:48) σ (cid:48) with each index having dimensionequal to 2. When we contract this tensor with the MPS,only the tensors around the bond ( k, k + 1) get trans-formed. More specifically, we have that the matrices C σ k σ k +1 = A [ k ] σ k Σ [ k ] B [ k +1] σ k +1 , (23)get updated into (cid:101) C σ k σ k +1 = (cid:88) σ (cid:48) k σ (cid:48) k +1 G σ k σ k +1 σ (cid:48) k σ (cid:48) k +1 A [ k ] σ (cid:48) k Σ [ k ] B [ k +1] σ (cid:48) k +1 . (24)We can now restore the MPS form by performing onefinal SVD of (cid:101) C σ k σ k +1 properly reashaped as a matrix (cid:101) C ( α k − σ k )( σ k +1 β k +2 ) = (cid:101) r k (cid:88) α k U ( α k − σ k ) α k Σ α k α k V † α k ( σ k +1 β k +2 ) , (25)with a possibly larger bond dimension (cid:101) r k . Now, by identi-fying (cid:101) A [ k ] σ k with U , (cid:101) B [ k +1] σ k +1 with V † and the singular vector matrices (cid:101) Σ [ k ] with Σ, we have a new MPS withdecomposition (cid:101) Ψ σ ··· σ N = A [1] σ · · · (cid:101) A [ k ] σ k (cid:101) Σ [ k ] (cid:101) B [ k +1] σ k +1 · · · B [ N ] σ N . (26)The SVD requires O (( r k − + r k +1 ) r k +1 ) operations andis efficient only if we are able to keep a small bond di-mension by truncating the MPS as shown above.With nearest-neighbor two-site operators at our dis-posal, we can now present the scheme used to simulatetime evolution under Hamiltonians of the form H = N (cid:88) k =1 (cid:126)B k · (cid:126)σ k + N (cid:88) i 1) gates. In order to replace longrange gates into nearest-neighbor ones, it is convenientto introduce the SWAP operation. This is defined asSWAP = . (30)When applied to a pair of sites, it exchanges the state asSWAP | φ (cid:105) ⊗ | ψ (cid:105) = | ψ (cid:105) ⊗ | φ (cid:105) . (31)Using this operation we can transport two sites next toeach other, apply the unitary T ij = exp( − iδh ij ) as anearest-neighbor operation and finally swap the two sitesback to their original position [67]. Given the all-to-allnature of the pair interaction in Eq. (28), a naive applica-tion of this strategy will result in O ( N ) nearest-neighboroperations to implement one step of time evolution usingthe approximation in Eq. (29) for the propagator.In this work we employ instead the swap-networkscheme introduced in Ref. [63] for studying neutrino dy-namics on a digital quantum computer with linear con-nectivity. The construction is inspired by earlier work0 FIG. 4. (Color online) Schematic representation of one step of time evolution for a system with N = 4 neutrinos. The label onthe top indicate the Hamiltonian terms h ij considered at each one of the 2 N = 8 sub-steps. on fermionic swap-networks [68] and relies on the obser-vation that, if we perform N layers of nearest-neighborswap operations by alternating the even bonds with theodd bonds, we can bring every site next to every othersite at least once using a total of N ( N − / N steps, the sites will be left in re-verse order. It is now clear that, if we define a modifiedpropagator S ij as follows S ij = SWAP × exp( − iδh ij ) , (32)and apply nearest-neighbor S ij operations instead ofswaps as before, at the end of the sequence we will haveimplemented the first product of pair operators in theapproximation Eq. (29) for the propagator. If we now ap-ply the sequence in reverse order we can complete a fulltime step and obtain a state where sites have the originalorder. A schematic depiction of a complete step for a sys-tem with N = 4 neutrinos is displayed in Fig. 4. In thiswork we perform a SVD truncation to MPS form aftereach application of the nearest neighbor operation usinga threshold of 10 for the discarded eigenvalues. Bettertruncation schemes are possible (see eg. [43]) and their potential advantages will be explored in future work.Finally, we decompose the full time evolution for a to-tal time t into steps of size δ (cid:28) t asexp ( − itH ) = [exp ( − iδH )] t/δ . (33)We then approximate each short-time propagator fortime δ using the swap-network scheme described above.For all the results presented in the main text, we useda fixed time-step δ = 0 . µ − which we found provideda good compromise between precision and efficiency forthe systems studied here. As already mentioned above,the overall scheme remains efficient provided the singularvalues on each bond decay sufficiently fast to zero that agood approximation to the MPS can be obtained usinga low bond dimension. Thanks to the (at most) loga-rithmic scaling of the entanglement entropy observed inthe systems studied in main text, the maximum bonddimension required to keep the truncation error below10 − was ≈ N/ N (cid:39)(cid:39)