# 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

###### Abstract

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 Hoffman  _et al._ (1997); Janka
(2012); Winteler  _et al._ (2012); Wanajo  _et al._ (2014). In these
environments with a large neutrino density $\rho_{\nu}$, neutrino flavor
evolution is substantially modified by neutrino-neutrino scattering processes
which can lead to self-sustained collective flavor oscillations Savage  _et
al._ (1991); Pantaleone (1992a, b); Pastor and Raffelt (2002); Balantekin and
Yüksel (2005); Fuller and Qian (2006); Duan  _et al._ (2006a); Friedland
(2010); Martin  _et al._ (2020); Wu and Tamborra (2017). Since neutrinos in
supernovae are emitted with fluxes and spectra that are strongly flavor
dependent Janka (2012), the presence of collective flavor oscillations could
then lead to important effects Fogli  _et al._ (2007); Qian  _et al._ (1993);
Qian and Fuller (1995). In some situations however (see eg. Duan and Friedland
(2011); Dasgupta  _et al._ (2012)) these could appear only at too large radii
to modify the explosion mechanism in supernovae.

For large neutrino densities near the proto-neutronstar, flavor oscillations
can occur at frequencies of the order of the neutrino-neutrino interaction
energy $\mu=\sqrt{2}G_{F}\rho_{\nu}$ instead of the much smaller vacuum
oscillation frequency $\omega=\Delta m^{2}/(2E)$, a phenomenon known as fast
flavor conversion Sawyer (2005, 2016); Chakraborty  _et al._ (2016a);
Izaguirre  _et al._ (2017) (see Chakraborty  _et al._ (2016b) for a review).
In these expressions $G_{F}$ is Fermi’s constant, $E$ the neutrino energy and
$\Delta m^{2}=m_{2}^{2}-m_{1}^{2}$ 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 forward-scattering neutrino-neutrino interactions where only
flavor can be exchanged among neutrinos (for the full description see eg.
Vlasenko  _et al._ (2014); Cirigliano  _et al._ (2015)). In this scheme, a
system with $N$ neutrinos is represented as a collection of $SU(2)$ flavor
isospins governed by the Hamiltonian Pehlivan  _et al._ (2011) 111Note the
difference in the definition of $\mu$ there.

$H_{fs}=\sum_{k=1}^{N}\frac{\Delta
m^{2}}{4E_{k}}\vec{B}\cdot\vec{\sigma}_{k}+\frac{\mu}{2N}\sum_{i<j}^{N}\mathcal{J}_{ij}\vec{\sigma}_{i}\cdot\vec{\sigma}_{j}\;,$
(1)

with $\vec{\sigma}_{i}=(\sigma^{x}_{i},\sigma^{y}_{i},\sigma^{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
$\\{\ket{\uparrow},\ket{\downarrow}\\}=\\{\ket{\nu_{e}},\ket{\nu_{x}}\\}$, the
orientation of the (normalized) vector
$\vec{B}=(\sin(2\theta),0,-\cos(2\theta))$ is related to the mixing angle
$\theta$. The geometry of the problem is encoded in the two-body coupling
matrix $\mathcal{J}_{ij}=\left(1-\hat{p}_{i}\cdot\hat{p}_{k}\right)$ with
$\hat{p}_{k}=\vec{p}_{k}/|\vec{p}_{k}|$ and $\vec{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.
Stirner  _et al._ (2018)). This Hamiltonian model is the base of all the past
studies of (two-flavor) collective neutrino oscillations in anisotropic but
homogeneous settings (see Pehlivan  _et al._ (2011) 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 Sawyer
(2004); Duan  _et al._ (2010); Izaguirre  _et al._ (2017). The role of
entanglement and quantum effects in the out-of-equilibrium dynamics Eisert
_et al._ (2015) of neutrinos has received renewed interest recently Cervia
_et al._ (2019); Rrapaj (2020) 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 Friedland and Lunardini
(2003a); Friedland  _et al._ (2006) or substantial flavor evolution over time
scales $\tau_{F}=\mu^{-1}\log(N)$, which can remain relevant for large systems
Bell  _et al._ (2003); Sawyer (2004). We will refer to oscillations on the
time scale $\tau_{F}$ to be “fast” in contrast to “slow” oscillations
occurring for $\tau_{L}=\mu^{-1}\sqrt{N}$. This convention is different from
the more common definition in the literature where “fast” and “slow”
oscillations are associated with time scales $\approx\mu^{-1}$ and
$\approx\sqrt{\mu\omega}$ 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 at the fast scale $\tau_{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) Vidal (2003)
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. Roggero (2021).
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 Izaguirre  _et al._ (2017) 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
$\ket{\Psi_{0}}=\left(\bigotimes_{n=1}^{N/2}\ket{\downarrow}\right)\otimes\left(\bigotimes_{m=1}^{N/2}\ket{\uparrow}\right)$.
This corresponds to choosing all the neutrinos in the $A$ beam to be
$\ket{\nu_{2}}$ and all the neutrinos in the $B$ beam to be $\ket{\nu_{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 $\mathcal{J}_{ij}$ and take the single angle approximation by setting
$\mathcal{J}_{ij}=1$ for all neutrinos Duan  _et al._ (2010); Pehlivan  _et
al._ (2011) (see also 222Since $[H,J^{2}_{A}]=[H,J^{2}_{B}]=0$ we are free to
change the intra-beam couplings at will thus changing the geometry. In the
general case we can substitute $\mu\to\mu(1-\cos(\theta_{AB}))$, with
$\theta_{AB}$ the angle between the momenta in the $A$ and $B$ beams, for all
the results in this work.). 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.
Cervia  _et al._ (2019), can be written as follows

$\begin{split}H&=-\frac{\omega_{A}}{2}\sum_{i\in\mathcal{A}}\sigma_{i}^{z}-\frac{\omega_{B}}{2}\sum_{i\in\mathcal{B}}\sigma_{i}^{z}+\frac{\mu}{2N}\sum_{i<j}\vec{\sigma}_{i}\cdot\vec{\sigma}_{j}\;,\end{split}$
(2)

where $\mathcal{A}$ and $\mathcal{B}$ are sets of indices for the neutrinos in
the $A$ and $B$ beam respectively while we set the one body energies as
$\omega_{k}=\Delta 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 Pehlivan  _et al._ (2011). This was used in Ref. Cervia
_et al._ (2019) to study a time-dependent generalization of the model for
small systems.

At this point it is convenient to introduce the total spin operator
$\vec{J}=\sum_{i}\vec{\sigma}_{i}/2$ and correspondingly $\vec{J}^{A}$ and
$\vec{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_{z}^{2}$. As the initial state $\ket{\Psi_{0}}$ is an eigenstate
of $J_{z}$, we can rewrite the relevant Hamiltonian for our setup as follows

$H=-\delta_{\omega}\left(J_{z}^{A}-J_{z}^{B}\right)+\frac{\mu}{N}J^{2}\;,$ (3)

where we dropped irrelevant constant factors and introduced a single parameter
$\delta_{\omega}=(\omega_{A}-\omega_{B})/2$ for the energy difference between
the two beams. This model is a variant of the Lipkin-Meshov-Glick (LMG) model
Lipkin  _et al._ (1965) explored extensively in the past Vidal  _et al._
(2004a, b, c); Latorre  _et al._ (2005); Ribeiro  _et al._ (2008). 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. Friedland and Lunardini (2003a); Friedland  _et al._ (2006). In
this model, collective flavor oscillation are shown to occur at time scales
much longer than $\tau_{F}$ with the scaling $\tau_{L}=\mu^{-1}\sqrt{N}$
expected from incoherent scattering from the background neutrinos Friedland
and Lunardini (2003a, b). Due to this strong scaling with system size, when
$N\gg 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 Friedland and Lunardini
(2003b, a).

As a way to track flavor oscillations, we compute the survival probability
$p(t)=(1-\langle\sigma^{z}_{1}(t)\rangle)/2$ for a neutrino (we take the
first) to be found in it initial flavor state (in this case
$\ket{\downarrow}$). By using the exchange symmetry for neutrinos within each
beam of both the Hamiltonian and the initial state $\ket{\Psi_{0}}$, together
with energy conservation, we can write the expectation value of the total spin
as

$\langle
J^{2}(t)\rangle=\frac{N}{2}\left(1+2\frac{\delta_{\omega}}{\mu}N\left(1-p(t)\right)\right)\;.$
(4)

When $\delta_{\omega}/\mu<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)\geq
1-\mu/(2N\delta_{\omega})$ and the neutrinos remain frozen in the initial
state in the thermodynamic limit. In the opposite limit, for
$\delta_{\omega}/\mu\geq 1/4$ flavor oscillations become parametrically small
as $p(t)\geq 1-\mu/(4\delta_{\omega})$ but remain finite in the $N\gg 1$
limit. In the intermediate regime $0<\delta_{\omega}/\mu\lesssim 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 Vidal (2003, 2004); Paeckel
_et al._ (2019). The approach has a controllable error and it’s
computationally efficient for systems with low levels of bipartite
entanglement Vidal (2003); Schuch  _et al._ (2008) 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 Vidal (2003) and in more recent reviews Schollwöck (2011);
Paeckel  _et al._ (2019). In order to implement long-range interactions
efficiently with a MPS, we use the same strategy introduced recently for
circuit based simulations in Hall  _et al._ (2021), a brief explanation of the
construction is presented in the last subsection. In the last section we also
compare directly the mean-field approximation to an MPS-based simulation with
unit bond dimension.

###  II.1 MPS Representation

Consider a general quantum state $\ket{\Psi}$ of a system with $N$ spin-$1/2$
particle, it’s wavefunction can be expressed as a $N$-component tensor as

$\ket{\Psi}=\sum_{\sigma_{1},\dots,\sigma_{N}}\Psi^{\sigma_{1}\cdots\sigma_{N}}\ket{\sigma_{1}\cdots\sigma_{N}}\;,$
(5)

with $\ket{\sigma_{1}\cdots\sigma_{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
$\Psi^{\sigma_{1}\cdots\sigma_{N}}$ as a product of $N$, smaller,
$3$-component tensors $A^{\sigma_{k}}_{\alpha,\beta}$. A constructive
procedure to find this representation, starting from the original tensor
$\Psi^{\sigma_{1}\cdots\sigma_{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}\times D_{B}$ as the following product of three matrices

$M_{ij}=\sum_{k}U_{ik}\Sigma_{kk}V^{\dagger}_{kj}\;,$ (6)

where $U$ is a $D_{A}\times\min(D_{A},D_{B})$ matrix with orthonormal columns,
$\Sigma$ is a diagonal matrix of dimension
$\min(D_{A},D_{B})\times\min(D_{A},D_{B})$ with either zero or positive
entries and $V^{\dagger}$ is a $\min(D_{A},D_{B})\times 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\leq\min(D_{A},D_{B})$, using a second matrix
$M^{\prime}$ with rank $r^{\prime}<r$ can be done optimally, in the Frobenius
norm, by truncating the sum in Eq. (6) so that we keep only the $r^{\prime}$
largest singular values $\Sigma_{kk}$. We will assume that the SVD is
performed with the convention that singular values in $\Sigma$ are ordered as
$\Sigma_{11}\geq\cdots\geq\Sigma_{NN}$. With this convention, the optimal
approximation corresponds to truncating the sum as

$M^{\prime}_{ij}=\sum_{k}^{r^{\prime}}U_{ik}\Sigma_{kk}V^{\dagger}_{kj}\;.$
(7)

Starting from the first spin on the left, we first reshape
$\Psi^{\sigma_{1}\cdots\sigma_{N}}$ into the rectangular matrix
$\Psi_{\sigma_{1}(\sigma_{2}\cdots\sigma_{N})}$ of dimension $2\times
2^{N-1}$. An explicit SVD of this matrix using Eq. (6) will then produce the
following decomposition

$\begin{split}\Psi_{\sigma_{1}(\sigma_{2}\cdots\sigma_{N})}&=\sum_{\alpha_{1}}^{r_{1}}U_{\sigma_{1}\alpha_{1}}\Sigma_{\alpha_{1}\alpha_{1}}V^{\dagger}_{\alpha_{1}(\sigma_{2}\cdots\sigma_{N})}\\\
&=\sum_{\alpha_{1}}^{r_{1}}U_{\sigma_{1}\alpha_{1}}W_{\alpha_{1}\sigma_{2}\cdots\sigma_{N}}\\\
&=\sum_{\alpha_{1}}^{r_{1}}A^{[1]\sigma_{1}}_{\alpha_{1}}\Psi^{[1]}_{(\alpha_{1}\sigma_{2})(\sigma_{3}\cdots\sigma_{N})}\;,\end{split}$
(8)

with $r_{1}\leq 2$ the rank of the original rectangular matrix and in the
second line we have reshaped the product of $\Sigma$ and $V^{\dagger}$ to a
vector $\vec{W}$. In the last line, we have further reshaped the matrix $U$
into two row vectors, with components
$A^{[1]\sigma_{1}}_{\alpha_{1}}=U_{\sigma_{1}\alpha_{1}}$, and the vector
$\vec{W}$ to a $2r_{1}\times 2^{N-2}$ matrix
$\Psi^{[1]}_{(\alpha_{1}\sigma_{2})(\sigma_{3}\cdots\sigma_{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

$\begin{split}\Psi^{\sigma_{1}\cdots\sigma_{N}}&=\sum_{\alpha_{1},\dots,\alpha_{N-1}}A^{[1]\sigma_{1}}_{\alpha_{1}}A^{[2]\sigma_{2}}_{\alpha_{1}\alpha_{2}}\cdots
A^{[N-1]\sigma_{N-1}}_{\alpha_{N-2}\alpha_{N-1}}A^{[N]\sigma_{N}}_{\alpha_{N-1}}\\\
&=\mathbf{A}^{[1]\sigma_{1}}\mathbf{A}^{[2]\sigma_{2}}\cdots\mathbf{A}^{[N-1]\sigma_{N-1}}\mathbf{A}^{[N]\sigma_{N}}\;,\end{split}$
(9)

where in the last line we suppressed the lower indices and replaced them by
matrix multiplications. Using this construction, a general state $\ket{\Psi}$
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
$\Psi^{\sigma_{1}\cdots\sigma_{N}}$ as the one described above, but now
starting from the right-most site. The first step leads to

$\begin{split}\Psi_{(\sigma_{1}\cdots\sigma_{N-1})\sigma_{N}}&=\sum_{\beta_{N}}^{r_{N}}U_{(\sigma_{1}\cdots\sigma_{N-1})\beta_{N}}\Sigma_{\beta_{N}\beta_{N}}V^{\dagger}_{\beta_{N}\sigma_{N}}\\\
&=\sum_{\beta_{N}}^{r_{N}}W^{\prime}_{\sigma_{1}\cdots\sigma_{N-1}\beta_{N}}V^{\dagger}_{\beta_{N}\sigma_{N}}\\\
&=\sum_{\beta_{N}}^{r_{N}}\Psi^{[N]}_{(\sigma_{1}\cdots\sigma_{N-2})(\sigma_{N-1}\beta_{N})}B^{[N]\sigma_{N}}_{\beta_{N}}\;,\end{split}$
(10)

where we have introduced the two column vectors
$B^{[N]\sigma_{N}}_{\beta_{N}}=V^{\dagger}_{\beta_{N}\sigma_{N}}$ and the
residual matrix has maximum dimensions given by $2^{N-2}\times 2r_{N}$.
Proceeding with the same procedure until the left edge, we find another matrix
product state with amplitudes

$\begin{split}\Psi^{\sigma_{1}\cdots\sigma_{N}}&=\mathbf{B}^{[1]\sigma_{1}}\mathbf{B}^{[2]\sigma_{2}}\cdots\mathbf{B}^{[N-1]\sigma_{N-1}}\mathbf{B}^{[N]\sigma_{N}}\;.\end{split}$
(11)

The two representations are related by a gauge transformation Schollwöck
(2011); Paeckel  _et al._ (2019), and they are characterized by different
normalization properties. In particular, we have

$\sum_{\sigma_{k}}\mathbf{A}^{[k]\sigma_{k}\dagger}\mathbf{A}^{[k]\sigma_{k}}=\mathbb{1}\quad\sum_{\sigma_{k}}\mathbf{B}^{[k]\sigma_{k}}\mathbf{B}^{[k]\sigma_{k}\dagger}=\mathbb{1}\;.$
(12)

For this reason, the state resulting from the decomposition in Eq. (9) in
terms of $\mathbf{A}$ tensors is called a left-canonical MPS, while the
construction in terms of $\mathbf{B}$ tensors in Eq. (11) is called right-
canonical Schollwöck (2011); Paeckel  _et al._ (2019).

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

$\begin{split}\Psi^{\sigma_{1}\cdots\sigma_{N}}&=\mathbf{A}^{[1]\sigma_{1}}\cdots\mathbf{A}^{[k]\sigma_{k}}\mathbf{\Sigma}^{[k]}\mathbf{B}^{[k+1]\sigma_{k+1}}\cdots\mathbf{B}^{[N]\sigma_{N}}\;,\end{split}$
(13)

with the diagonal matrix $\mathbf{\Sigma}^{[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

$\ket{\Psi}=\sum_{\alpha_{k}}^{r_{k}}\Sigma^{[k]}_{{\alpha_{k}}{\alpha_{k}}}\ket{\Phi^{L}_{\alpha_{k}}}\otimes\ket{\Phi^{R}_{\alpha_{k}}}\;,$
(14)

where we have defined two sets of states

$\ket{\Phi^{L}_{\alpha_{k}}}=\sum_{\sigma_{1},\dots,\sigma_{k}}\left[\mathbf{A}^{[1]\sigma_{1}}\cdots\mathbf{A}^{[k]\sigma_{k}}\right]_{\alpha_{k}}\ket{\sigma_{1}\cdots\sigma_{k}}\;,$
(15)

and

$\ket{\Phi^{R}_{\alpha_{k}}}=\sum_{\sigma_{k+1},\dots,\sigma_{N}}\left[\mathbf{B}^{[k+1]\sigma_{k+1}}\cdots\mathbf{B}^{[N]\sigma_{N}}\right]_{\alpha_{k}}\ket{\sigma_{k+1}\cdots\sigma_{N}}\;.$
(16)

Thanks to the normalization properties of the $\mathbf{A}$ and $\mathbf{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 $\ket{\Psi}$ on
the left-right 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 $\ket{\Psi}$ can be expressed as a
function of the singular values as

$S_{k}=-\sum_{\alpha_{k}=1}^{r_{k}}\Sigma^{[k]2}_{{\alpha_{k}}{\alpha_{k}}}\log_{2}\left(\Sigma^{[k]2}_{{\alpha_{k}}{\alpha_{k}}}\right)\;.$
(17)

The Schmidt rank $r_{k}$ gives the number of non-zero singular values so that
the entanglement is bounded by

$S_{k}\leq\log_{2}(r_{k})\;.$ (18)

We can now see that MPS with small bond dimensions correspond to quantum
states with little entanglement Vidal (2003). Before moving on to the time-
dependent simulation, it’s important to mention that the mixed-canonical form
in Eq. (13) is also particularly convenient to find more efficient
approximations to an initial MPS. Given an MPS $\ket{\Psi}$ with bond
dimension $r_{k}$ at the bond between site $k$ and $k+1$, we can find another
MPS $\ket{\Psi^{\prime}}$ with lower bond dimension $r^{\prime}_{k}<r_{k}$ so
that the error in $2$-norm is bounded from above by Verstraete and Cirac
(2006)

$\left\|\ket{\Psi}-\ket{\Psi^{\prime}}\right\|^{2}_{2}\leq
2\sum_{\alpha_{k}=r_{k}^{\prime}+1}^{r_{k}}\Sigma^{[k]2}_{{\alpha_{k}}{\alpha_{k}}}\;.$
(19)

Thanks to the bi-orthogonality of the Schmidt decomposition Eq. (14), we can
obtain $\ket{\Psi^{\prime}}$ from $\ket{\Psi}$ by keeping only the
$r_{k}^{\prime}$ 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 Vidal (2003).

###  II.2 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 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\times 4$ unitary
$U$ of this type can be represented as a $4$-component tensor
$G^{\sigma_{1}\sigma_{2}}_{\sigma_{1}^{\prime}\sigma_{2}^{\prime}}$ 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

$\mathbf{C}^{\sigma_{k}\sigma_{k+1}}=\mathbf{A}^{[k]\sigma_{k}}\mathbf{\Sigma}^{[k]}\mathbf{B}^{[k+1]\sigma_{k+1}}\;,$
(20)

get updated into

$\begin{split}\widetilde{\mathbf{C}}^{\sigma_{k}\sigma_{k+1}}&=\sum_{\sigma_{k}^{\prime}\sigma_{k+1}^{\prime}}G^{\sigma_{k}\sigma_{k+1}}_{\sigma_{k}^{\prime}\sigma_{k+1}^{\prime}}\mathbf{A}^{[k]\sigma_{k}^{\prime}}\mathbf{\Sigma}^{[k]}\mathbf{B}^{[k+1]\sigma_{k+1}^{\prime}}\;.\end{split}$
(21)

We can now restore the MPS form by performing one final SVD of
$\widetilde{\mathbf{C}}^{\sigma_{k}\sigma_{k+1}}$ properly reshaped as a
matrix

$\widetilde{C}_{(\alpha_{k-1}\sigma_{k})(\sigma_{k+1}\beta_{k+2})}=\sum_{\alpha_{k}}^{\widetilde{r}_{k}}\widetilde{U}_{(\alpha_{k-1}\sigma_{k})\alpha_{k}}\widetilde{\Sigma}_{\alpha_{k}\alpha_{k}}\widetilde{V}^{\dagger}_{\alpha_{k}(\sigma_{k+1}\beta_{k+2})}\;,$
(22)

with a possibly larger bond dimension $\widetilde{r}_{k}$. Now, by identifying
$\widetilde{\mathbf{A}}^{[k]\sigma_{k}}$ with $\widetilde{U}$,
$\widetilde{\mathbf{B}}^{[k+1]\sigma_{k+1}}$ with $\widetilde{V}^{\dagger}$
and the singular vector matrices $\widetilde{\mathbf{\Sigma}}^{[k]}$ with
$\widetilde{\Sigma}$, we have a new MPS with decomposition

$\begin{split}\widetilde{\Psi}^{\sigma_{1}\cdots\sigma_{N}}&=\mathbf{A}^{[1]\sigma_{1}}\cdots\widetilde{\mathbf{A}}^{[k]\sigma_{k}}\widetilde{\mathbf{\Sigma}}^{[k]}\widetilde{\mathbf{B}}^{[k+1]\sigma_{k+1}}\cdots\mathbf{B}^{[N]\sigma_{N}}\;.\end{split}$
(23)

The SVD requires $\mathcal{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.

Figure 1: (Color online) Schematic representation of one step of time
evolution for a system with $N=4$ neutrinos. The label on the top indicate the
Hamiltonian terms $h_{ij}$ considered at each one of the $2N=8$ sub-steps.

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

$H=\sum_{k=1}^{N}\vec{B}_{k}\cdot\vec{\sigma}_{k}+\sum_{i<j}^{N}\mathcal{V}_{ij}\vec{\sigma}_{i}\cdot\vec{\sigma}_{j}\;,$
(24)

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

$H=\sum_{i<j}^{N}\left[\frac{\vec{B}_{i}\cdot\vec{\sigma}_{i}+\vec{B}_{j}\cdot\vec{\sigma}_{j}}{N-1}+\mathcal{V}_{ij}\vec{\sigma}_{i}\cdot\vec{\sigma}_{j}\right]=\sum_{i<j}^{N}h_{ij}.$
(25)

The time evolution operator under this Hamiltonian can be approximated for
short times using

$e^{-i\delta H}=\prod_{i<j=1}^{N}e^{-i\delta
h_{ij}}\prod^{1}_{i<j=N}e^{-i\delta
h_{ij}}+\mathcal{O}\left(\delta^{3}\|H\|^{3}\right)\;,$ (26)

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 Suzuki (1991). 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

$\text{SWAP}=\begin{pmatrix}1&0&0&0\\\ 0&0&1&0\\\ 0&1&0&0\\\ 0&0&0&1\\\
\end{pmatrix}\;.$ (27)

When applied to a pair of sites, it exchanges the state as

$\text{SWAP}\ket{\phi}\otimes\ket{\psi}=\ket{\psi}\otimes\ket{\phi}\;.$ (28)

Using this operation we can transport two sites next to each other, apply the
unitary $T_{ij}=\exp(-i\delta h_{ij})$ as a nearest-neighbor operation and
finally swap the two sites back to their original position Stoudenmire and
White (2010). Given the all-to-all nature of the pair interaction in Eq. (25),
a naive application of this strategy will result in $\mathcal{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. Hall
_et al._ (2021) for studying neutrino dynamics on a digital quantum computer
with linear connectivity. The construction is inspired by earlier work on
fermionic swap-networks Kivlichan  _et al._ (2018) 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))

$S_{ij}=\text{SWAP}\times\exp(-i\delta h_{ij})\;,$ (29)

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

$\begin{split}Q^{e}_{k}&=\prod_{a=1}^{N/2}S^{2a-1,2a}_{i(2a-1,k,e),j(2a,k,e)}\\\
Q^{o}_{k}&=\prod_{a=1}^{N/2-1}S^{2a,2a+1}_{i(2a,k,o),j(2a+1,k,o)}\;,\end{split}$
(30)

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

$\prod_{i<j=1}^{N}e^{-i\delta h_{ij}}\prod^{1}_{i<j=N}e^{-i\delta
h_{ij}}=\prod_{k=1}^{N/2}\left[Q^{e}_{k}Q^{o}_{k}\right]\prod_{k=N/2}^{1}\left[Q^{o}_{k}Q^{e}_{k}\right]\;.$
(31)

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.
Paeckel  _et al._ (2019)) 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 $\delta\ll t$ as

$\exp\left(-itH\right)=\left[\exp\left(-i\delta H\right)\right]^{t/\delta}\;.$
(32)

We then approximate each short-time propagator for time $\delta$ using the
swap-network scheme described above.

For all the results presented in the main text, we used a fixed time-step
$\delta=0.05\mu^{-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
$\approx N/2$ allowing us to easily explore systems with $N\gtrapprox 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 Fishman  _et al._ (2020).

##  III Results

Figure 2: (Color online) Flavor polarization per particle $\langle
J^{A}_{z}(t)\rangle/(N/4)$ of neutrinos in the A beam as a function of time
for six values of the energy asymmetry parameter $\delta_{\omega}/\mu$ (from
top to bottom): $-0.5,0.0,0.125,0.25,0.5,1.0$.  Figure 3: (Color online) Panel
(a) shows the survival probability $p(t)$ in a system of $N=96$ neutrinos
initialized in $\ket{\Psi_{0}}$ and evolved with different values of the one
body asymmetry parameter. Panel (b) shows the evolution of the minimum time
$t_{P}$ as a function of system size and panel (c) shows the (normalized)
expectation value of the total angular momentum for different values of
$\delta_{\omega}$ in a system of $N=96$ neutrinos.

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 $\langle J^{A}_{z}(t)\rangle/(N/4)=1-2p(t)$, for a
system of $N=96$ neutrinos and different values for the asymmetry energy
$\delta_{\omega}$. Due to the conservation of $J_{z}$, the polarization of
beam B is simply the inverse. We can see strong flavor oscillations for
$0<\delta_{\omega}\lesssim\mu$ which for $\delta_{\omega}<\mu/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 $\delta_{\omega}/\mu$ (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 $\delta_{\omega}$ 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<\delta_{\omega}\lesssim\mu$ the survival probability
experiences instead fast oscillations on time scales much shorter than in the
high density limit $\delta_{\omega}=0$. If, in the latter case, the first
minimum of the survival probability is reached for a time $t_{P}$ that scales
with $\tau_{L}\propto\sqrt{N}$, in the former case we have instead evolution
at the fast scale $\tau_{F}\propto\log(N)$ found in Refs. Bell  _et al._
(2003); Sawyer (2004), but now for a physical model with $SU(2)$ invariance
(see also Roggero (2021)). 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 $\delta_{\omega}\neq 0$ a good fit to data is obtained with the
ansatz $\mu t_{P}(N)=a_{t}\log(N)+c_{t}$. The best fit values for the
parameters extracted from our simulations are reported in Tab. 1, together
with the extrapolated value of the minimum of the survival probability
$p_{min}$ in the $N\gg 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 $\delta_{\omega}$ outside the
critical region $(0.0,1.0)$, the latter converges to $p_{min}\approx 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}\propto\sqrt{\mu/\delta_{\omega}}$ similarly to the more familiar
bipolar oscillations Kostelecký and Samuel (1995); Hannestad  _et al._ (2006);
Duan  _et al._ (2007). 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 $\omega_{B}$ Duan  _et al._ (2010)
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 $\delta_{\omega}<0$ irrespective of
the energies, for normal hierarchy we have always positive energy differences
and oscillations can occur when $\omega_{A}+|\omega_{B}|\lesssim\mu$. 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

$0\leq\pm\frac{E_{B}-E_{A}}{E_{A}E_{B}}\lesssim\frac{4\sqrt{2}G_{F}\rho_{\nu}}{\Delta
m^{2}}\;,$ (33)

with $+$ for $\nu$ and $-$ for $\bar{\nu}$ in the normal hierarchy 333Note
that this condition can be modified by the presence of an MSW-like matter term
which couple differently to the two flavors..

$\delta_{\omega}/\mu$ | $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 1: Optimal parameters of the model discussed in the main text $\mu
t_{P}(N)=a_{t}\log(N)+b_{t}\sqrt{N}+c_{t}$ for time to reach the minimum of
the survival probability $t_{P}$. The minimum of the survival probability in
the thermodynamic limit is extracted using the ansatz
$p_{min}(N)=p_{min}+p_{c}/N^{\gamma}$ with exponents $\gamma=1$ for
$\delta_{\omega}=0$ and $\gamma=0.5$ otherwise.

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 $\delta_{\omega}=\mu/2$ (blue line) where it reaches a maximum
of $\approx 80\%$.

As an alternative way to quantify quantum correlations in the evolved state,
we compute the half-chain entanglement entropy (see eg. Eisert  _et al._
(2010)) defined explicitly as

$S_{N/2}(t)=-Tr\left[\rho_{B}(t)\log_{2}\left(\rho_{B}(t)\right)\right]\;.$
(34)

Here the density matrix $\rho_{B}=Tr_{A}\left[\rho(t)\right]$ is obtained by
tracing the full density matrix
$\rho(t)=\rvert\Psi(t)\rangle\langle\Psi(t)\lvert$ of the system at time $t$
over the first $N/2$ neutrinos.

Figure 4: (Color online) Half chain entanglement entropy $S_{N/2}$ for $N=96$
neutrinos and different values of the one-body energy difference: the orange
line is $\delta_{\omega}=\mu$, the blue line is for $\delta_{\omega}=\mu/2$,
the red line for $\delta_{\omega}=-\mu/2$ and the black line has
$\delta_{\omega}=0$. Panel (b) shows the maximum entropy $S_{max}$ for three
system sizes $N=16,32,64,128$ and different values of $\delta_{\omega}$. Panel
(c) shows the time $t_{ent}$ to reach the maximum.

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 $\delta_{\omega}$ 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
$\delta_{\omega}<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 $\delta_{\omega}\gtrsim\mu$ 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 $\delta_{\omega}=\mu$
(orange line), the oscillations in the entanglement entropy are much less
regular than in the $\delta_{\omega}<0$ phase and do not return close to zero.
At the transition point $\delta_{\omega}=0$ we find that, after an initial
rapid growth phase, the entropy $S_{N/2}(t)$ reaches a peak at
$S_{max}\approx\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 Vidal  _et al._ (2003); Refael and Moore (2004) 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 Schuch  _et al._ (2008). We find a
similar behavior for the maximum entropy reached in the dynamics for the
intermediate regime $0<\delta_{\omega}/\mu\lesssim 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 Roggero
(2021). 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 Heyl (2018); Haldar  _et al._ (2020).

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 $\delta_{\omega}<0$ is approximately independent from system
size, at the critical point we recover the familiar $t_{P}\approx\sqrt{N}$ and
for positive $\delta_{\omega}>0$ we find a logarithmic evolution
$t_{P}\approx\log(N)$ which seems to approach a constant again for large value
of $\delta_{\omega}$. In this regime the system seems to enter a second
”gapped” phase, similar to the one for $\delta_{\omega}<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 Vidal (2003); Schuch  _et al._ (2008).

##  IV Summary and Perspective

Entanglement is a fundamental feature of interacting quantum many body systems
Vidal  _et al._ (2003); Eisert  _et al._ (2010) which, somewhat surprisingly,
can be shown to play no important role in determining ground-state properties
of systems with all-to-all interactions like for the neutrino forward-
scattering Hamiltonian of Eq. (1) (see eg. Brandão and Harrow (2016)). 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 Cervia  _et al._ (2019); Rrapaj (2020).

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 Hannestad  _et al._ (2006); Duan  _et
al._ (2010). In order to overcome the exponential complexity of a direct
solution of the problem we use a simple technique, based on a MPS
representation Vidal (2003), 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 Vidal (2008) which naturally describe states with
logarithmically increasing entropy.

As shown in an accompanying paper Roggero (2021), 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 Roggero (2021). 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 Zhang  _et al._ (2017) or future circuit based quantum
computers Hall  _et al._ (2021), 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.

###### Acknowledgements.

I want to thank Joseph Carlson, Vincenzo Cirigliano, Huaiyu Duan, Joshua
Martin, Amol Patwardhan, Ermal Rrapaj and Martin Savage for the many useful
discussions about the subject of this work. This work was supported by the
InQubator for Quantum Simulation under U.S. DOE grant No. DE-SC0020970 and by
the Quantum Science Center (QSC), a National Quantum Information Science
Research Center of the U.S. Department of Energy (DOE).

## References

  * Hoffman  _et al._ (1997) R. D. Hoffman, S. E. Woosley, and Y.-Z. Qian, “Nucleosynthesis in neutrino-driven winds. II. implications for heavy element synthesis,” The Astrophysical Journal 482, 951–962 (1997).
  * Janka (2012) Hans-Thomas Janka, “Explosion mechanisms of core-collapse supernovae,” Annual Review of Nuclear and Particle Science 62, 407–451 (2012).
  * Winteler  _et al._ (2012) C. Winteler, R. Käppeli, A. Perego, A. Arcones, N. Vasset, N. Nishimura, M. Liebendörfer, and F.-K. Thielemann, “MAGNETOROTATIONALLY DRIVEN SUPERNOVAE AS THE ORIGIN OF EARLY GALAXY r-PROCESS ELEMENTS?” The Astrophysical Journal 750, L22 (2012).
  * Wanajo  _et al._ (2014) Shinya Wanajo, Yuichiro Sekiguchi, Nobuya Nishimura, Kenta Kiuchi, Koutarou Kyutoku, and Masaru Shibata, “PRODUCTION OF ALL THE r-PROCESS NUCLIDES IN THE DYNAMICAL EJECTA OF NEUTRON STAR MERGERS,” The Astrophysical Journal 789, L39 (2014).
  * Savage  _et al._ (1991) Martin J. Savage, Robert A. Malaney, and George M. Fuller, “Neutrino Oscillations and the Leptonic Charge of the Universe,” Astrophys. J. 368, 1 (1991).
  * Pantaleone (1992a) James Pantaleone, “Dirac neutrinos in dense matter,” Phys. Rev. D 46, 510–523 (1992a).
  * Pantaleone (1992b) James Pantaleone, “Neutrino oscillations at high densities,” Physics Letters B 287, 128 – 132 (1992b).
  * Pastor and Raffelt (2002) Sergio Pastor and Georg Raffelt, “Flavor oscillations in the supernova hot bubble region: Nonlinear effects of neutrino background,” Phys. Rev. Lett. 89, 191101 (2002).
  * Balantekin and Yüksel (2005) A B Balantekin and H Yüksel, “Neutrino mixing and nucleosynthesis in core-collapse supernovae,” New Journal of Physics 7, 51–51 (2005).
  * Fuller and Qian (2006) George M. Fuller and Yong-Zhong Qian, “Simultaneous flavor transformation of neutrinos and antineutrinos with dominant potentials from neutrino-neutrino forward scattering,” Phys. Rev. D 73, 023004 (2006).
  * Duan  _et al._ (2006a) Huaiyu Duan, George M. Fuller, J. Carlson, and Yong-Zhong Qian, “Coherent development of neutrino flavor in the supernova environment,” Phys. Rev. Lett. 97, 241101 (2006a).
  * Friedland (2010) Alexander Friedland, “Self-refraction of supernova neutrinos: Mixed spectra and three-flavor instabilities,” Phys. Rev. Lett. 104, 191102 (2010).
  * Martin  _et al._ (2020) Joshua D. Martin, Changhao Yi, and Huaiyu Duan, “Dynamic fast flavor oscillation waves in dense neutrino gases,” Physics Letters B 800, 135088 (2020).
  * Wu and Tamborra (2017) Meng-Ru Wu and Irene Tamborra, “Fast neutrino conversions: Ubiquitous in compact binary merger remnants,” Phys. Rev. D 95, 103007 (2017).
  * Fogli  _et al._ (2007) Gianluigi Fogli, Eligio Lisi, Antonio Marrone, and Alessandro Mirizzi, “Collective neutrino flavor transitions in supernovae and the role of trajectory averaging,” Journal of Cosmology and Astroparticle Physics 2007, 010–010 (2007).
  * Qian  _et al._ (1993) 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 significant neutrinos and heavy element nucleosynthesis in supernovae,” Phys. Rev. Lett. 71, 1965–1968 (1993).
  * Qian and Fuller (1995) Yong-Zhong Qian and George M. Fuller, “Neutrino-neutrino scattering and matter-enhanced neutrino flavor transformation in supernovae,” Phys. Rev. D 51, 1479–1494 (1995).
  * Duan and Friedland (2011) Huaiyu Duan and Alexander Friedland, “Self-induced suppression of collective neutrino oscillations in a supernova,” Phys. Rev. Lett. 106, 091101 (2011).
  * Dasgupta  _et al._ (2012) Basudeb Dasgupta, Evan P. O’Connor, and Christian D. Ott, “Role of collective neutrino flavor oscillations in core-collapse supernova shock revival,” Phys. Rev. D 85, 065008 (2012).
  * Sawyer (2005) R. F. Sawyer, “Speed-up of neutrino transformations in a supernova environment,” Phys. Rev. D 72, 045003 (2005).
  * Sawyer (2016) R. F. Sawyer, “Neutrino cloud instabilities just above the neutrino sphere of a supernova,” Phys. Rev. Lett. 116, 081101 (2016).
  * Chakraborty  _et al._ (2016a) S. Chakraborty, R. S. Hansen, I. Izaguirre, and G.G. Raffelt, “Self-induced neutrino flavor conversion without flavor mixing,” Journal of Cosmology and Astroparticle Physics 2016, 042–042 (2016a).
  * Izaguirre  _et al._ (2017) Ignacio Izaguirre, Georg Raffelt, and Irene Tamborra, “Fast pairwise conversion of supernova neutrinos: A dispersion relation approach,” Phys. Rev. Lett. 118, 021101 (2017).
  * Chakraborty  _et al._ (2016b) Sovan Chakraborty, Rasmus Hansen, Ignacio Izaguirre, and Georg Raffelt, “Collective neutrino flavor conversion: Recent developments,” Nuclear Physics B 908, 366 – 381 (2016b), neutrino Oscillations: Celebrating the Nobel Prize in Physics 2015.
  * Vlasenko  _et al._ (2014) Alexey Vlasenko, George M. Fuller, and Vincenzo Cirigliano, “Neutrino quantum kinetics,” Phys. Rev. D 89, 105004 (2014).
  * Cirigliano  _et al._ (2015) Vincenzo Cirigliano, George M. Fuller, and Alexey Vlasenko, “A new spin on neutrino quantum kinetics,” Physics Letters B 747, 27–35 (2015).
  * Pehlivan  _et al._ (2011) Y. Pehlivan, A. B. Balantekin, Toshitaka Kajino, and Takashi Yoshida, “Invariants of collective neutrino oscillations,” Phys. Rev. D 84, 065008 (2011).
  * Note (1) Note the difference in the definition of $\mu$ there.
  * Stirner  _et al._ (2018) Tobias Stirner, Günter Sigl, and Georg Raffelt, “Liouville term for neutrinos: flavor structure and wave interpretation,” Journal of Cosmology and Astroparticle Physics 2018, 016–016 (2018).
  * Sawyer (2004) R.F. Sawyer, “’Classical’ instabilities and ’quantum’ speed-up in the evolution of neutrino clouds,” (2004), arXiv:hep-ph/0408265 .
  * Duan  _et al._ (2010) Huaiyu Duan, George M. Fuller, and Yong-Zhong Qian, “Collective neutrino oscillations,” Annual Review of Nuclear and Particle Science 60, 569–594 (2010), https://doi.org/10.1146/annurev.nucl.012809.104524 .
  * Eisert  _et al._ (2015) J. Eisert, M. Friesdorf, and C. Gogolin, “Quantum many-body systems out of equilibrium,” Nature Physics 11, 124–130 (2015).
  * Cervia  _et al._ (2019) Michael J. Cervia, Amol V. Patwardhan, A. B. Balantekin, S. N. Coppersmith, and Calvin W. Johnson, “Entanglement and collective flavor oscillations in a dense neutrino gas,” Phys. Rev. D 100, 083001 (2019).
  * Rrapaj (2020) Ermal Rrapaj, “Exact solution of multiangle quantum many-body collective neutrino-flavor oscillations,” Phys. Rev. C 101, 065805 (2020).
  * Friedland and Lunardini (2003a) Alexander Friedland and Cecilia Lunardini, “Do many-particle neutrino interactions cause a novel coherent effect?” Journal of High Energy Physics 2003, 043–043 (2003a).
  * Friedland  _et al._ (2006) Alexander Friedland, Bruce H. J. McKellar, and Ivona Okuniewicz, “Construction and analysis of a simplified many-body neutrino model,” Phys. Rev. D 73, 093002 (2006).
  * Bell  _et al._ (2003) Nicole F. Bell, Andrew A. Rawlinson, and R.F. Sawyer, “Speed-up through entanglement—many-body effects in neutrino processes,” Physics Letters B 573, 86 – 93 (2003).
  * Note (2) Since $[H,J^{2}_{A}]=[H,J^{2}_{B}]=0$ we are free to change the intra-beam couplings at will thus changing the geometry. In the general case we can substitute $\mu\to\mu(1-\mathop{cos}\nolimits(\theta_{AB}))$, with $\theta_{AB}$ the angle between the momenta in the $A$ and $B$ beams, for all the results in this work.
  * Lipkin  _et al._ (1965) H.J. Lipkin, N. Meshkov, and A.J. Glick, “Validity of many-body approximation methods for a solvable model: (i). exact solutions and perturbation theory,” Nuclear Physics 62, 188 – 198 (1965).
  * Vidal  _et al._ (2004a) Julien Vidal, Rémy Mosseri, and Jorge Dukelsky, “Entanglement in a first-order quantum phase transition,” Phys. Rev. A 69, 054101 (2004a).
  * Vidal  _et al._ (2004b) Julien Vidal, Guillaume Palacios, and Claude Aslangul, “Entanglement dynamics in the lipkin-meshkov-glick model,” Phys. Rev. A 70, 062304 (2004b).
  * Vidal  _et al._ (2004c) Julien Vidal, Guillaume Palacios, and Rémy Mosseri, “Entanglement in a second-order quantum phase transition,” Phys. Rev. A 69, 022107 (2004c).
  * Latorre  _et al._ (2005) José I. Latorre, Román Orús, Enrique Rico, and Julien Vidal, “Entanglement entropy in the lipkin-meshkov-glick model,” Phys. Rev. A 71, 064101 (2005).
  * Ribeiro  _et al._ (2008) Pedro Ribeiro, Julien Vidal, and Rémy Mosseri, “Exact spectrum of the lipkin-meshkov-glick model in the thermodynamic limit and finite-size corrections,” Phys. Rev. E 78, 021106 (2008).
  * Friedland and Lunardini (2003b) Alexander Friedland and Cecilia Lunardini, “Neutrino flavor conversion in a neutrino background: Single- versus multi-particle description,” Phys. Rev. D 68, 013007 (2003b).
  * Vidal (2003) Guifré Vidal, “Efficient classical simulation of slightly entangled quantum computations,” Phys. Rev. Lett. 91, 147902 (2003).
  * Vidal (2004) Guifré Vidal, “Efficient simulation of one-dimensional quantum many-body systems,” Phys. Rev. Lett. 93, 040502 (2004).
  * Paeckel  _et al._ (2019) Sebastian Paeckel, Thomas Köhler, Andreas Swoboda, Salvatore R. Manmana, Ulrich Schollwöck, and Claudius Hubig, “Time-evolution methods for matrix-product states,” Annals of Physics 411, 167998 (2019).
  * Schuch  _et al._ (2008) Norbert Schuch, Michael M. Wolf, Frank Verstraete, and J. Ignacio Cirac, “Entropy scaling and simulability by matrix product states,” Phys. Rev. Lett. 100, 030504 (2008).
  * Schollwöck (2011) Ulrich Schollwöck, “The density-matrix renormalization group in the age of matrix product states,” Annals of Physics 326, 96 – 192 (2011), january 2011 Special Issue.
  * Hall  _et al._ (2021) Benjamin Hall, Alessandro Roggero, Alessandro Baroni, and Joseph Carlson, “Simulation of Collective Neutrino Oscillations on a Quantum Computer,” arXiv e-prints , arXiv:2102.12556 (2021), arXiv:2102.12556 [quant-ph] .
  * Verstraete and Cirac (2006) F. Verstraete and J. I. Cirac, “Matrix product states represent ground states faithfully,” Phys. Rev. B 73, 094423 (2006).
  * Suzuki (1991) Masuo Suzuki, “General theory of fractal path integrals with applications to many‐body theories and statistical physics,” Journal of Mathematical Physics 32, 400–407 (1991).
  * Stoudenmire and White (2010) E M Stoudenmire and Steven R White, “Minimally entangled typical thermal state algorithms,” New Journal of Physics 12, 055026 (2010).
  * Kivlichan  _et al._ (2018) Ian D. Kivlichan, Jarrod McClean, Nathan Wiebe, Craig Gidney, Alán Aspuru-Guzik, Garnet Kin-Lic Chan, and Ryan Babbush, “Quantum simulation of electronic structure with linear depth and connectivity,” Phys. Rev. Lett. 120, 110501 (2018).
  * Fishman  _et al._ (2020) Matthew Fishman, Steven R. White, and E. Miles Stoudenmire, “The ITensor software library for tensor network calculations,” (2020), arXiv:2007.14822 .
  * Roggero (2021) Alessandro Roggero, “Dynamical Phase Transitions in models of Collective Neutrino Oscillations,” arXiv e-prints , arXiv:2103.11497 (2021), arXiv:2103.11497 [hep-ph] .
  * Kostelecký and Samuel (1995) V. Alan Kostelecký and Stuart Samuel, “Self-maintained coherent oscillations in dense neutrino gases,” Phys. Rev. D 52, 621–627 (1995).
  * Hannestad  _et al._ (2006) Steen Hannestad, Georg G. Raffelt, Günter Sigl, and Yvonne Y. Y. Wong, “Self-induced conversion in dense neutrino gases: Pendulum in flavor space,” Phys. Rev. D 74, 105010 (2006).
  * Duan  _et al._ (2007) Huaiyu Duan, George M. Fuller, J. Carlson, and Yong-Zhong Qian, “Analysis of collective neutrino flavor transformation in supernovae,” Phys. Rev. D 75, 125005 (2007).
  * Note (3) Note that this condition can be modified by the presence of an MSW-like matter term which couple differently to the two flavors.
  * Eisert  _et al._ (2010) J. Eisert, M. Cramer, and M. B. Plenio, “Colloquium: Area laws for the entanglement entropy,” Rev. Mod. Phys. 82, 277–306 (2010).
  * Vidal  _et al._ (2003) G. Vidal, J. I. Latorre, E. Rico, and A. Kitaev, “Entanglement in quantum critical phenomena,” Phys. Rev. Lett. 90, 227902 (2003).
  * Refael and Moore (2004) G. Refael and J. E. Moore, “Entanglement entropy of random quantum critical points in one dimension,” Phys. Rev. Lett. 93, 260602 (2004).
  * Heyl (2018) Markus Heyl, “Dynamical quantum phase transitions: a review,” Reports on Progress in Physics 81, 054001 (2018).
  * Haldar  _et al._ (2020) Asmi Haldar, Krishnanand Mallayya, Markus Heyl, Frank Pollmann, Marcos Rigol, and Arnab Das, “Signatures of quantum phase transitions after quenches in quantum chaotic one-dimensional systems,” arXiv e-prints , arXiv:2004.02905 (2020), arXiv:2004.02905 [cond-mat.stat-mech] .
  * Brandão and Harrow (2016) Fernando G. S. L. Brandão and Aram W. Harrow, “Product-state approximations to quantum states,” Communications in Mathematical Physics 342, 47–80 (2016).
  * Vidal (2008) G. Vidal, “Class of quantum many-body states that can be efficiently simulated,” Phys. Rev. Lett. 101, 110501 (2008).
  * Zhang  _et al._ (2017) 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 dynamical phase transition with a 53-qubit quantum simulator,” Nature 551, 601–604 (2017).
  * Childs  _et al._ (2021) Andrew M. Childs, Yuan Su, Minh C. Tran, Nathan Wiebe, and Shuchen Zhu, “Theory of trotter error with commutator scaling,” Phys. Rev. X 11, 011020 (2021).
  * Duan  _et al._ (2006b) Huaiyu Duan, George M. Fuller, and Yong-Zhong Qian, “Collective neutrino flavor transformation in supernovae,” Phys. Rev. D 74, 123004 (2006b).
  * Haegeman  _et al._ (2011) Jutho Haegeman, J. Ignacio Cirac, Tobias J. Osborne, Iztok Pižorn, Henri Verschelde, and Frank Verstraete, “Time-dependent variational principle for quantum lattices,” Phys. Rev. Lett. 107, 070601 (2011).

##  Appendix A Analysis of time-step errors

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 approximate evolution operator used in our simulations

$U(\delta)=\prod_{i<j=1}^{N}e^{-i\delta h_{ij}}\prod^{1}_{i<j=N}e^{-i\delta
h_{ij}}\;,$ (35)

and the very conservative error bound (see e.g. Childs  _et al._ (2021) for
tighter ones) presented in Eq. (26) of the main text

$\left\|e^{-i\delta H}-U(\delta)\right\|<\delta^{3}\|H\|^{3}\;,$ (36)

where $\|\cdot\|$ 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 $\delta\ll t$ as described in Eq. (32) of the main text. Using the union
bound, one can straightforwardly generalize Eq. (37) to arbitrary long times
as follows

$\left\|e^{-itH}-U(\delta)^{\frac{t}{\delta}}\right\|<\frac{t}{\delta}\delta^{3}\|H\|^{3}=t\delta^{2}\|H\|^{3}\;.$
(37)

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
$\epsilon$ it is sufficient to choose

$\delta\leq\sqrt{\frac{\epsilon}{t\|H\|^{3}}}\quad\text{or}\quad\frac{t}{\delta}\geq\frac{(t\|H\|)^{\frac{3}{2}}}{\sqrt{\epsilon}}\;,$
(38)

where $t/\delta$ is the total number of steps in the simulation.

Figure 5: (Color online) Top panel shows the half-chain entanglement entropy
in a system of $N=96$ neutrino amplitudes for two values of the one body
asymmetry $\delta_{\omega}=0$ and $\delta_{\omega}=\mu/2$ and two different
time steps $\delta=0.05\mu^{-1}$ (continuous lines) and $\delta=0.2\mu^{-1}$
(dots). The bottom panel shows the relative percent error caused by time step
errors.

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 $\delta_{\omega}=0$ and $\delta_{\omega}=\mu/2$ using a time
step $\delta=0.2\mu^{-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 $\delta_{\omega}$ 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 $\delta=0.05\mu^{-1}$ (the choice used for
the results in the main text) and the dots indicate the results with
$\delta=0.2\mu^{-1}$ which are expected to have $\approx 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 $\approx 1\%$ for the entanglement
entropy and $\approx 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.

Figure 6: (Color online) Top panel shows the survival probability $p(t)$ in a
system of $N=96$ neutrino amplitudes for two values of the one body asymmetry
$\delta_{\omega}=0$ and $\delta_{\omega}=\mu/2$ and two different time steps
$\delta=0.05\mu^{-1}$ (continuous lines) and $\delta=0.2\mu^{-1}$ (dots). The
bottom panel shows the relative percent error caused by time step errors.

##  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\to\infty$ 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

$H_{f}=\delta_{\omega}\vec{B}\cdot\left(\vec{J}^{A}-\vec{J}^{B}\right)+\frac{\mu}{N}\left(\vec{J}^{A}+\vec{J}^{B}\right)^{2}$
(39)

with $\vec{J}^{A}$,$\vec{J}^{B}$ the total spin operators in the two beams and
$\vec{B}=(\sin(2\theta),0,-\cos(2\theta))$ is related to the mixing angle
$\theta$. In the following, we will also keep the same initial condition
$\ket{\Psi_{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

$\begin{split}\frac{d}{dt}\vec{J}_{A}&=\left[\delta_{\omega}\vec{B}+\frac{2\mu}{N}\left(\vec{J}_{A}+\vec{J}_{B}\right)\right]\times\vec{J}_{A}\\\
\frac{d}{dt}\vec{J}_{B}&=\left[-\delta_{\omega}\vec{B}+\frac{2\mu}{N}\left(\vec{J}_{A}+\vec{J}_{B}\right)\right]\times\vec{J}_{B}\;.\end{split}$
(40)

One can eliminate the explicit dependence on the number of amplitudes $N$ by
expressing these equations in terms of polarization vectors
$\vec{P}_{A}=\frac{4}{N}\langle\vec{J}_{A}\rangle$ and
$\vec{P}_{B}=\frac{4}{N}\langle\vec{J}_{B}\rangle$. In the mean field
approximation we have then

$\begin{split}\frac{d}{dt}\vec{P}_{A}&=\left[\delta_{\omega}\vec{B}+\frac{\mu}{2}\left(\vec{P}_{A}+\vec{P}_{B}\right)\right]\times\vec{P}_{A}\\\
\frac{d}{dt}\vec{P}_{B}&=\left[-\delta_{\omega}\vec{B}+\frac{\mu}{2}\left(\vec{P}_{A}+\vec{P}_{B}\right)\right]\times\vec{P}_{B}\;.\end{split}$
(41)

The solution of these coupled equations can be simplified with the change of
variables proposed in Refs. Duan  _et al._ (2006b); Hannestad  _et al._
(2006).

In Fig. 7 we compare the result for the time evolution of the survival
probability

$P(t)=\frac{1}{2}\left(1-\vec{P}_{A}\right)=\frac{1}{2}\left(1+\vec{P}_{B}\right)\;,$
(42)

as obtained by either solving the equations in Eq. (41) (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. Duan  _et al._ (2006b), upon rescaling of the time variable). The
different colored curves correspond to different choices for the time-step
$\delta$ used to simulate time-evolution with Eq. (32). The results show a
discrepancy with the mean-field 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.2 is one of the main motivation behind alternatives like the Time Dependent
Variational Principle(TDVP) approach Haegeman  _et al._ (2011) which, instead,
approximates globally optimal tensors remaining directly on the manifold of
matrix product states with fixed bond dimension. As the operators
$G^{\sigma_{1}\sigma_{2}}_{\sigma_{1}^{\prime}\sigma_{2}^{\prime}}$ 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 $\delta$ (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.

Figure 7: (Color online) Survival probability $P(t)$ as a function of time for
a system with $N=24$, $\delta_{\omega}/\mu=0.1$ and bond dimension fixed at
$r=1$ for different values of the time-step $\delta$. The solid black curve
corresponds to the solution of the mean-field equations Eq. (41).