Entanglement and collective flavor oscillations in a dense neutrino gas

TL;DR

This paper explores how quantum entanglement influences collective neutrino oscillations in dense gases, highlighting the limitations of mean-field approximations and employing many-body physics tools.

Contribution

It introduces the use of entanglement measures to analyze neutrino oscillations, advancing understanding beyond traditional mean-field models.

Findings

Entanglement significantly affects neutrino oscillation dynamics.

Mean-field approximation can overlook important quantum correlations.

Comparisons show differences between many-body and mean-field evolution.

Abstract

We investigate the importance of going beyond the mean-field approximation in the dynamics of collective neutrino oscillations. To expand our understanding of the coherent neutrino oscillation problem, we apply concepts from many-body physics and quantum information theory. Specifically, we use measures of nontrivial correlations (otherwise known as "entanglement") between the constituent neutrinos of the many-body system, such as the entanglement entropy and the Bloch vector of the reduced density matrix. The relevance of going beyond the mean field is demonstrated by comparisons between the evolution of the neutrino state in the many-body picture vs the mean-field limit, for different initial conditions.

Tables

Table 1. Table 1: Summary of parameters used in our calculations.

	Value
Parameter	Number	Unit
$R_{\nu }$	$60$	km
${\omega }_0$	$1.055\times {10}^{-16}$	MeV
$\mu \left(R_{\nu }\right)$	$3.62\times {10}^4$	${\omega }_0$
${\sin }^2\left(2\theta \right)$	$0.10$

Equations

$$H=\sum _{\mathbf{p}}{\omega }_{\mathbf{p}}\vec{B}\cdot {\vec{J}}_{\mathbf{p}}+\sum _{\mathbf{p},\mathbf{q}}{\mu }_{\mathbf{pq}}{\vec{J}}_{\mathbf{p}}\cdot {\vec{J}}_{\mathbf{q}},$$

$${\mu }_{\mathbf{pq}}=\frac{\sqrt{2}{G}_F}{V}(1-\cos {\theta }_{\mathbf{pq}}),$$

$$\mu(r)=\frac{G_{F}}{\sqrt{2}V}\bigg{[}1-\bigg{(}1-\frac{R_{\nu}^{2}}{r^{2}}\bigg{)}^{1/2}\bigg{]}^{2},$$

$$H=-\sum _{p=1}^M{\omega }_p{J}_p^z+\mu \vec{J}\cdot \vec{J},$$

$$h_p=-J_p^z+2\mu \sum _{\begin{array}{c}q=1\\ q\ne p\end{array}}^M\frac{{\vec{J}}_p\cdot {\vec{J}}_q}{{\omega }_p-{\omega }_q}.$$

$$|q⟩=\sum _{i_1,\dots ,i_N=1}^2{V}_{i_1,\dots ,i_N}(\mu )|{\nu }_{i_1},\dots ,{\nu }_{i_N}⟩$$

$$i\text{diff}|\Psi ⟩t=H|\Psi ⟩.$$

$$i\text{diff}|\Phi ⟩t+i{V}^{†}\text{diff}Vt|\Phi ⟩=\Sigma |\Phi ⟩.$$

$$|\Phi (t)⟩=e^{-i{\int }^t\Sigma d{t}^{\prime }}|{\Phi }_0⟩.$$

$$|\Psi (t)⟩=V{e}^{-i{\int }^t\Sigma d{t}^{\prime }}{V}_0^{†}|{\Psi }_0⟩,$$

$$|\Psi (\mu )⟩=$$

$$\times V_0^T|\psi ({\mu }_0)⟩$$

$${\rho }_A$$

$$S_{(A|B)}=-\text{Tr}[{\rho }_A\log {\rho }_A]=-\text{Tr}[{\rho }_B\log {\rho }_B],$$

$${\rho }_A=\sum _{i_1,\dots ,i_{N-1}=1,2}⟨{\nu }_{i_1}\dots {\nu }_{i_{N-1}}⟩\rho {\nu }_{i_1}\dots {\nu }_{i_{N-1}}$$

$${\rho }_A=\frac{1}{2}(\mathbf{I}+\vec{\sigma }\cdot {\vec{P}}_A),$$

$$S=-\frac{1-P}{2}\log \left(\frac{1-P}{2}\right)-\frac{1+P}{2}\log \left(\frac{1+P}{2}\right).$$

$$\Delta P_z({\omega }_p)\equiv |P_z^{\mathrm{MF}}({\omega }_p)-P_z^{\mathrm{MB}}({\omega }_p)|$$

$$\Delta P_z({\omega }_p)=1-P(S({\omega }_p)),$$

$$H=\left(\begin{array}{cccc}-\Omega +2\mu & 0& 0& 0\\ 0& -\eta +\mu & \mu & 0\\ 0& \mu & \eta +\mu & 0\\ 0& 0& 0& \Omega +2\mu \end{array}\right)$$

$$\Sigma =\left(\begin{array}{cccc}-\Omega +\mu & 0& 0& 0\\ 0& \sqrt{{\eta }^2+{\mu }^2}& 0& 0\\ 0& 0& -\sqrt{{\eta }^2+{\mu }^2}& 0\\ 0& 0& 0& \Omega +\mu \end{array}\right),$$

$$V=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& \cos \xi & -\sin \xi & 0\\ 0& \sin \xi & \cos \xi & 0\\ 0& 0& 0& 1\end{array}\right),$$

$$\cos 2\xi =\frac{-\eta }{\sqrt{{\mu }^2+{\eta }^2}}\mathrm{and}\sin 2\xi =\frac{\mu }{\sqrt{{\mu }^2+{\eta }^2}}.$$

$$\alpha (t)=$$

$$\beta (t)=$$

$$-i\sin Q[{\beta }_0\cos (\xi +{\xi }_0)+{\gamma }_0\sin (\xi +{\xi }_0)]$$

$$\gamma (t)=$$

$$-i\sin Q[{\beta }_0\sin (\xi +{\xi }_0)-{\gamma }_0\cos (\xi +{\xi }_0)]$$

$$\delta (t)=$$

$${\rho }_2^{\text{red}}=\left(\begin{array}{cc}|\alpha {|}^2+|\gamma {|}^2& \alpha {\beta }^{\ast }+\gamma {\delta }^{\ast }\\ {\alpha }^{\ast }\beta +{\gamma }^{\ast }\delta & |\beta {|}^2+|\delta {|}^2\end{array}\right).$$

$$\begin{array}{rl}{P}_z=& |\alpha {|}^2-|\beta {|}^2+|\gamma {|}^2-|\delta {|}^2\\ =& {\alpha }_0^2-{\delta }_0^2-({\beta }_0^2-{\gamma }_0^2)\\ & \times [\cos 2\xi \cos 2{\xi }_0+\sin 2\xi \sin 2{\xi }_0\cos 2Q]\\ & -2{\beta }_0{\gamma }_0[\cos 2\xi \sin 2{\xi }_0-\sin 2\xi \cos 2{\xi }_0\cos 2Q].\end{array}$$

$$P_{\perp }=2|{\alpha }^{\ast }\beta +{\gamma }^{\ast }\delta |.$$

$$P^2=$$

Full text

Suggested keywords

††preprint: AAPM/123-QED

I Introduction

Collective neutrino oscillations, which result from coherent forward scattering of neutrinos off other free-streaming neutrinos, is a topic of long-standing interest, particularly in connection with neutrinos coming from core-collapse supernovae and merging binary neutron stars, as well as neutrinos in the early universe Fuller et al. (1987); Nötzold and Raffelt (1988); Pantaleone (1992); Sigl and Raffelt (1993); Raffelt et al. (1993); Samuel (1993); Kostelecky and Samuel (1994); Kostelecký and Samuel (1995); Dolgov et al. (2002); Ho et al. (2005); McKellar and Thomson (1994); Enqvist et al. (1991); Pastor et al. (2002); Abazajian et al. (2002); Hannestad et al. (2006); Johns et al. (2016); Duan and Kneller (2009); Duan et al. (2010); Chakraborty et al. (2016); Balantekin (2018); Duan et al. (2006a, b, c, 2007a, 2007b, 2007c, 2008); Raffelt and Smirnov (2007, 2007); Fogli et al. (2007); Dasgupta and Dighe (2008); Dasgupta et al. (2008, 2009); Friedland (2010); Dasgupta et al. (2010); Galais and Volpe (2011); Mirizzi (2013); Chakraborty and Mirizzi (2014); Pehlivan et al. (2017); Tian et al. (2017); Das et al. (2017); Birol et al. (2018); Cirigliano et al. (2018); Malkus et al. (2012, 2014); Wu et al. (2016); Stapleford et al. (2016); Malkus et al. (2016); Frensel et al. (2017); Väänänen and McLaughlin (2016); Zhu et al. (2016); Johns et al. (2016); Shalgar (2018); Vlasenko and McLaughlin (2018); Sawyer (2005); Chakraborty et al. (2016); Capozzi et al. (2016); Wu and Tamborra (2017); Sen (2018); Dighe and Sen (2018); Abbar and Duan (2018); Wu et al. (2017); Dasgupta and Sen (2018); Dasgupta et al. (2018); Capozzi et al. (2019a); Abbar and Volpe (2019); Shalgar and Tamborra ; Capozzi et al. (2019b) (see, e.g., reviews in Refs. Duan et al. (2010); Duan and Kneller (2009); Chakraborty et al. (2016)). Terrestrial experiments such as DUNE, with the capability to detect supernova neutrinos, could potentially be used to discern the signatures of collective oscillations in the observed energy spectra and flavor content of the neutrinos and thereby yield new insight into the nature of neutrinos and the physics of core-collapse supernovae.

It is not practical to determine the exact evolution of $10^{58}$ or so neutrinos and antineutrinos emitted during the core collapse of a typical supernova. For example, for a two-flavor system of $N$ neutrinos (with no antineutrinos), the dynamical variable, i.e., the $N$-particle wave function, occupies a $2^{N}$-dimensional Hilbert space that is the direct product of the individual one-body Hilbert spaces. Therefore it is common to adopt simplifying approaches such as the mean-field approximation. This approach, however, is rather drastic; in the mean-field approximation, all quantum correlations between particles are ignored, ensuring that the dynamical variable is always separable into $N$ single-particle wave functions, thereby effectively reducing the size of the Hilbert space to just $2N$. It is therefore crucial to test the limits of validity of the mean-field approximation in different scenarios.

The size of the Hilbert space in an interacting many-body system scales exponentially with the number of particles; it is therefore computationally expensive to analyze differences between many-body and mean-field approaches in realistic supernova models with large numbers of neutrinos. In this work, we therefore adopt a toy model consisting of a few neutrinos, each having discrete well-defined momenta. This discretization necessitates considering the neutrinos as plane waves, and as a result the neutrino-neutrino interactions in our model are not localized in space.

Some previous treatments have used this approach, albeit with additional simplifications such as ignoring vacuum oscillations. For instance, it has been argued that, when the neutrinos start from an initial configuration that is asymmetrically distributed in momentum space, collisions can cause them to become uniformly distributed faster than one would expect from single-particle cross sections, due to many-body entanglement Bell et al. (2003). On the other hand, it was shown that, for certain geometries (i.e., orthogonal neutrino beams) in the absence of vacuum oscillations, the contribution from entangled states becomes insignificant if successive neutrino-neutrino scatterings are treated as independent events Friedland and Lunardini (2003a, b). More recently, there have been attempts to characterize the neutrino many-body Hamiltonian in terms of its eigenvalues and eigenstates Balantekin and Pehlivan (2007); Pehlivan et al. (2011); Birol et al. (2018); Patwardhan et al. (2019), and this insight has been applied to study the evolution of such a system for a particular initial condition Birol et al. (2018).

As an aside, previous authors have also explored entanglement of neutrinos with charged leptons in particle decays Kayser et al. (2010); Boyanovsky (2011); Floerchinger and Schwindt and in electron capture Pavlichenkov (2011), flavor entanglement during both free propagation and interaction with the background environment Blasone et al. (2009); Akhmedov and Smirnov (2011); Wu et al. (2011); Blasone et al. (2015); Evslin et al. (2019), as well as entanglement between a two-level system (which could model two neutrino flavors) and a thermal bath (which could represent the environment) Bell et al. (2002). However, in these cases, neutrino propagation is treated as a one-body problem from the perspective of the neutrino, in contrast with collective neutrino oscillations which arise from neutrino-neutrino interactions and represent a many-body problem.

In this work we focus on quantifying the entanglement that can develop in a many-body neutrino system as described by our toy model. For this purpose, we calculate the time evolution of the entanglement entropy of the different neutrinos in systems with varying sizes and initial configurations. Since entanglement is absent in mean-field treatments, an entanglement measure can serve as a quantifier for the extent to which many-body systems can deviate from the mean-field approximation. Indeed, we observe that the entropy of entanglement appears to be correlated with the magnitude of the differences in the flavor evolution between many-body and mean-field results, with significant deviations from the mean field observed in some cases.

In our analysis we consider adiabatic evolution of the many-neutrino system, with at most one neutrino in each energy bin, in the single-angle approximation. The equations defining the model that we study are presented in Sec. II. Adiabatic eigenvalues and eigenstates of this many-neutrino Hamiltonian are obtained Pehlivan et al. (2011); Birol et al. (2018); Patwardhan et al. (2019) using the Richardson-Gaudin algebraic approach Richardson (1966); Gaudin (1976); Balantekin (2018). This solution is outlined briefly in Sec. III. In Sec. IV we define the entanglement measures that we have used in our analysis. Section V presents an overview of the main results, in particular the relationship between entanglement entropy and the deviation of the many-body results from their mean-field counterparts, for small-$N$ systems up to $N=9$. Section VI presents the solution for a two-neutrino system, which can be obtained analytically. A short summary of the mean-field evolution equations for two neutrinos is also given. Section VII contains detailed numerical solutions for the cases where the number of neutrinos is greater than two. Finally, Sec. VIII contains a brief discussion of our conclusions. In Appendix A, we summarize our procedure for obtaining eigenvalues and eigenvectors of the single-angle neutrino Hamiltonian, and in Appendix B we explore the validity of using the adiabatic approximation for evolving the many-body system.

II The Neutrino Hamiltonian

We consider the evolution of an ensemble of $N$ free-streaming neutrinos undergoing vacuum and collective oscillations in two flavors: $e$ and $x$. We neglect inelastic collisions between these neutrinos, an approximation that may not be appropriate under all circumstances Cherry et al. (2012) but has the advantage of yielding a more tractable model. This system is described by the many-body Hamiltonian Pehlivan et al. (2011)

[TABLE]

where the vacuum oscillation frequencies are $\omega_{\mathbf{p}}=\delta m^{2}/2|\mathbf{p}|$ and $\mathbf{p}$ are the neutrino momenta. The strength of collective interactions for each pair of neutrinos is

[TABLE]

obtained from the leading-order effective four-point Fermi weak interaction diagrams of neutrinos in two flavors; $G_{F}$ is the Fermi coupling constant, $V$ is the volume in a box quantization, and $\theta_{\mathbf{pq}}$ is the angle between the momenta of interacting neutrinos. The weak isospin vectors $\vec{J}_{\mathbf{p}}$ for each neutrino are defined in terms of creation and annihilation operators of individual neutrinos in the mass basis: $J_{\mathbf{p}}^{z}=\frac{1}{2}(a_{1,\mathbf{p}}^{\dagger}a_{1,\mathbf{p}}^{\phantom{\dagger}}-a_{2,\mathbf{p}}^{\dagger}a_{2,\mathbf{p}}^{\phantom{\dagger}})$, $J_{\mathbf{p}}^{+}=a_{1,\mathbf{p}}^{\dagger}a_{2,\mathbf{p}}^{\phantom{\dagger}}$, and $J_{\mathbf{p}}^{-}=a_{2,\mathbf{p}}^{\dagger}a_{1,\mathbf{p}}^{\phantom{\dagger}}$. Here we define isospins up and down as corresponding to the neutrino vacuum mass eigenstates $\nu_{1}$ and $\nu_{2}$. The vector $\vec{B}$, which characterizes the oscillations of the individual neutrinos and plays a role similar to that of an external magnetic field if the $J$’s are interpreted as spins, takes the value $\vec{B}=(0,0,-1)$ in the mass basis. Note that each of the isospin operators and the vector $\vec{B}$ may also be expressed in the flavor basis, using a two-flavor mixing angle, $\theta$, unrelated to the $\theta_{\mathbf{pq}}$ in Eq. (2).

We study the tractable, but nevertheless interesting problem that is obtained by taking a geometric average over the angles $\theta_{\mathbf{pq}}$, called the single-angle approximation, which results in a single collective interaction strength as a function of position or time. For example, in a spherical neutrino bulb model Duan et al. (2006a, 2010),

[TABLE]

where $r$ is the distance from the center of a neutrino sphere of radius $R_{\nu}$. Here we adopt this form of $\mu(r)$ for the rest of our study. At this point, we can replace the $3$-vector momenta $\mathbf{p}$ with simple indices $p=1,\ldots,M$ to label the discrete energy bins in our model. For brevity, we henceforth also refer to $\mu(r)$ as simply $\mu$, except when the $r$ dependence is relevant for calculations. Thus, our Hamiltonian reduces to the simpler form

[TABLE]

where $\vec{J}=\sum_{p=1}^{M}\vec{J}_{p}$. This simpler Hamiltonian has several conserved quantities that may be helpful to consider, such as the total $z$-component isospin $J^{z}=\sum_{p=1}^{M}J_{p}^{z}$ and a set of $M$ conserved charges Pehlivan et al. (2011):

[TABLE]

Interestingly, the Hamiltonians in Eqs. (4) and (5) are also of interest in condensed matter physics, e.g., Refs. Owusu et al. (2008); Faribault et al. (2011); Claeys et al. (2017).

For simplicity, we choose a system where the $N$ neutrinos are distributed evenly across $M$ oscillation frequencies with $N=M$; therefore $j_{p}=1/2$ for $p=1,\ldots,N$. For this case, methods were developed to obtain all $2^{N}$ solutions, with eigenvalues and eigenstates of the Hamiltonian and these conserved quantities determined numerically for $N\geq 2$ (for details, refer to Appendix A or Refs. Claeys (2018); Cervia et al. (2019); Patwardhan et al. (2019)). These results are recorded as functions of $\mu/\omega_{0}$, with the oscillation frequencies fixed, $\omega_{p}=p\omega_{0}$; our value of $\omega_{0}$ is given in Sec. VII, Table 1.

III Adiabatic eigenstates and eigenvalues of the Hamiltonian

The Hamiltonian in Eq. (4) has eigenstates $\ket{q}$ where $q=0,1,\ldots,2^{N}-1$, which we may determine as superpositions of Kronecker products of mass basis states:

[TABLE]

where the $2^{N}$ coefficients $V_{i_{1},\ldots,i_{N}}(\mu)$ for each state are functions of $\omega_{p}$ and $\mu$.111Here, we make the dependence of these coefficients on $\mu$ explicit, since we change its value as time progresses. In contrast, we fix our values of $\omega_{p}$ at the beginning of our problem, and they do not change as time progresses. These coefficients have simple forms at $\mu=0$: $V_{i_{1},\ldots,i_{N}}(0)=\prod_{j=1}^{N}\delta_{i_{j}q_{j}}$, where $q=\sum_{j=1}^{N}(q_{j}-1)\,2^{N-j}$ is the binary representation of $q$ (i.e., $q_{j}-1$ is the $j$th digit from the left). That is, the energy eigenstates completely reduce to products of mass eigenstates of the individual neutrinos in the $\mu\to 0$ limit of the Hamiltonian in Eq. (4), as one would expect. For generic $\mu>0$, we begin with these solutions and numerically solve a system of $N$ coupled quadratic equations that constrain these coefficients. The relationship between these equations and coefficients [as well as eigenvalues of the charges in Eq. (5)] is described in Appendix A. From these equations one can find that all of these coefficients must be real valued. Therefore, the $2^{N}$ sets of $2^{N}$ coefficients compose a $2^{N}\times 2^{N}$ real matrix $V$ mapping energy states to mass states for a fixed $\mu$.

Suppose $\ket{\Psi}$ is the wave function of a many-body system with a Hamiltonian $H$, both in the mass basis. Then the time evolution of the system is given by

[TABLE]

The Hamiltonian may be diagonalized as $H=V\Sigma V^{\dagger}$ by a unitary transformation $V$, and one may define $\ket{\Psi}=V\ket{\Phi}$ to rewrite the equation as

[TABLE]

Assuming adiabaticity,222For a few initial conditions, we compared our adiabatic many-body results with the corresponding solutions obtained without assuming adiabaticity, and confirmed that they agree up to a reasonable level of precision (comparison not shown here). i.e., neglecting the second term on the left-hand side, one obtains an evolution equation for $\ket{\Phi}$, which may be integrated to obtain

[TABLE]

Since $\Sigma$ is a diagonal matrix, it is easy to exponentiate. In Appendix B we argue the time evolution operator used here does not need to be a Dyson series, as the commutator of Hamiltonians at different time values is $[H(t),H(t^{\prime})]\sim\mathcal{O}(t-t^{\prime})$ and therefore is numerically insignificant for small time steps. The equation for $\ket{\Psi}$ then becomes

[TABLE]

where $V$ and $V_{0}$ are the unitary transformations that diagonalize the Hamiltonian at the values $\mu=\mu(r)$ and the initial value $\mu_{0}\equiv\mu(r_{0})$, respectively, with $r_{0}\geq R_{\nu}$, where $R_{\nu}$ is the neutrino sphere radius as introduced in Eq. (3). Starting from an initial state $\ket{\Psi(\mu_{0})}$ of our many-body neutrino system,the adiabatic evolution of the system is then given by

[TABLE]

where $\Sigma_{qq}\equiv\braket{q}{\Sigma}{q}$. Note that here we have used the position $r$ as a proxy for time $t$ in the evolution equation, as is justified if we assume a steady-state configuration wherein the interaction strength $\mu$ depends explicitly only on position and not time.

It is worth pointing out, as observed in Refs. Birol et al. (2018); Patwardhan et al. (2019), that energy eigenvalues of the Hamiltonian in Eq. (4) exhibit numerous energy level crossings. However, in our previous study Patwardhan et al. (2019), it was determined that the conserved charge operators in Eq. (5) cannot all be simultaneously degenerate for any given value of $\mu$. Therefore, as is argued in Appendix B, we propose that one can adiabatically evolve this state without encountering nontrivial level crossings to be dealt with, as these charges break any degeneracies that may arise in our Hamiltonian.

IV Quantifying entanglement

In this section we review the entanglement entropy and the Bloch vector, which are well-established useful concepts in the context of quantum information theory Everett (1956); Nielsen and Chuang (2011). For an $N$-particle wave function $\ket{\Psi}$ defined on a Hilbert space $\mathcal{H}=\mathcal{H}_{1}\otimes\cdots\otimes\mathcal{H}_{N}$, one can define the density matrix $\rho=\ket{\Psi}\!\bra{\Psi}$. Partitioning the Hilbert space as $\mathcal{H}=\mathcal{H}_{A}\otimes\mathcal{H}_{B}$, one can obtain the reduced density matrix for partition $A$ by taking the partial trace over partition $B$,

[TABLE]

A suitable measure of entanglement for nonmixed states is the bipartite entanglement entropy. For partitions $A$ and $B$, the entanglement entropy is defined as

[TABLE]

In this work, we choose one of the neutrinos to be in partition $A$, and the remaining neutrinos to be in partition $B$, so that the reduced density matrix $\rho_{A}$ is simply a $2\times 2$ matrix. For example, if we choose the neutrino at the highest $\omega$ to be in partition $A$, then $\rho_{A}$ is given by

[TABLE]

One may define the “polarization vector” for each particle as $\vec{P}_{\omega}=2\braket{\Psi}{\vec{J}_{\omega}}{\Psi}$. If $\mathcal{H}_{A}$ contains a single neutrino, its polarization vector $\vec{P}_{A}$ is related to its reduced density matrix as

[TABLE]

where $\mathbf{I}$ is the $2\times 2$ identity matrix and $\vec{\sigma}$ are the Pauli spin matrices. For the sake of brevity, we drop the subscripts “$A$” and “$(A|B)$” from the polarization vector and entropy, respectively, in the remainder of the text. In this scenario, the polarization vector is simply the standard $SU(2)$ Bloch vector. It is straightforward to check that the eigenvalues of $\rho_{A}$ are $(1\pm P)/2$, where $P\equiv|\vec{P}|$. One can then relate the magnitude $P$ of the polarization vector to the entropy of entanglement via the relation

[TABLE]

In this sense, $P$ may be used as a measure of entanglement; $P=1$ indicates a pure state, whereas $P<1$ implies that the reduced density matrix $\rho_{A}$ is mixed, signaling entanglement between the two partitions. Additionally, one can observe that $P$ is related to the concurrence of the state across the partitions $A$ and $B$, via $C=\sqrt{1-P^{2}}$ Wootters (1998); Hill and Wootters (1997).

V Overview of the results

Here we examine the differences between the predictions of mean-field theory and many-body results obtained for small $N$ and also show that the magnitude of these differences is strongly correlated with measures of quantum entanglement in the many-body solutions. To explore differences between the mean-field approximation and the adiabatic many-body solutions we calculated several quantities in both approaches. Details of these calculations are given in subsequent sections. In this section we provide an overview of our results. The key results that we obtain are the following:

The magnitude of the differences between the results obtained using the mean-field approximation and those obtained using the many-body approach are correlated with measures of entanglement of the asymptotic configuration reached in the limit of long times, and 2. 2.

Deviations between the many-body solutions and mean-field theory for asymptotic values of observables depend strongly on the initial condition and often are substantial.

To quantify physically relevant differences between the mean-field and many-body solutions, in Fig. 1a we compare asymptotic values (i.e., at large radius $r$ where $\mu\ll\omega_{0}$) of the $z$ component of the polarization vector, $P_{z}$, for the neutrino with the highest frequency $\omega_{N}=N\omega_{0}$. The asymptotic values are calculated in both the mean-field and the many-body approaches, for different values of $N$, the number of neutrinos. In this figure the initial state is a state with $N$ electron neutrinos. For this initial configuration, the deviation between the many-body and the mean-field results increases with the number of neutrinos. Furthermore, in Fig. 1b we present the evolution of the entropy of entanglement $S$ of the highest frequency neutrino, with the rest of the ensemble [calculated using Eqs. (13) or (16)]. Note that the entanglement entropy is always zero in the mean-field approximation. The growth in asymptotic values of $S$ with $N$ demonstrates that the entanglement entropy and the magnitude of the discrepancy of $P_{z}$ between the exact and mean-field solutions are correlated.

We next present results demonstrating that the many-body quantum correlations can vary greatly between different initial configurations. To illustrate this in Fig. 2 we display asymptotic ($r\gg R_{\nu}$) values of the difference

[TABLE]

$p=1,\ldots,N$, where MF and MB denote mean-field and many-body results, respectively, for the $z$ component of the polarization vector of the neutrino with frequency $\omega_{p}=p\omega_{0}$, versus the entanglement entropy for this neutrino with the rest of the ensemble in the many-body calculation. The left panel is for $N=4$ neutrinos for all 16 ($2^{4}$) possible initial state configurations with neutrinos of definite flavor.333A global inversion of both flavor and the ordering of $\omega$ values of the evolved state results in another solution for a different initial configuration, with $\vec{P}(\omega_{p})\mapsto-\vec{P}(\omega_{N-p})$. Therefore, each point appearing in Fig. 2a, in fact, represents the same data for two different initial configurations of overall opposite flavor and reversed ordering; for example, initial states $\ket{\nu_{e}\nu_{x}\nu_{x}\nu_{x}}$ and $\ket{\nu_{e}\nu_{e}\nu_{e}\nu_{x}}$ correspond to the same four data points $\{S(\omega_{p}),\Delta P_{z}(\omega_{p})\}$ for $p=1,\ldots,4$. The right panel is for an $N=8$ neutrino system, for a subset of configurations where four of the neutrinos, placed at either $\omega_{1},\ldots,\omega_{4}$ or $\omega_{5},\ldots,\omega_{8}$ are iterated through the same 16 initial configurations as in the left panel, whereas the remaining ones are all taken to be in the $\nu_{e}$ flavor initially. This figure illustrates that, whenever a large deviation from mean-field theory is exhibited, the entanglement entropy tends to be large. However, the converse is not necessarily true—in some cases, even when the entanglement entropy is large, the deviations from mean-field theory can nevertheless be small. The lines drawn in Fig. 2 are of the quantity

[TABLE]

where $P(S)$ is the inverse of the function in Eq. (16). We observe that a substantial number of $\Delta P_{z}(\omega_{p})$ values cluster around this line, particularly for $N=8$. This alignment could occur, for instance, when $P_{x}$ and $P_{y}$ components of these vectors are vanishingly small.

VI Analytic results for a two-neutrino system

Before we elaborate further on the results overviewed in Sec. V, it is instructive to illustrate the example of the two-neutrino system. This system is simple enough to be examined analytically, but it nevertheless brings to light some key features that are generically present in many-body systems. In Sec. VI.1, we describe the adiabatic evolution of the $N=2$ neutrino system starting from the different possible initial conditions, and in Sec. VI.2, for comparison, we present the mean-field evolution equations for the same system.

VI.1 Many-body results in the adiabatic limit

For two neutrinos, the Hamiltonian in the product mass basis can be written as

[TABLE]

where $\Omega\equiv(\omega_{1}+\omega_{2})/2$ and $\eta\equiv(\omega_{1}-\omega_{2})/2$. This can be diagonalized into

[TABLE]

where we have also subtracted a term $\mu\mathbb{I}_{4\times 4}$, which contributes only a global phase to the evolved state and therefore is not present in calculations with the resulting density matrix of the state. Here, the unitary transformation matrix can be parametrized as

[TABLE]

where one has

[TABLE]

One can thus observe that only the second and the third components of the wave function mix. Writing the wave function as $\ket{\Psi(t)}=(\alpha(t),\beta(t),\gamma(t),\delta(t))^{T}$, one obtains the following evolution equations in the adiabatic limit:

[TABLE]

where $\alpha_{0},\ldots,\delta_{0}$ and $\xi_{0}$ are the respective initial values, and where $R=\int^{t}_{t_{0}}\mu\,dt^{\prime}$ and $Q=\int^{t}_{t_{0}}\sqrt{\mu^{2}+\eta^{2}}\,dt^{\prime}$. The reduced density matrix for the second particle can be written as

[TABLE]

For this subsystem, one can use trigonometric identities to obtain an expression for $P_{z}=(\rho_{2}^{\text{red}})_{11}-(\rho_{2}^{\text{red}})_{22}$,

[TABLE]

One can, in principle, also obtain analytic expressions for $P_{x}=2\,\mathfrak{Re}[(\rho_{2}^{\text{red}})_{12}]$ and $P_{y}=-2\,\mathfrak{Im}[(\rho_{2}^{\text{red}})_{12}]$. Of particular interest, however, is the magnitude $P$, since it is directly related to the entanglement entropy, as described previously with Eq. (16). First, using Eq. (27), one can write down an expression for $P_{\perp}=|P_{x}-iP_{y}|=2|(\rho_{2}^{\text{red}})_{12}|$, given by

[TABLE]

These may then be combined with $P_{z}$ to obtain a succinct expression for the magnitude $P$, given by

[TABLE]

Alternatively, this expression can be obtained by recognizing that the determinant of $\rho_{B}$ is the product of its eigenvalues, and is therefore given by $(1-P^{2})/4$ (as per our discussion in Sec. IV).

One can now examine the evolution of this system starting from different initial conditions. For a system where the initial state is $|\nu_{e}\nu_{e}\rangle$ or $\ket{\nu_{x}\nu_{x}}$, one has $\beta_{0}=\gamma_{0}=\pm\cos\theta\sin\theta$ and $\{\alpha_{0},\delta_{0}\}=\{\cos^{2}\theta,\sin^{2}\theta\}$ or $\{\sin^{2}\theta,\cos^{2}\theta\}$ (respectively for the two distinct initial conditions). Using this condition and the expressions for $\cos 2\xi$ and $\sin 2\xi$ from Eq. (22) leads to

[TABLE]

In particular, if $\mu,\mu_{0}\gg|\eta|$, then for small $\theta$ one can observe that $P_{z}$ oscillates with an amplitude of $\frac{\eta}{2\mu_{0}}\sin^{2}2\theta\ll 1$, therefore oscillations are insignificant in this case. Interestingly, if we further take the limit $\eta\to 0$ while $\mu,\mu_{0}>0$, we find that $P^{2}\to 1$ and so the entanglement entropy $S\to 0$. In contrast, for an initial state $|\nu_{e}\nu_{x}\rangle$ or $\ket{\nu_{x}\nu_{e}}$, one has $\alpha_{0}=-\delta_{0}=-\cos\theta\sin\theta$, and $\{\beta_{0},\gamma_{0}\}=\{\cos^{2}\theta,-\sin^{2}\theta\}$ or $\{-\sin^{2}\theta,\cos^{2}\theta\}$ (respectively for the two distinct initial conditions), which gives

[TABLE]

Here, in the limit $\mu,\mu_{0}\gg|\eta|$, the oscillation amplitude approaches $\cos 2\theta$, which is $\sim 1$ for small $\theta$. Moreover, the instantaneous oscillation frequency is $2\,dQ/dt=2\sqrt{\mu^{2}+\eta^{2}}\approx 2\mu$, which indicates that the oscillations are fast (compared to the oscillations observed in the mean-field limit, which have a frequency $\sim\sqrt{|\eta|\mu}$—see, e.g., Ref. Hannestad et al. (2006)). Fast oscillations of neutrinos are also observed in certain mean-field calculations in the presence of spatial asymmetries or nonstandard neutrino-neutrino interactions Sen (2018); Dighe and Sen (2018); however, in this model, the presence of many-body effects appears to render these conditions unnecessary. Fast oscillations arising from many-body effects have also been pointed out in Ref. Rrapaj , and we plan to investigate this feature in a future study. In a core-collapse supernova envelope, fast oscillations can have implications for nucleosynthesis as well as for the explosion mechanism itself.

Additionally, we may consider the asymptotic forms of $P_{z}$ and $P$ for $\mu\to 0$ in the cases of each initial condition. Here, our $z$-component polarization reduces to

[TABLE]

for an initial state of $\ket{\nu_{e}\nu_{e}}$ or $\ket{\nu_{x}\nu_{x}}$. Meanwhile, from the initial states $\ket{\nu_{e}\nu_{x}}$ or $\ket{\nu_{x}\nu_{e}}$ we find

[TABLE]

Interestingly, the forms for $P_{z}$ depend on only the initial and the final values of $\mu$ ($\mu_{0}$ and [math], in this case), while the magnitudes $P^{2}$ also depend on all the intermediate values of $\mu$, i.e, $0<\mu<\mu_{0}$, through terms proportional to $\cos 2R$. If we also take $\mu_{0}\gg|\eta|$, this additional dependence vanishes for the monoflavor initial states (i.e., $\ket{\nu_{e}\nu_{e}}$ and $\ket{\nu_{x}\nu_{x}}$) but not for the mixed-flavor initial states.

Finally we note that, for a two-neutrino system, the results obtained for $P^{2}$ must be independent of our choice of particle to take the partial trace over; $P$ shares a one-to-one relationship with entanglement entropy $S$, as per Eq. (16), which for a bipartite pure-state system does not depend on which of the two partitions is being traced over. However, this invariance is not shared by $P_{z}$, which will differ between the two neutrinos’ subsystems. To explicitly see this fact, one can show that the first particle’s $z$ component of polarization is given by $P_{z}=|\alpha|^{2}+|\beta|^{2}-|\gamma|^{2}-|\delta|^{2}$, in contrast with the form for the second particle in Eq. (VI.1).

VI.2 An example of the mean-field limit

To illustrate the mean-field limit of the Hamiltonian in Eq. (4), here we show a simple example with two neutrinos. It is instructive to subtract constant terms proportional to $\vec{J}_{1}^{2}$ and $\vec{J}_{2}^{2}$, to obtain

[TABLE]

Expressing the average of an operator $X$ over the mean field as $\langle X\rangle_{M}$, we can write the quadratic term in the above equation as

[TABLE]

where one observes that the Hamiltonian in Eq. (4) becomes, up to a constant term, a direct sum of Hamiltonians describing individual neutrinos

[TABLE]

With this approximation, the Hamiltonian can be written as a direct sum of the individual one-body Hamiltonians. Consequently, the eigenstates are direct products of the single neutrino states, and therefore unentangled. The time evolution of these operators takes the form

[TABLE]

and a similar equation for $\vec{J}_{2}$, where the indices $i,j,k$ represent the Cartesian components. Note that a consistent mean field requires that Eq. (41) imply

[TABLE]

and a similar equation for $\vec{J}_{2}$. Defining the mean-field polarization vectors as $\vec{P}_{1}=2\langle\vec{J}_{1}\rangle_{M}$, $\vec{P}_{2}=2\langle\vec{J}_{2}\rangle_{M}$, and $\vec{P}=\vec{P}_{1}+\vec{P}_{2}$, we obtain

[TABLE]

and

[TABLE]

From these equations we see that both $|\vec{P}_{1}|$ and $|\vec{P}_{2}|$ are constants in time (and equal to one for this choice of normalization), as expected for unentangled particles. Additionally, by adding the two equations, it is easy to see that $\vec{B}\cdot\vec{P}$ is also conserved, which is consistent with the commutation of total $J^{z}$ with the Hamiltonian in Eq. (4).

In Figs. 3 and 4, we show the evolution of $P_{z}$ of neutrino 2 as a function of the radius, for an $N=2$ system with initial conditions $\ket{\nu_{e}\nu_{e}}$ and $\ket{\nu_{e}\nu_{x}}$, respectively. In addition, we also show the dependence of entropy on $\mu$ for each of the systems, alongside the corresponding $P_{z}$ evolution plots. Here, $\mu$ is taken to vary with the radius in accordance with the single-angle bulb model expression, given in Eq. (3). A comparison of the two figures reveals that for the system with a $\ket{\nu_{e}\nu_{x}}$ initial condition, the entanglement entropy is quite high, and this is associated with the evolution of $P_{z}$ differing substantially in the many-body vs mean-field treatments. In contrast, for the system starting from a $\ket{\nu_{e}\nu_{e}}$ configuration, the entropy never grows too large, and the mean-field evolution does not deviate much from that of the many-body solution.

VII Numerical results for more than two neutrinos

In this section we quantify the entanglement between neutrinos in time-evolved many-body states described in Section III. For this purpose, we choose systems with different initial configurations, wherein each neutrino is either in the $\nu_{e}$ or in the $\nu_{x}$ flavor, and we apply the methods described in Sec. IV to calculate measures of entanglement for these systems.

To determine this entanglement entropy at each $\mu<\mu_{0}$ (or equivalently $r>r_{0}$), we first construct a density matrix from our evolved state, in the mass basis.444One can just as well carry out this process using the evolved state in the flavor basis. The only change to our results will be in the orientation of the Bloch vector between the two bases, for $\theta\neq 0$. For calculating entanglement between the neutrino with the highest value of $\omega$ and the other $N-1$ neutrinos, we compute a reduced density matrix by taking a partial trace over the subspaces of the $N-1$ other neutrinos, as per Eq. (14). With this effective one-body density matrix we can calculate the entropy of entanglement between this neutrino and the rest. Results found from this procedure are displayed in Fig. 1b. As we expect, $S\sim 0$ in the limit of small $r\sim R_{\nu}$. Further, as $r$ grows ($\mu$ decreases), we find that entropy values eventually level off to constant values as the collective oscillation strength becomes much smaller than vacuum oscillation frequencies (i.e., $\mu\ll\omega_{0}$). For the results displayed in this paper, we use typical values for the parameters mentioned thus far, summarized in Table 1.

We investigated other measures of entanglement for these evolved states in addition to the entanglement entropy. For example, we calculated logarithmic negativity of the subsystem for the neutrino with vacuum oscillation frequency $\omega_{N}$, as prescribed in Ref. Vidal and Werner (2002). Additionally, we calculated the entanglement of formation, as defined by Ref Wootters (1998), of the effective mixed, two-neutrino state corresponding with oscillation frequencies $\omega_{N-1}$ and $\omega_{N}$ after taking the partial trace over all $N-2$ other neutrinos. We verified that both measures yield qualitatively similar behavior to that of entanglement entropy and are therefore not presented in this article.

Moreover, we seek to determine whether many-body effects can drastically alter the evolution of a multineutrino system, as compared with the evolution observed in the mean-field limit. To this end, we plot the evolution of $P_{z}$ of each neutrino in both the many-body and the mean-field picture.

As described in Sec. VI, one can obtain analytic closed-form expressions for the evolving quantities for an $N=2$ system, which may then be used for comparison with the corresponding mean-field solution. Such a comparison is presented in Figs. 3 and 4. We can repeat this comparison for systems with larger $N$ using numerical methods as outlined in Sec. III and Appendix A. Shown in Figs. 5 and 6 are $P_{z}$ (in both many-body and mean-field calculations) and $S$ results of each neutrino for a system with $N=8$, with two different initial states. In Fig. 5, we depict a system starting from an initial state with neutrinos with frequencies $\omega_{1},\ldots,\omega_{7}$ in the $\nu_{e}$ flavor, and the $\omega_{8}$ neutrino in the $\nu_{x}$ flavor, whereas in Fig. 6, we pick an initial state consisting of $\nu_{e}$ at frequencies $\omega_{1},\ldots,\omega_{4}$, and $\nu_{x}$ at $\omega_{5},\ldots,\omega_{8}$. Furthermore, we illustrate in Figs. 7 and 8 how entanglement entropy can reflect the discrepancies between many-body and mean-field predictions of time-evolved mass-basis spectra of the ensemble, for the same initial states as in Figs. 5 and 6 respectively. In the final spectra plots comparing mean-field and many-body predictions, Figs. 7a and 8a, we observe a weaker spectral swap when including many-body effects, as asymptotic values of $P_{z}$ at the frequencies surrounding the swap do not approach $\pm 1$, unlike their counterparts in the mean-field theory. Further, as previously outlined in Sec. V, there appears to be a correlation, illustrated in Fig. 2, between the entanglement entropy of each neutrino with the ensemble and the deviation of its polarization vector $z$ component, $P_{z}$, relative to the mean-field limit.

VIII Conclusions

Collective neutrino oscillations represent a nonlinear quantum problem with a significant impact on astrophysics and cosmology. However, calculating collective oscillations of a very large number of neutrinos present in these settings is a challenging many-body problem. To simplify the problem one often resorts to the mean-field approximation albeit at the cost of losing entanglement between individual neutrinos. Hence, measures of entanglement could be used to quantify the deviation of the mean-field approximation from many-body results. We showed that the entanglement entropy, which is zero within the mean-field approximation but nonzero in the many-body calculations, is a useful quantity for identifying regimes where mean-field theory may be insufficient.

In this paper, we examined collective neutrino oscillations for systems with small numbers of neutrinos, for which we can easily obtain many-body solutions using the methods described in Ref. Patwardhan et al. (2019), and we showed that the mean-field approximation can be incomplete. We also demonstrated that, for systems consisting initially of a single flavor, the third component of the polarization vector for the neutrino with the highest frequency deviates more and more from the mean-field value as the number of neutrinos is increased from $N=2$ to $9$. Correspondingly the asymptotic value of the entanglement entropy of those neutrinos increases as the number of neutrinos increases.

Our results exhibit a strong dependence on the initial conditions. We also examine various initial conditions, consisting of both $\nu_{e}$ and $\nu_{x}$ flavors, where the evolution in the mean field results in a spectral swap. In the many-body calculations, we observe that the swap is largely preserved, but appears weaker because of the $P_{z}$ values surrounding the swap frequency being further from $\pm 1$ compared to the mean-field case. This deviation is also correlated with the entanglement entropy of those neutrinos being large.

Calculations performed using the mean-field approximation have revealed a lot of interesting physics about collective behavior of neutrinos in astrophysical environments. Here we have explored possible scenarios where further interesting features can arise by going beyond this approximation.

Acknowledgements.

We thank Y. Pehlivan, E. Rrapaj, and P. Claeys for helpful conversations. This work was supported in part by the U.S. Department of Energy, Office of Science, Office of High Energy Physics, under Award No. DE-SC0019465. It was also supported in part by the U.S. National Science Foundation Grants No. PHY-1630782 and No. PHY-1806368.

Appendix A Determination of Eigenvalues and Eigenstates

In this section, we outline how to determine the eigenvalues and eigenstates of the Hamiltonian and the conserved charge operators described in Sec. II. A more detailed exposition can be found in our previous study, Ref. Patwardhan et al. (2019)—however, here we point out, at the end of our summary, a symmetry in the eigenvalue equations that can be used to reduce the computational time. Though arguments in this appendix can be generalized to the higher dimensional (i.e., $j_{p}>1/2$) representations, in this discussion we mainly focus on the case of $j_{p}=1/2$ for all $p$, which corresponds to having a single neutrino at each $\omega_{p}$.

Because the operator $J^{z}$ (the $z$ component of the total isospin) commutes with the Hamiltonian, each energy eigenstate will have a definite eigenvalue $m=-N/2,-N/2+1,\ldots,+N/2$ of $J^{z}$. For a given $N$ and $m$ it can be shown that there are ${}^{N}C_{N/2-m}$ energy states with total $z$-component isospin number $m$. Since this quantity is conserved even in the event of degeneracies in energy eigenvalues Birol et al. (2018), it is helpful to categorize solutions to our system by their corresponding value of $m$. Additionally, the solutions for eigenvalues of conserved charges $h_{p}$ allow us to easily determine the $m$ value of a given state.

It is possible to obtain an algebraic system of coupled equations for the eigenvalues $\epsilon_{p}$ of the charge operators $h_{p}$, for any $j_{p}$. For $j_{p}=1/2$ we find there are $N$ coupled quadratic equations constraining these values Dimo and Faribault (2018):

[TABLE]

The eigenvalues $\epsilon_{p}$ can then be used to obtain the energy eigenvalues of the Hamiltonian:

[TABLE]

Further, these charge eigenvalues are related to another system of $N$ parameters $\Lambda_{p}$ Balantekin (2018):

[TABLE]

These parameters similarly satisfy a system of $N$ coupled quadratic equations for $j_{p}=1/2$ Babelon and Talalaev (2007); Faribault et al. (2011):

[TABLE]

where $\vec{\Lambda}\equiv(\Lambda_{1},\ldots,\Lambda_{N})$. The choice of alternative sign in Eqs. (47) and (48) is determined by the choice of formalism (i.e., raising or lowering, respectively) for defining these parameters as per the Bethe ansatz equations Claeys et al. (2017); Claeys (2018); Birol et al. (2018); Patwardhan et al. (2019). With these $N$ parameters one is equipped to determine the overlaps (up to normalization) of all energy eigenstates $\ket{q}$ with the mass product basis elements $\ket{\nu_{i_{1}},\ldots,\nu_{i_{N}}}$:

[TABLE]

where $\mathcal{J}$ is a $\kappa\times\kappa$ matrix, with $\kappa=N/2\pm m$, defined by

[TABLE]

with indices $p_{\alpha}$ only spanning values satisfying $i_{p_{\alpha}}=1$ ($2$) Claeys et al. (2017); Claeys (2018). That is, this matrix only considers the frequency bins $\omega_{p_{\alpha}}$ occupied by a neutrino with isospin-up (-down) in the state $\ket{\nu_{i_{1}},\ldots,\nu_{i_{N}}}$. Moreover, elements of this matrix $\mathcal{J}$ are evaluated at the particular solution $\vec{\Lambda}$ to Eqs. (48) that satisfies the condition $\pm\mu\Lambda_{j}\to q_{j}$ as $\mu\to 0$ for $j=1,\ldots,N$, where $q=\sum_{j=1}^{N}(q_{j}-1)\,2^{N-j}$ is the binary representation of $q$. Alternatively, the eigenstates (again, up to a normalization factor) can be constructed from the identities

[TABLE]

where $y=2$ ($1$) in the raising (lowering) formalism, $\mathrm{Sym}(\kappa)$ is the symmetry group of $\kappa$ letters, $\mathrm{Tr}_{\sigma}(A)=\mathrm{Tr}(A^{f_{1}})\cdots\mathrm{Tr}(A^{f_{n}})$ for a permutation $\sigma$ of cycle type $(f_{1},\ldots,f_{n})$, and the individual traces are

[TABLE]

Here, $A$ is a diagonal matrix of Gaudin raising (lowering) operators, as defined in Ref. Patwardhan et al. (2019). One can determine the corresponding $m$ for a solution using either $\vec{\Lambda}$ or $\vec{\epsilon}$, as $\sum_{p}\mu\Lambda_{p}=\pm\kappa$ and $\sum_{p}\epsilon_{p}=-m$. For this reason as well, values of $m$ are convenient for categorizing solutions.

Note that from a solution $\vec{\Lambda}$ obtained in a raising (lowering) formalism one can obtain another solution in the lowering (raising) formalism, given by $\vec{\Lambda}\mp(1/\mu,\ldots,1/\mu)$ Claeys et al. (2017); Patwardhan et al. (2019). Furthermore, by selecting either the raising or the lowering formalism for determining an eigenstate so that $\kappa=\min\{N/2+m,N/2-m\}$, considerable computation time may be reduced in calculating these states in the mass basis.

Using the methods outlined in Sec. IV, we may compute the entanglement entropy between any given neutrino and the rest of the ensemble, in any particular energy eigenstate $\ket{q}$, for generic values of $0<\mu<\mu_{0}$. First, we may express the state in either the mass or the flavor basis—here, we use the mass-basis density matrix $\rho(\mu)=V\ket{q}\!\bra{q}V^{T}$, and then reduce it by taking the partial trace over the remaining $N-1$ neutrinos. For example, to obtain an entropy of entanglement between the neutrino with the highest frequency and the rest of the ensemble for a given energy eigenstate, we may use Eq. (14) for the reduced density matrix and Eq. (13) for the entropy. Results for energy eigenstates of a $N=5$ system are displayed in Fig. 9. As one might expect, in the limit of $\mu\to 0$, one obtains $S\to 0$, as the energy eigenstates reduce simply to direct products of individual neutrino mass eigenstates. Notably, states with $m=\pm N/2$ (not shown in the figure) remain as direct products even for generic $\mu>0$, and therefore continue to have $S=0$.

Appendix B Conditions for Adiabaticity

It was noted in Refs. Birol et al. (2018); Patwardhan et al. (2019) that the Hamiltonian in Eq. (4) has an energy spectrum exhibiting numerous level crossings. Some of these crossings are between energies of eigenstates with differing eigenvalues $m$ of the total $z$-component isospin $J^{z}$. Since $J^{z}$ is a time-independent conserved quantity of the Hamiltonian, one can show using Ehrenfest’s theorem that the expectation value of $J^{z}$ with respect to any state is preserved as the state evolves through such a level crossing. Therefore, there can be no mixing between states with different eigenvalues $m$, thus rendering these crossings unimportant in the evolution of the system.

However, there are additional level crossings between energy eigenvalues of states with identical values of $m$. In this appendix, we

justify approximating our time evolution operator as $\mathcal{T}\{e^{-i\int_{0}^{t}H(t^{\prime})dt^{\prime}}\}\approx e^{-i\int_{0}^{t}H(t^{\prime})dt^{\prime}}$ from the analyticity of the Hamiltonian in time, and then 2. 2.

use the behavior of the commutators $[h_{p}(\mu),h_{q}(\mu^{\prime})]$ to show that differing eigenvalues in $\vec{h}\equiv(h_{1},\ldots,h_{N})$ prevent mixing of degenerate states at the level crossings, including the states with the same $m$ values.

It is helpful in demonstrating both claims listed above to first outline calculations for the commutators $[h_{p}(\mu),h_{q}(\mu^{\prime})]$. Using the definition of these charges in Eq. (5) and the fact

[TABLE]

it can be shown that our desired commutators are given by

[TABLE]

Furthermore, using the identity

[TABLE]

it can also be shown that

[TABLE]

Notably, in the limit $\delta\mu\to 0$, Eqs. (54) and (56) reduce to the well-known commutation relations between the conserved charges.

Recall that a Magnus expansion of the Dyson series is given by

[TABLE]

for $\int_{t}^{t+\delta t}||H(t^{\prime})||_{2}dt^{\prime}<\pi$ with the matrix norm ||X||_{2}\equiv\big{[}\sum_{i,j}|X_{ij}|^{2}\big{]}^{1/2} Blanes et al. (2009). Since our Hamiltonian is analytic in time, we can use the fact that $[H(t),H(t^{\prime})]\sim\mathcal{O}(\delta t)$, for sufficiently small time steps $\delta t\equiv t^{\prime}-t$, to approximate this expansion by the lowest order term in $\delta t$—as all terms with commutators will be $\lesssim\mathcal{O}[(\delta t)^{3}]$. This commutator relationship is also corroborated by Eq. (57). Combining these approximate terms $\mathcal{T}\left\{e^{-i\int_{t}^{t+\delta t}H(t^{\prime})dt^{\prime}}\right\}\approx e^{-i\int_{t}^{t+\delta t}H(t^{\prime})dt^{\prime}+\mathcal{O}[(\delta t)^{3}]}$ for all time steps between [math] and $t$, we obtain

[TABLE]

where $N_{t}\sim t/\delta t$. This argument can similarly be applied to taking small steps in $\delta\mu$, as these steps can be related to small time steps, given a function $\mu(t)$ such as Eq. (3) with a well-behaved derivative, $d\mu/dt$.555While our form for $d\mu/dt$ implies that $\delta t$ will grow as $\mu\to 0$ if step sizes $\delta\mu$ are held fixed, in this limit we find that the nonlinear term of our Hamiltonian in Eq. (4) becomes negligible, in which case our Hamiltonian would commute with itself at varying times and make our approximate time evolution operator exact regardless. Thus, we propose the approximation of the time evolution operator used in Eq. (9), justifying claim 1 stated earlier.

Next we turn our attention to claim 2. Consider an energy eigenstate with charge eigenvalues $\vec{\epsilon}(t)$ for times $0\leq t\leq t_{d}$, where $t_{d}$ is a time at which this state shares the same energy eigenvalue with another eigenstate of $\vec{h}$. The time evolution of this state to values beyond $t_{d}$ over a small time interval $(t_{d}-\delta,t_{d}+\delta)$ can be described by \exp\big{[}-i\int_{t_{d}-\delta}^{t_{d}+\delta}H(t^{\prime})dt^{\prime}\big{]}\ket{\vec{\epsilon}(t_{d}-\delta)}; moreover, applying the operator $h_{p}$ to this state yields

[TABLE]

where we use the approximation $h_{p}(t_{d}+\delta)\ket{\vec{\epsilon}(t_{d}-\delta)}\approx\epsilon_{p}(t_{d}-\delta)\ket{\vec{\epsilon}(t_{d}-\delta)}+\mathcal{O}(\delta)$.

Furthermore, we can observe for the latter term of Eq. (60) that

[TABLE]

Using our earlier calculation in Eq. (56), we can see that this term is $\mathcal{O}(\delta^{2})$ and therefore may be neglected for $\delta\to 0$. Taking this limit, we thus explicitly see that states with different eigenvalues $\vec{\epsilon}$ of $\vec{h}$ are prevented from mixing with one another as they evolve through the crossing.