Eigenvalues and eigenstates of the many-body collective neutrino oscillation problem

TL;DR

This paper introduces a novel numerical method based on the Richardson-Gaudin framework to systematically compute eigenvalues and eigenstates of a many-body Hamiltonian describing collective neutrino oscillations, simplifying the solution process.

Contribution

The paper develops an auxiliary variable approach that transforms Bethe Ansatz equations into more tractable forms for analyzing many-body neutrino oscillations.

Findings

Method successfully computes eigenvalues and eigenstates.

Approach improves numerical stability and tractability.

Applicable to two-flavor, single-angle neutrino systems.

Abstract

We demonstrate a method to systematically obtain eigenvalues and eigenstates of a many-body Hamiltonian describing collective neutrino oscillations. The method is derived from the Richardson-Gaudin framework, which involves casting the eigenproblem as a set of coupled nonlinear "Bethe Ansatz equations", the solutions of which can then be used to parametrize the eigenvalues and eigenvectors. The specific approach outlined in this paper consists of defining auxiliary variables that are related to the Bethe-Ansatz parameters, thereby transforming the Bethe-Ansatz equations into a different set of equations that are numerically better behaved and more tractable. We show that it is possible to express not only the eigenvalues, but also the eigenstates, directly in terms of these auxiliary variables without involving the Bethe Ansatz parameters themselves. In this paper, we limit ourselves to a two-flavor, single-angle neutrino system.

Figures

Tables

Table 1. Table 1: Algebraic solutions for Λ p subscript Λ 𝑝 \Lambda_{p} from Eq. ( 22 ) with M = 2 𝑀 2 M=2 .

	${\Lambda }_1$	${\Lambda }_2$
$\kappa =0$	0	0
$\kappa =1$	$-\frac{1}{2\mu }+\frac{1}{{\eta }_{12}}-\text{sgn}\left({\eta }_{12}\right)\sqrt{\frac{1}{{\eta }_{12}^2}+\frac{1}{4{\mu }^2}}$	$-\frac{1}{2\mu }-\frac{1}{{\eta }_{12}}+\text{sgn}\left({\eta }_{12}\right)\sqrt{\frac{1}{{\eta }_{12}^2}+\frac{1}{4{\mu }^2}}$
	$-\frac{1}{2\mu }+\frac{1}{{\eta }_{12}}+\text{sgn}\left({\eta }_{12}\right)\sqrt{\frac{1}{{\eta }_{12}^2}+\frac{1}{4{\mu }^2}}$	$-\frac{1}{2\mu }-\frac{1}{{\eta }_{12}}-\text{sgn}\left({\eta }_{12}\right)\sqrt{\frac{1}{{\eta }_{12}^2}+\frac{1}{4{\mu }^2}}$
$\kappa =2$	$-\frac{1}{\mu }$	$-\frac{1}{\mu }$

Table 2. Table 2: Energy eigenvalues and eigenstates determined from Eq. ( 22 ) M = 2 𝑀 2 M=2 .

	$E$	$|E⟩$
$\kappa =0$	$2\mu -\frac{1}{2}\left({\omega }_1+{\omega }_2\right)$	$|{\nu }_1,{\nu }_1⟩$
$\kappa =1$	$-\left({\omega }_1+{\omega }_2\right)-\mu \left[1+\sqrt{1+\frac{{\left({\omega }_1-{\omega }_2\right)}^2}{4{\mu }^2}}\right]$	$\frac{1}{{𝒩}_{1,+}}\left({\Lambda }_{1,+}|{\nu }_2,{\nu }_1⟩+{\Lambda }_{2,+}|{\nu }_1,{\nu }_2⟩\right)$
	$-\left({\omega }_1+{\omega }_2\right)-\mu \left[1-\sqrt{1+\frac{{\left({\omega }_1-{\omega }_2\right)}^2}{4{\mu }^2}}\right]$	$\frac{1}{{𝒩}_{1,-}}\left({\Lambda }_{1,-}|{\nu }_2,{\nu }_1⟩+{\Lambda }_{2,-}|{\nu }_1,{\nu }_2⟩\right)$
$\kappa =2$	$2\mu +\frac{1}{2}\left({\omega }_1+{\omega }_2\right)$	$|{\nu }_2,{\nu }_2⟩$

Equations

$$H=\sum _{\mathbf{p}}{\omega }_{\mathbf{p}}\vec{B}\cdot {\vec{J}}_{\mathbf{p}}+\frac{\sqrt{2}{G}_F}{V}\sum _{\mathbf{p},\mathbf{q}}\left(1-\cos {\vartheta }_{\mathbf{pq}}\right){\vec{J}}_{\mathbf{p}}\cdot {\vec{J}}_{\mathbf{q}},$$

$$J_{\mathbf{p}}^{+}=a_1^{†}(\mathbf{p})a_2(\mathbf{p}),$$

$$J_{\mathbf{p}}^{-}=a_2^{†}(\mathbf{p})a_1(\mathbf{p}),$$

$$J_{\mathbf{p}}^z=\frac{1}{2}\left(a_1^{†}(\mathbf{p})a_1(\mathbf{p})-a_2^{†}(\mathbf{p})a_2(\mathbf{p})\right).$$

$$a_e(\mathbf{p})$$

$$a_x(\mathbf{p})$$

$$[J_{\mathbf{p}}^{+},J_{\mathbf{q}}^{-}]=2{\delta }_{\mathbf{pq}}{J}_{\mathbf{p}}^z,[J_{\mathbf{p}}^z,J_{\mathbf{q}}^{\pm }]=\pm {\delta }_{\mathbf{pq}}{J}_{\mathbf{p}}^{\pm }.$$

$${\vec{J}}_{\omega }=\sum _{|\mathbf{p}|=\frac{\delta m^2}{2\omega }}{\vec{J}}_{\mathbf{p}}.$$

$$\begin{array}{rl}{H}_{\nu \nu }& \approx \frac{\sqrt{2}{G}_F}{V}⟨1-\cos {\vartheta }_{\mathbf{pq}}⟩\vec{J}\cdot \vec{J}\\ & \equiv \mu (r)\vec{J}\cdot \vec{J},\end{array}$$

$$\mu (r)\propto \frac{1}{2}{\left(1-\sqrt{1-\frac{R_{\nu }^2}{r^2}}\right)}^2,$$

$$H=\sum _p{\omega }_p\vec{B}\cdot {\vec{J}}_p+\mu (r)\vec{J}\cdot \vec{J},$$

$$|j,+j⟩$$

$$|j,-j⟩$$

$$E_{\pm N/2}=\mp \sum _p{\omega }_p\frac{N_p}{2}+\mu \frac{N}{2}\left(\frac{N}{2}+1\right).$$

$$\vec{S}({\zeta }_{\alpha })\equiv \sum _p\frac{{\vec{J}}_p}{{\omega }_p-{\zeta }_{\alpha }}$$

$$\ket{\zeta_{1},\ldots,\zeta_{\kappa}}=\mathcal{N}(\zeta_{1},\ldots,\zeta_{\kappa})\bigg{(}\prod_{\alpha=1}^{\kappa}S^{-}_{\alpha}\bigg{)}\ket{j,+j}$$

$$J^z|{\zeta }_1,\dots ,{\zeta }_{\kappa }⟩=\left(j-\kappa \right)|{\zeta }_1,\dots ,{\zeta }_{\kappa }⟩.$$

$$-\frac{1}{2\mu }-\sum _{p=1}^M\frac{j_p}{{\omega }_p-{\zeta }_{\alpha }}=\sum _{\begin{array}{c}\beta =1\\ \beta \ne \alpha \end{array}}^{\kappa }\frac{1}{{\zeta }_{\alpha }-{\zeta }_{\beta }},$$

$$E=E_{+N/2}+\sum _{\alpha =1}^{\kappa }{\zeta }_{\alpha }-\mu \kappa (N-\kappa +1).$$

$$\Lambda (\lambda )=\sum _{\alpha =1}^{\kappa }\frac{1}{\lambda -{\zeta }_{\alpha }},$$

$$\Lambda (\lambda {)}^2+{\Lambda }^{\prime }(\lambda )+\frac{1}{\mu }\Lambda (\lambda )=\sum _{q=1}^{M}2j_q\frac{\Lambda (\lambda )-\Lambda ({\omega }_q)}{\lambda -{\omega }_q},$$

$${\Lambda }_p^2+(1-2j_p){\Lambda }_p^{\prime }+\frac{1}{\mu }{\Lambda }_p=\sum _{\begin{array}{c}q=1\\ q\ne p\end{array}}^{M}2j_q\frac{{\Lambda }_p-{\Lambda }_q}{{\omega }_p-{\omega }_q},$$

$${\Lambda }_p^2+\frac{1}{\mu }{\Lambda }_p=\sum _{\begin{array}{c}q=1\\ q\ne p\end{array}}^M\frac{{\Lambda }_p-{\Lambda }_q}{{\omega }_p-{\omega }_q},$$

$$\sum _{p=1}^M{j}_p{\omega }_p{\Lambda }_p=\kappa \sum _{p=1}^M{j}_p-\frac{1}{2\mu }\sum _{\alpha =1}^{\kappa }{\zeta }_{\alpha }-\frac{\kappa (\kappa -1)}{2},$$

$$E=E_{N/2}-2\mu \sum _{p=1}^M{j}_p{\omega }_p{\Lambda }_p.$$

$$\begin{array}{rl}\sum _{p=1}^M{j}_p{\omega }_p^k{\Lambda }_p& =\sum _{\alpha =1}^{\kappa }\sum _{p=1}^M{j}_p\frac{{\omega }_p^k}{{\omega }_p-{\zeta }_{\alpha }}\\ & =\sum _{\alpha }\sum _p{j}_p\frac{{\omega }_p^k-{\zeta }_{\alpha }^k}{{\omega }_p-{\zeta }_{\alpha }}+\sum _{\alpha }{\zeta }_{\alpha }^k\sum _p\frac{j_p}{{\omega }_p-{\zeta }_{\alpha }}\\ & =\sum _{\alpha }\sum _p{j}_p\left(\sum _{l=1}^k{\omega }_p^{k-l}{\zeta }_{\alpha }^{l-1}\right)\\ & -\sum _{\alpha }{\zeta }_{\alpha }^k\sum _{\beta \ne \alpha }\frac{1}{{\zeta }_{\alpha }-{\zeta }_{\beta }}-\sum _{\alpha }\frac{{\zeta }_{\alpha }^k}{2\mu }.\end{array}$$

$$\begin{array}{rl}\sum _{p=1}^M{j}_p{\omega }_p^k{\Lambda }_p& =\sum _{l=1}^k\left(\sum _{\alpha =1}^{\kappa }{\zeta }_{\alpha }^{l-1}\right)\left(\sum _{p=1}^M{j}_p{\omega }_p^{k-l}\right)\\ & -\frac{1}{2\mu }\sum _{\alpha =1}^{\kappa }{\zeta }_{\alpha }^k-\frac{1}{2}\sum _{\begin{array}{c}\alpha ,\beta =1\\ \alpha \ne \beta \end{array}}^{\kappa }\frac{{\zeta }_{\alpha }^k-{\zeta }_{\beta }^k}{{\zeta }_{\alpha }-{\zeta }_{\beta }}\\ & =\sum _{l=1}^k\left(\sum _{\alpha }{\zeta }_{\alpha }^{l-1}\right)\left(\sum _p{j}_p{\omega }_p^{k-l}\right)\\ & -\frac{1}{2\mu }\sum _{\alpha }{\zeta }_{\alpha }^k-\frac{1}{2}\sum _{\begin{array}{c}\alpha ,\beta \\ \alpha \ne \beta \end{array}}\sum _{l=1}^k{\zeta }_{\alpha }^{k-l}{\zeta }_{\beta }^{l-1}.\end{array}$$

$$\begin{array}{r}\sum _{p=1}^M{j}_p{\omega }_p^k{\Lambda }_p=-\frac{1}{2\mu }{Q}_k+\sum _{l=1}^k{Q}_{l-1}\left(\sum _{p=1}^M{j}_p{\omega }_p^{k-l}\right)\\ -\frac{1}{2}\left[\sum _{l=1}^k{Q}_{k-l}{Q}_{l-1}-k{Q}_{k-1}\right].\end{array}$$

$$\sum _p{j}_p{\Lambda }_p$$

$$\sum _p{j}_p{\omega }_p{\Lambda }_p$$

$$\sum _p{j}_p{\omega }_p^2{\Lambda }_p$$

$$-(\kappa -1)Q_1,$$

$$k{e}_k({\zeta }_1,\dots ,{\zeta }_{\kappa })=\sum _{i=1}^k(-1{)}^{i-1}{e}_{k-i}{Q}_i.$$

$$e_1$$

Full text

I Introduction

Experimental evidence has ascertained that neutrinos can undergo flavor oscillations, as a result of mass-differences between propagation eigenstates, which are distinct from the eigenstates of the weak interaction Fukuda et al. (1998); Ahmad et al. (2001); An et al. (2012). Additionally, it has been shown that the mixing between the two sets of eigenstates can be modified in the presence of coherent forward scattering of neutrinos on charged leptons Wolfenstein (1978); Mikheyev and Smirnov (1985); Haxton (1986); Bethe (1986), as well as coherent neutrino-neutrino forward scattering Fuller et al. (1987); Nötzold and Raffelt (1988); Pantaleone (1992); Sigl and Raffelt (1993); Raffelt et al. (1993); Samuel (1993). The latter contribution is of particular interest on account of its nonlinear nature, and becomes relevant in environments with high neutrino fluxes, such as the hot and dense early universe 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), as well as in compact-object systems, such as the neutrino emission accompanying a core-collapse supernova explosion or black hole-neutron star or binary neutron star merger. The interplay between the linear and nonlinear terms can result in various forms of interesting collective flavor oscillation phenomena Duan and Kneller (2009); Duan et al. (2010); Chakraborty et al. (2016); Balantekin (2018), including spectral splits/swaps 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), matter-neutrino resonances 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), and fast flavor oscillations arising from spatial or temporal instabilities 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. (2019); Abbar and Volpe (2019); Shalgar and Tamborra (2019).

Collective oscillations of neutrinos have been investigated in literature predominantly using the “mean field” approach, wherein each neutrino is considered to be interacting with a background mean field composed of all other neutrinos. However, it has been pointed out that this problem also lends itself to a full many-body description Bell et al. (2003); Friedland and Lunardini (2003a, b); Balantekin and Pehlivan (2007); Pehlivan et al. (2011); Birol et al. (2018); Cervia et al. (2019); Balantekin (2018). Such a description is much more complete than the mean field approach, in the sense that it operates within a larger Hilbert space, and captures exclusive many-body effects such as the formation and the evolution of entangled neutrino states. The many-body Hamiltonian describing neutrino oscillations exhibits a mathematical analogy to that of a one-dimensional spin chain with one-particle and “long-range” two-particle interactions (in momentum space), as well as the separable pairing Hamiltonians describing nucleon pairs present in different Shell Model orbitals Pan et al. (1998); Balantekin et al. (2007). Eigenvalues and eigenvectors of the latter class of Hamiltonians were already constructed in the 1960s by Richardson Richardson (1966). This solution was cast into an algebraic form by Gaudin Gaudin (1976), whose formalism we shall use extensively in this article.

The article is organized as follows. In Sec. II, we present the many-body neutrino oscillation Hamiltonian with vacuum and self-interaction terms, and simplify it using the single-angle approximation. Sec. III introduces the Bethe Ansatz method, also known as the Richardson-Gaudin diagonalization technique, and demonstrate how the eigenvalues and eigenvectors of the Hamiltonian may be systematically expressed using this method. In Sec. IV, we re-cast the Bethe Ansatz equations into a different set of simpler, more tractable equations, and also outline a procedure for calculating the eigenvalues and eigenstates of the Hamiltonian in terms of the solutions to these simpler equations. In Sec. V, we describe the analytic and numerical solutions for specific cases. We conclude in Sec. VI. Many of the mathematical details and derivations are presented in the Appendices.

II The many-body neutrino Hamiltonian

In scenarios where the neutrino scattering rates are sufficiently low for the neutrinos to be essentially free-streaming across the relevant physical scales, the flavor evolution of neutrinos is dominated by coherent forward-scattering processes. In that case, the interacting neutrinos may be described as a many-body Hamiltonian system. Generally speaking, such a Hamiltonian shall consist of terms that represent neutrino oscillations in vacuum, as well as interactions of neutrinos with ordinary background matter and with other neutrinos, along with the corresponding terms for anti-neutrinos.

In this article, for ease of discussion and to reduce numerical complexity, we make a number of simplifying assumptions. First and foremost, we restrict the discussion to two flavor/mass states of neutrinos, rather than the full three-flavor picture. Additionally, we pick a system consisting only of neutrinos, and no antineutrinos. Furthermore, keeping the environments where neutrino-neutrino interactions are dominant in mind, we ignore the interactions between neutrinos and ordinary (non-neutrino) background matter. With these assumptions, the Hamiltonian may be written as a sum of vacuum and self-interaction terms, $H_{\text{vac}}$ and $H_{\nu\nu}$, given by

[TABLE]

where $\vec{B}=(0,0,-1)_{\mbox{\tiny mass}}=(\sin 2\theta,0,-\cos 2\theta)_{\mbox{\tiny flavor}}$ is a unit vector indicating the direction of the mass basis in isospin space, with $\theta$ being the vacuum mixing angle, and $\omega_{\mathbf{p}}={\delta m^{2}}/{(2|\mathbf{p}|)}$ are the vacuum oscillation frequencies. Here, $\vartheta_{\mathbf{p}\mathbf{q}}$ is the intersection angle between the trajectories of neutrinos with 3-momenta $\mathbf{p}$ and $\mathbf{q}$, $V$ is the quantization volume, and $G_{F}$ is the Fermi coupling constant. Here, we have defined the neutrino mass-basis isospin operators $\vec{J}_{\mathbf{p}}$ in terms of the Fermionic creation and annihilation operators Balantekin and Pehlivan (2007)

[TABLE]

An analogous set of weak isospin operators may be defined in the flavor basis, using the corresponding creation/annihilation operators, which are related to their mass basis counterparts via the following unitary transformation

[TABLE]

The isospin operators form an SU(2) algebra, obeying the usual commutation relations

[TABLE]

At this point, it is important to note that the strength of the interaction between any two neutrinos is dependent on the intersection angle between their momenta, as can seen from the second term of the Hamiltonian in Eq. (1). This geometric dependence makes the collective neutrino oscillation problem extremely complex, even in the mean field limit. To avoid these complexities, and to facilitate a qualitative understanding of various collective flavor evolution phenomena, the so-called “single-angle” approximation has been frequently adopted in the literature, wherein the angle-dependent interaction is instead replaced by an overall, angle-averaged coupling strength. The implication of this approximation is that it removes any trajectory-dependence from the flavor evolution of neutrinos. It has been demonstrated that the single-angle calculations are able to capture many of the qualitative behaviors observed in the more sophisticated multi-angle calculations. To this end, we define a direction-independent weak isospin operator $\vec{J}_{\omega}$ as follows:

[TABLE]

The self-interaction term may then be approximated as

[TABLE]

where we have defined $\vec{J}=\sum_{\omega}J_{\omega}=\sum_{\mathbf{p}}J_{\mathbf{p}}$. The radial dependence of the coupling strength $\mu$ arises due to the geometric dilution of the neutrino fluxes and the narrowing of the intersection angles as one moves further from the source. For example in a “neutrino bulb” model, where neutrinos are assumed to be emitted isotropically from a single spherical emission surface, the dependence is given by

[TABLE]

where $R_{\nu}$ is the radius of the neutrinosphere.

III Diagonalization and the Bethe Ansatz equations

In the single-angle picture, the Hamiltonian from Eq. (1) may be re-written as

[TABLE]

where $p$ is now just an index denoting the oscillation frequencies present in the system. It has been shown that this particular Hamiltonian is amenable to diagonalization using the Richardson-Gaudin procedure. To begin with, one may observe that, among the common eigenstates $\ket{j,m}$ of the total weak isospin operators $\vec{J}^{2}$ and $J^{z}$, some can be written as direct products of eigenstates $\ket{j_{p},\pm j_{p}}$ of $\vec{J}_{p}^{2}$ and $J_{p}^{z}$. For instance, the states

[TABLE]

where $j=\sum_{p}j_{p}$, can be shown to be simultaneous eigenstates of $\vec{J}^{2}$ and each $J_{p}^{z}$, and are therefore eigenstates of the Hamiltonian111Note that, for $j<N/2$, not every state of the form $\ket{j,\pm j}$ is an eigenstate. $j=\sum_{p}j_{p}$ is a necessary condition.. In particular, the highest and lowest weight states $\ket{\frac{N}{2},\frac{N}{2}}$ and $\ket{\frac{N}{2},-\frac{N}{2}}$, where $N$ is the total number of neutrinos in the system, are eigenstates. These represent the states where all the neutrinos are either isospin-up or isospin-down, i.e., $\ket{\nu_{1},\ldots,\nu_{1}}$ and $\ket{\nu_{2},\ldots,\nu_{2}}$, respectively. The corresponding eigenvalues for these two states can be easily shown to be

[TABLE]

where $N_{p}$ is the number of neutrinos at the oscillation frequency $\omega_{p}$, and where we have suppressed the radial dependence in the notation of $\mu$ for convenience. The remaining eigenstates of the Hamiltonian may then be systematically constructed by first defining the Gaudin algebra with the operators

[TABLE]

where we have introduced a sequence of Bethe Ansatz variables $\zeta_{\alpha}$, which are yet to be determined, following the formalism introduced in Pehlivan et al. (2011).

Our Bethe Ansatz is the claim that the eigenstates of our system are states of the form

[TABLE]

where $\ket{j,+j}$ is an eigenstate, $S_{\alpha}^{-}=S^{-}(\zeta_{\alpha})$, and $\mathcal{N}(\zeta_{1},\ldots,\zeta_{\kappa})$ is a normalization factor. An eigenstate of the Hamiltonian, such as the one defined in Eq. (15), is typically a linear superposition of several mass basis states (i.e., Kronecker-products of mass states of individual neutrinos), each consisting of $N/2+j-\kappa$ neutrinos in the $\nu_{1}$ state and $N/2-j+\kappa$ neutrinos in the $\nu_{2}$ state. The fact that each eigenstate has a well-defined number of $\nu_{1}$ and $\nu_{2}$ is a consequence of the total $J^{z}=\sum_{p}J_{p}^{z}$ commuting with the Hamiltonian $H$ from Eq. (11). Eigenstates of the Hamiltonian are therefore also eigenstates of $J^{z}$:

[TABLE]

A similar procedure may be carried out in which we instead apply $S^{+}_{\alpha}$ operators to $\ket{j,-j}$ to obtain such eigenstates. In fact, these two procedures yield identical results. This fact can be shown by defining a rotation operator $T=e^{i\frac{\pi}{2}(J^{+}+J^{-})}$, for which $TJ^{z}T^{-1}=-J^{z}$ and $TJ^{\pm}T^{-1}=J^{\mp}$, and so $T\ket{j,-j}=\ket{j,+j}$. Moreover, one can see that, if $\ket{\psi}$ is an eigenstate obtained from Eq. (15), then $T\ket{\psi}$ is an eigenstate of a Hamiltonian $THT^{-1}=-\sum_{p}\omega_{p}\vec{B}\cdot\vec{J}_{p}+\mu(r)\vec{J}\cdot\vec{J}$ that may be obtained from a “raising formalism” using $S^{+}_{\alpha}$ instead.

For a particular eigenstate consisting of $\kappa$ neutrinos in the isospin-down ($\nu_{2}$) configuration, we will have $\kappa$ different ansatz variables $\zeta_{\alpha}$ to determine, and one can show Birol et al. (2018); Balantekin (2018) using the commutation relations between the Gaudin operators that this requirement is equivalent to the condition

[TABLE]

for $\alpha=1,\ldots,\kappa$. Here, $M$ is the total number of energy (or equivalently, $\omega$) values in the system, and where $j_{p}$ is related to the eigenvalue of the Casimir operator $\vec{J}_{p}^{2}$ (each $j_{p}$ can take values from [math] or $1/2$ to $N_{p}/2$). For a particular choice of $\kappa$, and $j_{p}$s, Eq. (17) can admit multiple solutions of the form $\{\zeta_{1},\ldots,\zeta_{\kappa}\}$. Obtaining the complete set of eigenstates and eigenvalues for an $N$-neutrino system therefore involves solving multiple sets of equations of the form Eq. (17), for $\kappa=1,\ldots,N$ and each $j_{p}=0\text{ or }1/2,\ldots,N_{p}/2$. For the specific case where $j_{p}=1/2$ for all $p$, the number of solutions of Eq. (17) for a given $\kappa$ is equal to ${}^{N}C_{\kappa}$, and the total number of solutions is $2^{N}$, across all values of $\kappa$.

One can always choose $j=N/2$ as the starting point for the construction of eigenstates using the Bethe Ansatz method. With this choice, the energy eigenvalue for the eigenstate defined in Eq. (15) can be shown to be

[TABLE]

IV The Lambda method

The Bethe Ansatz equations, Eq. (17), constitute a set of coupled algebraic equations in $\kappa$ variables. In principle, numerical solutions to this set of equations may be sought—however, in their original form the equations are unwieldy. For starters, the variables $\zeta_{\alpha}$ admit complex values, and the equations contain singularities for certain values of parameters $\omega_{p}$ and $\mu$ where the different $\zeta_{\alpha}$ approach each other. Moreover, if each equation were to be converted into a coupled polynomial form (by cross-multiplying all the denominators), then the order of each polynomial would be $M+\kappa-2$. Therefore, it is worthwhile to explore the possibilities of recasting the Bethe Ansatz equations into a different, more tractable form.

In order to accomplish this, one may introduce certain auxiliary functions which depend on the Bethe Ansatz variables. For instance, following Refs. Faribault et al. (2011); Babelon and Talalaev (2007) we may define the function

[TABLE]

and transform Eq. (17) into a first-order ordinary differential equation

[TABLE]

where the prime represents derivative with respect to $\lambda$. Eq. (20) is not straightforward to integrate because of the presence of the parameters $\Lambda(\omega_{q})$, whose values are not known a priori, and are in fact dependent on the equation itself. These parameters can be determined by taking the limit of Eq. (20) as $\lambda\to\omega_{p}$ for each $p=1,\ldots,M$. Doing so yields the following system of equations:

[TABLE]

where $\Lambda_{p}=\Lambda(\omega_{p})$ and $\Lambda_{p}^{\prime}=\Lambda^{\prime}(\omega_{p})$, for $p=1,\ldots,M$. In particular, if $j_{p}=1/2$ for all $p$, then our equations for $\Lambda_{p}$ reduce simply to the form

[TABLE]

yielding a system of coupled algebraic equations of quadratic order in the parameters $\Lambda_{p}$. Physically, this represents the case where the system is composed of neutrinos that all have pairwise distinct momenta, allowing us to choose a discrete set of $\omega_{p}$ bins in which each bin includes exactly one neutrino. Details of the derivation of Eqs. (20) and (22) are given in Appendix A.

Eq. (22) is manifestly much simpler than the original Bethe Ansatz equations Eq. (17). Moreover, unlike Eq. (17) where each value of $\kappa=1,\ldots,N$ requires solving a separate set of Bethe Ansatz equations, Eq. (22) represents just a single set of equations that can be solved to yield all the solutions corresponding to different values of $\kappa$. In Sec. IV.1, we show that the following relation holds between the variables $\Lambda_{p}$ and $\zeta_{\alpha}$:

[TABLE]

which may be used to express our energy eigenvalues from Eq. (18) in terms of $\Lambda_{p}$ instead of the Bethe Ansatz variables $\zeta_{\alpha}$:

[TABLE]

Thus, we seek to determine the parameters $\Lambda_{1},\ldots,\Lambda_{M}$ from Eq. (22). Following that, one may use Eq. (19), which may then be inverted to obtain the Bethe Ansatz variables $\zeta_{\alpha}$. The $\zeta_{\alpha}$ may then be used to reconstruct the states and their energies that solve our model using Eqs. (14) and (18). Alternatively, we show in this paper that for the system in which $j_{p}=1/2$ for all $p$, one may instead directly reconstruct the eigenstates and their energies using the variables $\Lambda_{p}$, without involving the Bethe Ansatz variables. These two procedures are detailed in Secs. IV.1 and IV.2, respectively.

IV.1 Obtaining the Bethe Ansatz roots $\zeta_{\alpha}$ from the variables $\Lambda_{p}$

After solving Eq. (22) for $\Lambda_{p}$, it is possible to reduce the problem of obtaining the Bethe Ansatz variables $\zeta_{\alpha}$ to that of solving a single polynomial equation of order $\kappa$. This process involves two steps, the first being deriving a set of constraint relations between the $\Lambda_{p}$s and the power sums of $\zeta_{\alpha}$s. Using the definition Eq. (19), one has

[TABLE]

Here, in the final step, we have used the Bethe Ansatz equation Eq. (17) to replace the second inner sum. Using symmetry arguments and changing the order of the summations, the above equation becomes

[TABLE]

Here, it is useful to define the power sums of the Bethe Ansatz variables, $Q_{k}\equiv\sum_{\alpha=1}^{\kappa}\zeta_{\alpha}^{k}$. This allows us to write the above expression as

[TABLE]

For the first few values of $k$, the above equation takes the following forms:

[TABLE]

for $k=0,1,\text{ and }2$, respectively. In particular, the equation for $k=0$ may be treated as a constraint relation for the solutions $\Lambda_{p}$ obtained by solving Eq. (22), and can also be used to classify those solutions according to $\kappa$. Also note that the equation for $k=1$ is identical to Eq. (23), which can be used to express the energy eigenvalues in terms of the $\Lambda_{p}$, as shown earlier.

Using Eq. (27), one can successively calculate the power sums $Q_{k}$ to any desired value of $k$, once all the $\Lambda_{p}$ are known. And the first $\kappa$ power sums, $Q_{1},\ldots,Q_{\kappa}$, can be used to obtain all the roots $\zeta_{\alpha}$. This involves first calculating the elementary symmetric polynomials of the $\zeta_{\alpha}$s from the power sums. The elementary symmetric polynomials are $e_{0}=1$, $e_{1}=\sum_{\alpha}\zeta_{\alpha}$, $e_{2}=\sum_{{\alpha,\beta<\alpha}}\zeta_{\alpha}\zeta_{\beta}$, and so on. For $k\leq\kappa$, these can be calculated recursively from the power sums using Newton’s identities:

[TABLE]

Therefore, one has

[TABLE]

and so on. Once the $\kappa$ elementary symmetric polynomials $e_{1},\ldots,e_{\kappa}$ are evaluated, then the polynomial $P(\lambda)\equiv\prod_{\alpha=1}^{\kappa}(\lambda-\zeta_{\alpha})$, whose roots are ${\zeta_{\alpha}:\alpha=1,\ldots,\kappa}$ may be constructed as

[TABLE]

Any one-dimensional polynomial root-finding algorithm can be employed to numerically obtain the roots $\zeta_{\alpha}$ of this polynomial. Once the $\zeta_{\alpha}$ are determined, then the eigenstates of the Hamiltonian may be explicitly constructed using Eq. (15). As an aside, it is also interesting to note that $P(\lambda)$ is related to the function $\Lambda(\lambda)$ as $\Lambda(\lambda)={P^{\prime}(\lambda)}/{P(\lambda)}$.

IV.2 Constructing the eigenstates directly using $\Lambda_{p}$

Alternatively, we show that it is possible to directly compute the eigenstates of the Hamiltonian in terms of the auxiliary variables, without needing to calculate the Bethe Ansatz variables first. Unlike the procedure described in Sec. IV.1, this method does not involve any additional numerical root-finding after obtaining the $\Lambda_{p}$, and therefore eliminates a potential source of numerical error. We can use the fact $[J_{p}^{-},J_{q}^{-}]=0$ to find that the right hand side of Eq. (15) can be rewritten directly in terms of $\Lambda_{p}$ (without involving the $\zeta_{\alpha}$) for any $\kappa$. Explicit derivations of these identities for $\kappa=2,3$ are provided in Appendices C and D, and here we present the argument for a generic $\kappa$. Let us define the $\kappa\times\kappa$ matrix $A$ with the matrix elements $A_{ij}=S_{i}^{-}\delta_{ij}$, where $\delta_{ij}$ is the Kronecker delta. Disregarding the normalization for the time being, one can rewrite Eq. (15) as

[TABLE]

where $e_{\kappa}$ is the $\kappa$-th elementary symmetric polynomial and where we have used the fact that $e_{\kappa}(x_{1},\ldots,x_{\kappa})=\prod_{j=1}^{\kappa}x_{j}$ for any variables $x_{j}$. Since $A$ is a square matrix over a commutative ring, we can use the Cayley-Hamilton Theorem to infer that

[TABLE]

where $S(\kappa)$ is the symmetry group of $\kappa$ letters and $\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})$ Bourbaki (1998). The traces may be expressed in terms of the $\Lambda_{p}$s via the following steps:

[TABLE]

Here we have used the Heaviside cover-up rule between the second and third equalities, and the definition from Eq. (19) between the third and fourth equalities (after exchanging the order of summation). An alternative derivation for the individual overlaps of energy eigenstates with the mass basis states is given in Ref. Claeys et al. (2017); Claeys (2018).

Since the form of $\mathrm{sgn}(\sigma)\mathrm{tr}_{\sigma}(A)$ depends on the cycle type—which have multiplicities $c_{\sigma}$—but not the particular $\sigma$, we can reduce Eq. (36) to the following:

[TABLE]

Lastly, we note that reversing the order of $l$ and $m$ in the denominator of Eq. (37) for each factor $\mathrm{tr}(A^{f})$ in $\mathrm{tr}_{\sigma}(A)$ produces a factor of $\mathrm{sgn}(\sigma)$, which cancels with the same factor already written in each term of Eq. (36). In fact, using a decomposition into cycle-types, $\sigma=(f_{1},\ldots,f_{n})$, we see that $\mathrm{sgn}(\sigma)=\prod_{i=1}^{n}(-1)^{f_{i}+1}=(-1)^{\kappa+n}$, while a factor of $\prod_{i=1}^{n}(-1)^{f_{i}-1}=(-1)^{\kappa-n}$ is produced from reversing the differences mentioned in the previous sentence.222Interestingly, the resulting form of Eq. (36) with this cancellation suggests that the module spanned by $S_{1}^{-},\ldots,S_{\kappa}^{-}$ forms a $\kappa$-th-power symmetric (as opposed to exterior) algebra, whose character can be described with the same formula.

Alternatively, since $\mathrm{tr}(A^{f})$ is simply the power sum $\sum_{i=1}^{\kappa}(S_{i}^{-})^{f}$, one may use Newton’s identities to systematically construct the elementary symmetric polynomial $e_{\kappa}(S_{1}^{-},\ldots,S_{\kappa}^{-})$ using the traces $\mathrm{tr}(A),\ldots,\mathrm{tr}(A^{\kappa})$. This was the basis of our numerical approach to calculating the eigenstates, the results of which are shown in Sec. V.

V Solutions of the Bethe Ansatz equations

Having obtained a set of algebraic equations in the auxiliary variables $\Lambda_{p}$, we are now in a position to discuss the solutions to these equations. To begin with, one can examine the solutions in the limit $\mu\to 0$, which we henceforth shall also refer to as the “asymptotic limit.” At this point, it is instructive to define the quantities $\widetilde{\Lambda}_{p}\equiv\mu\Lambda_{p}$ to rewrite Eq. (22) as

[TABLE]

For convenience, we have defined $\eta_{pq}=\omega_{p}-\omega_{q}$. Taking the limit $\mu\to 0$ decouples the various modes from each other, yielding the solutions $\widetilde{\Lambda}_{p}=0$ or $-1$, independently for each $p$. These asymptotic solutions can serve as a starting point for efficiently calculating numerical solutions at a generic $\mu>0$. This also makes it easy to see that for the case where $j_{p}=1/2$ for all $p$, $2^{M}$ solutions exist in total. In fact, one can infer that the number of solutions is $2^{M}$ even when $\mu>0$, by noting that Eq. (39) is a set of $M$ mutually co-prime coupled quadratic equations in $M$ variables, and therefore has a finite solution set. One can then invoke Homotopy continuation to argue that each solution for a generic $\mu>0$ is continuously connected to a unique solution in the $\mu\to 0$ limit Garcia and Zangwill (1979) (see also, Appendix B).

V.1 Algebraic Solutions for $M=2$

For the purposes of gaining some mathematical and physical intuition, we first show the analytic solutions for a simple system consisting only of two interacting neutrinos, at frequencies $\omega_{1}$ and $\omega_{2}$. This corresponds to solving Eq. (22) for $M=2$ and can serve as a test bed for the numerical technique that we later implement for dealing with the $M>2$ cases, for which algebraic solutions either do not exist or are difficult to obtain.

We thus aim to determine each $\Lambda_{p}$ for the model with $M=N=2$ with $j_{1}=j_{2}=1/2$. These $\Lambda_{p}$ can easily be found analytically from Eq. (22). The four solutions are listed in Table 1, and categorized by their corresponding $\kappa$ values. For a particular solution $\{\Lambda_{p}:p=1,\ldots,M\}$, $\kappa$ can be determined using the identity

[TABLE]

which was derived in Sec. IV.1 (Eq. (28)). Note that, across a complete set of solutions, $\kappa$ takes all values from $0,1,\ldots,N$, and Eq. (40) holds for arbitrary $M=N$.

Furthermore, using the algebraic $\Lambda_{p}$ solutions from Table 1, we may compute the energy eigenvalues using Eqs. (13) and (24) and eigenstates using Eqs. (36) and (37) (or more specifically, Eq. (62) for $M=2$). These eigenvalues and eigenstates are aggregated in Table 2. Note that the $\Lambda_{1,\pm}$ and $\Lambda_{2,\pm}$ referenced in the latter table for $\kappa=1$ are respectively from the two solutions in Table 1 for $\kappa=1$. Additionally, normalization coefficients calculated for states with $\kappa=1,2$ are given by

[TABLE]

respectively. Note that even though these normalization coefficients are singular as $\mu\rightarrow 0$, the eigenstates themselves are not.

V.2 Numerical Solutions for $M>2$

Extending from our argument of the previous section that demonstrated there must be finitely many solutions to Eq. (22) (or, equivalently, to Eq. (39)), we can further observe that the problem with $M>2$ essentially necessitates solving a polynomial of degree higher than $4$ in one variable. From this observation we expect that we may not be able to algebraically solve Eq. (22) or (39). Nevertheless, we may solve these equations numerically, as we will outline in this section.

Previously, we observed that we may obtain simple limits to the solutions of Eq. (39) for any $M\geq 2$ as $\mu\to 0$, which we refer to as the “asymptotic solutions” to our system. Moreover, we can show that the algebraic system described by Eq. (39) is, in a sense, homotopic (continuously deformable) to the “asymptotic system” in which we take the limit $\mu\to 0$ prior to determining its solutions. By applying a homotopy method for numerically solving algebraic systems of equations as studied in Ref. Garcia and Zangwill (1979), we may determine the solutions to Eq. (39) by tracking our solutions as they vary from the $2^{M}$ solutions in the limit $\mu\to 0$ of our system to our desired final $\mu>0$. In fact, we build upon this existing homotopy method for our purposes: if $0\leq t\leq 1$ is our homotopy parameter, then we observe that not only our solutions obtained for $t=1$, but also all solutions determined for intermediate values $0<t<1$ are physical solutions to Eq. (39). A proof of the applicability of the homotopy method to our system is given in Appendix B.

Numerical solutions can then be obtained using a generalization of the Newton-Raphson method by the homotopy method. We can modify this method to accommodate our homotopy process in estimating solutions to Eq. (39) for generic $\mu$. We apply Newton’s method to solve our system of equations at the $n^{\text{th}}$ homotopy step:

[TABLE]

for each $t^{(n)}>0$, where $\vec{\mathcal{F}}$ is defined in Eq. (57). At the $n$-th Homotopy step, the initial guess for Newton’s method is taken to be the numerical solution $\widetilde{\Lambda}_{p}^{(n-1)}$ obtained by solving the system of equations at the $(n-1)^{\text{th}}$ step, i.e.,

[TABLE]

Given the form of $\vec{\mathcal{F}}$ in Eq. (57), we can easily compute the analytic form of the Jacobian $\partial_{\vec{z}}\vec{\mathcal{F}}(\vec{z}_{0},t)$, to used in each step of Newton’s method. Finally, all $2^{M}$ solutions to our system can be obtained by repeating the procedure, starting from each $\widetilde{\Lambda}_{p}$ independently taking values [math] or $-1$ at $\mu=0$.

To test this method, we computed the numerical solutions to Eq. (22) for $M=2$, and compared them with the analytic results described in Sec. V.1. The root-mean-squared relative error between the analytic and numerical values in comparing each component of each solution was $\lesssim 10^{-14}$, suggesting that the method was able to attain a high level of numerical accuracy.

A potential issue with the utilization of homotopy continuation with large-$M$-dimensional algebraic systems is that the close proximity of solutions in $\mathbb{R}^{M}$ may cause numerical solvers to jump between distinct solutions as $t$ varies between 0 and 1. For larger $M$ values (such as $M=10$), introducing error-correcting algorithms during each Newton-Raphson step helps in obtaining the expected continuity of solutions as $\mu$ is varied.

We applied this method to systems of interacting neutrinos with equally spaced oscillation frequencies given by when $\omega_{p}=p\omega_{0}$, just that only a single neutrino resides at each oscillation frequency (therefore, $j_{p}=1/2$ for all $p$ and so $M=N$). Such a system of $N$ neutrinos admits $2^{N}$ solutions. Solutions of the system for up to $N\sim 15$ can be computed within a reasonable time frame using a personal computer. In Fig. 1, we show the evolution with $\mu$ of one of the solutions to such a system with $N=10$ neutrinos. Shown in the figure are the quantities $\Lambda_{p}=\widetilde{\Lambda}_{p}/\mu$, for $p=1,\ldots,10$.

We also used these solutions to calculate the eigenvalues of our Hamiltonian using Eq. (11). In Fig. 2 we show all the energy eigenvalues corresponding to solutions with $\kappa=0,\ldots,5$, for the same system as in Fig. 1. We also note that results for $\kappa>5$ may be obtained more efficiently by implementing a “raising formalism” analogous to our procedure from Eqs. (15)–(24), in which we apply $S^{+}_{\alpha}$ operators from our Gaudin algebra to the state $\ket{j,-j}$. In analogy to $\widetilde{\Lambda}_{p}$, one can define the variables $\widetilde{\Lambda}_{p}^{(\uparrow)}$ in the raising formalism, which obey the coupled quadratic equations

[TABLE]

Similar to Eq. (40), the solutions $\widetilde{\Lambda}_{p}^{(\uparrow)}$ obey the constraint relations

[TABLE]

where ${\Lambda_{p}^{(\uparrow)}}=\widetilde{\Lambda}_{p}^{(\uparrow)}/\mu$. In particular, a solution $\widetilde{\Lambda}^{(\uparrow)}_{p}$ of Eq. (45) can be shown to correspond to a solution $\widetilde{\Lambda}^{(\downarrow)}_{p}$ of Eq. (39), via the identity

[TABLE]

Using Eq. (40), it can be shown that the corresponding $\kappa^{(\uparrow)}$ of a solution in the raising formalism can be related to a given $\kappa^{(\downarrow)}$ via $\kappa^{(\uparrow)}=M-\kappa^{(\downarrow)}$, further reinforcing this duality. This correspondence of solutions for a given $\kappa$ and $M-\kappa$ is also discussed in Refs. Claeys et al. (2017); Birol et al. (2018); Pehlivan et al. (2011). In particular, Ref. Claeys et al. (2017) notes that, in the limit $\mu\to\infty$, $\Lambda_{p}^{(\uparrow)}\to\Lambda_{p}^{(\downarrow)}$.

For each value of $\kappa$, one can observe that the energy eigenvalues form $b+1$ distinct branches, where $b\equiv\min\{\kappa,M-\kappa\}$, in the limit of large $\mu$, as previously noticed in Birol et al. (2018). One may also observe that the self-interaction term in the Hamiltonian from Eq. (11) becomes dominant as $\mu\to\infty$, and therefore the eigenstates of $H$ approximately align with the eigenstates $\ket{j,m}$ of $\vec{J}\cdot\vec{J}$, with eigenvalues $E\approx\mu j(j+1)$. For an energy eigenstate $\ket{\psi}$ with a given $\kappa$, we can see that $J^{z}\ket{\psi}=m\ket{\psi}$ where $m=M/2-\kappa$; therefore as $\mu\to\infty$, the quantum number $j$ of all the eigenstates for that $\kappa$ may take values in the range $M/2,\ldots,M/2-b$. So, as $\mu\to\infty$, we expect that the energy eigenvalues of states with a given $\kappa$ will split into $b+1$ branches with the approximate energies given above.

In Fig. 3a, we zoom in on the small-$\mu$ regions of Fig. 2f (all eigenvalues of the $\kappa=5$ eigenstates) to better illustrate the various energy-level crossings that are present in this region. To examine the nature of these level crossings in better detail, we show in Fig. 3b the energy eigenvalues of all the states (spanning all permitted $\kappa$ values) in a system with $N=M=3$. Eigenstates corresponding to each $\kappa$ are assigned a particular color, as described in the figure caption. Ref. Birol et al. (2018) observes that the highest-energy states from among the solutions of each $\kappa$ (highlighted here as dashed lines) do not have energy level crossings with other states of the same $\kappa$, which is consistent with our observation. Here, we would additionally like to point out that the highest energy eigenstate of a particular $\kappa$ can cross with other eigenstates of a different $\kappa$. However, since $J^{z}$ is a symmetry of the Hamiltonian, and each $\kappa$ corresponds to a unique $J^{z}$ eigenvalue (see Eq. (16)), these crossings do not result in mixing between these eigenstates.

In addition to our energy eigenvalues, we may calculate energy eigenstates for generic $\mu$ using Eqs. (35)–(37). These eigenstates may be encoded as a sequence of their $2^{M}$ coefficients in the mass basis. Particular examples of eigenstates for systems with for $M=3,4$ are shown in Fig. 4. A complete knowledge of the eigenvalues and eigenstates of the system enables calculating the time-evolution (i.e., evolution with $\mu$) of an arbitrary initial state, in the adiabatic limit. This can be used for the purposes of comparing the results of a many-body calculation to the corresponding results in the mean-field limit, as well as for studying exclusive many-body effects like the emergence of quantum entanglement between the various modes as the system evolves.

VI Conclusions

Adiabatic evolution of a many-neutrino system in the single-angle approximation is exactly solvable in the sense that such an evolution can be completely characterized by the solutions of the appropriate Bethe ansatz equations. However, solving those nonlinear Bethe ansatz equations is a highly nontrivial problem. In this paper we presented a technique to evaluate the exact adiabatic eigenvalues and eigenstates of the collective neutrino oscillation Hamiltonian by casting Bethe ansatz equations first into a differential form, then into a set of algebraic equations which are numerically more tractable than the original Bethe ansatz equations. With our outlined procedure, to determine solutions for up to $N\sim 15$ for the interesting region of $\mu\leq 5$ a personal computer is sufficient. Going to higher values of $N$ would require more computational power. In the future we will explore the behavior of the solutions as $N$ gets larger.

An immediate benefit of obtaining such exact solutions is the ability to explore the limits of applicability of the commonly used mean-field solutions. Notably, our exact problem (with $j_{p}=1/2$) has $2^{N}$ solutions, while the mean-field problem has only $2N$ solutions in total. Clearly, in the mean-field case we are either losing or combining many states. Earlier studies of the question of entangled neutrinos also explored modifications of the collective oscillations due to the many-body effects present when one goes beyond the one-body description inherent to the mean-field approximation Bell et al. (2003); Friedland and Lunardini (2003a, b). These papers primarily investigate part of the Hamiltonian which survives in the large $\mu$ limit and explore if the oscillations speed up due to many-particle entanglements. The growth rates they obtain differ by a factor of $\sqrt{N}$ relative to each other, depending on the setup. In future publications we plan to explore many-body entanglement effects using our approach, for generic values of $\mu$ where both the one-body and two-body terms in the Hamiltonian play a role. Clearly exploring these issues is a step which needs to be taken before moving on to astrophysics applications.

Many of the elements heavier than iron were formed by rapidly capturing neutrons on seed nuclei (r-process nucleosynthesis). In the astrophysical site of the r-process nucleosynthesis many reactions take place during a rather short duration, a feature which is typically associated with explosive phenomena. Currently the leading candidates for the sites of r-process nucleosynthesis are core-collapse supernovae and binary neutron star mergers. It is known that these nucleosynthesis yields are sensitive to the neutrino flavor evolution inside the supernova envelope. Here we presented an exact procedure to describe neutrino propagation including neutrino-neutrino interactions. It would be also interesting to explore whether any possible differences between the mean-field approximation and the exact adiabatic many-body approach would impact nucleosynthesis yields as well as supernova neutrino detection in terrestrial experiments.

Acknowledgements.

We would like to thank E. Armstrong, S. Birol, S. Coppersmith, G. Fuller, E. Grohs, A. Hashimoto, C. Johnson, E. Rrapaj, J. Schmidt, M. Sen, and S. Shalgar for valuable conversations. This work was supported in part by the U.S. National Science Foundation Grants PHY-1630782 and PHY-1806368.

Appendix A Deriving the Lambda equations from the Bethe Ansatz

In this appendix, we shall present the derivation of the coupled quadratic equations in Lambda, Eq. (22), from the Bethe Ansatz equations Eq. (17). These results are dispersed through the condensed-matter physics literature. Here for the convenience of the reader we gather them in one place. Using the function $\Lambda(\lambda)$ defined in Eq. (19), one can write

[TABLE]

Using partial fraction decomposition and symmetry arguments, this may be rewritten as

[TABLE]

Replacing $W(\zeta_{\alpha})$ with the left hand side of the Bethe Ansatz equations Eq. (17), one obtains

[TABLE]

Using the definition of $\Lambda(\lambda)$ from Eq. (19) after changing the order of summation and using partial fraction decomposition, this equation reduces to

[TABLE]

It can be shown that the ordinary differential equation, Eq. (51), is exactly equivalent to the Bethe Ansatz equations, i.e., every solution of Eq. (51) corresponds to a unique solution of the Bethe Ansatz equations, and vice versa. This equivalence can be proven using the fact that every step of the above derivation is reversible—one could just as easily start with Eq. (51) and derive the Bethe Ansatz equations, Eq. (17).

The ordinary differential equation in Eq. (51) may be converted to a set of algebraic equations in $\Lambda(\omega_{p})$, for each $\omega_{p}$ in the system. In order to do so, one can Taylor-expand $\Lambda(\lambda)$ around $\lambda=\omega_{q}$. Following this Taylor expansion, the right-hand side of Eq. (51) may be written as

[TABLE]

where $\Lambda^{\prime}(\omega_{q}),\Lambda^{\prime\prime}(\omega_{q}),\ldots$ are the successive derivatives of $\Lambda(\lambda)$ with respect to $\lambda$, as evaluated at $\lambda=\omega_{q}$. Using this expansion and then taking the limit $\lambda\rightarrow\omega_{q}$, one obtains the following set of equations (one for each $\omega_{q}$)

[TABLE]

If $j_{q}=1/2$ for each $q$ (corresponding to there being just one neutrino per bin), then the $\Lambda^{\prime}$ term vanishes and the equations become purely algebraic (Eq. (22)), and can in principle be solved to obtain the solutions $\{\Lambda(\omega_{q}):q=1,\ldots,M\}$. If $j_{q}>1/2$, then one can take successive derivatives of Eq. (51) with respect to $\lambda$, as noted in Refs. Babelon and Talalaev (2007); Faribault et al. (2011); Araby et al. (2012); Claeys et al. (2015). For example, taking the first derivative of Eq. (51) followed by the $\lambda\rightarrow\omega_{q}$ limit gives us

[TABLE]

Now, if $j_{q}=1$, then the $\Lambda^{\prime\prime}$ term vanishes, and one may use Eqs. (53) and (54) to eliminate $\Lambda^{\prime}(\omega_{q})$ from the system of equations. Likewise, if $j_{q}=3/2$, one may obtain an additional equation by taking the second derivative of Eq. (51), and use it along with Eqs. (53) and (54) to eliminate $\Lambda^{\prime}(\omega_{q})$ and $\Lambda^{\prime\prime}(\omega_{q})$. And so on, for higher values of $j_{q}$.

Appendix B Proof of the applicability of the homotopy continuation method

To carry out this homotopy method, let us first rewrite Eq. (39) as the system $\vec{F}(\vec{z})=\vec{0}$ using the functions

[TABLE]

for $p=1,\ldots,M$. In particular, since we know each value of $\widetilde{\Lambda}_{p}$ must be real333Since the eigenvalues of a Hermitian operator are real (cf. Eqs. (23) and (24)), any complex-valued solutions $\zeta_{\alpha}$ of Eq. (17) must come in complex conjugate pairs Pehlivan et al. (2011); Birol et al. (2018), implying by Eq. (19) that $\Lambda(\omega_{p})$ and therefore $\widetilde{\Lambda}(\omega_{p})$ must be real. , we restrict $\vec{z}\in\mathbb{R}$. Additionally, we can define a set of functions $\vec{G}$ with simpler solutions by

[TABLE]

and the family of functions $\vec{\mathcal{F}}(\vec{z},t)=(1-t)\vec{G}(\vec{z})+t\vec{F}(\vec{z})$ for $0\leq t\leq 1$. Then, from the explicit form of

[TABLE]

we see that $\vec{\mathcal{F}}$ is an analytic (and therefore a $C^{2}$) mapping of $\mathbb{R}^{M}\times[0,1]\to\mathbb{R}^{M}$. By Sard’s Theorem, we find that $\vec{\mathcal{F}}$ is a regular function. Recall that regular, in the context of algebraic geometry, means that: (i) $\vec{\mathcal{F}}(\vec{z}_{0},t)=\vec{0}$ implies $\partial_{\vec{z}}\vec{\mathcal{F}}(\vec{z}_{0},t)$ has rank $M$ for fixed $t=0,1$ at almost all $\vec{z}_{0}\in\mathbb{R}^{M}$ and (ii) $\vec{\mathcal{F}}(\vec{z_{0}},t)=\vec{0}$ implies $\partial_{\vec{z},t}\vec{\mathcal{F}}(\vec{z_{0}},t)$ has rank $M$ for fixed $0<t<1$ at almost all $\vec{z}_{0}\in\mathbb{R}^{M}$. (To clarify, $\partial_{\vec{z},t}\vec{\mathcal{F}}$ is the total Jacobian of $\vec{H}$ while $\partial_{\vec{z}}\vec{\mathcal{F}}$ is the Jacobian of $\vec{\mathcal{F}}$ stripped of its final column, $\partial_{t}\vec{\mathcal{F}}$.) Further, a theorem by Garcia and Zangwill Garcia and Zangwill (1979) states that if, in addition, for all sequences with $||\vec{z}^{(n)}||\to\infty$ there is a $p$ such that $F_{p},G_{p}$ satisfy

[TABLE]

for some $\alpha>0$ and a sub-sequence of points $\vec{z}^{(n_{m})}$ at which $G_{p}(\vec{z}^{(n_{m})})\neq 0$, then all roots of $\vec{F}$ can be found one-to-one from the known roots of $\vec{G}$ (i.e., the asymptotic solutions as $\mu\to 0$) by homotopy curves as we vary $t=0$ to $t=1$. The symbol $\nless$ implies that the real part of the limit (which generically may be complex if our chosen $\vec{F}$ were not strictly real) must be greater than or equal to a positive $\alpha$ that we are free to choose. We can see that this additional hypothesis is true for our system by first calculating for an arbitrary $p$ that

[TABLE]

for any sequence of points $\vec{z}^{(n)}$ after eliminating any of the points in the sequence solving $G_{p}(\vec{z})=0$ for some $p$. (We know we may choose such a $p$ since we know there are finitely many points solving $G_{p}(\vec{z})=0$ for all $p$.) Moreover, we may observe that this sequence must have at least one $r$ for which $|z_{r}^{(n)}|\to\infty$ at the fastest rate amongst the $\vec{z}^{(n)}$. Selecting $p=r$ we find that $F_{p}(\vec{z}^{(n)})/G_{p}(\vec{z}^{(n)})\to 1$ on this sub-sequence. Thus, we have shown our desired system $\vec{F}$ is homotopic to $\vec{G}$ as we vary $t\mu$ from [math] to $\mu$, in the sense described by Ref. Garcia and Zangwill (1979), and so all our solutions to $\vec{F}(\vec{z})=\vec{0}$ can be found from those of $\vec{G}(\vec{z})=\vec{0}$ via a single homotopy parameter. Thus we justify our use of the algebraic homotopy method for finding our $\widetilde{\Lambda}(\omega_{p})$ solutions.

Additionally, we confirm from this process that there must be exactly $2^{M}$ solutions to Eq. (39) for $\mu\neq 0$ as well. Moreover, for a given $\kappa=-\sum_{p=1}^{M}\widetilde{\Lambda}_{p}$, there are ${}^{M}C_{\kappa}$ solutions. We can observe that there are equally as many ways for us to permute the values $0,-1$ amongst our different $\widetilde{\Lambda}_{p}$ for a given $\kappa$. Thus, we expect that by considering all possible permutations of $0,-1$ for each $\widetilde{\Lambda}_{p}$ at $\mu=0$ we may find all $\sum_{\kappa=0}^{M}{}^{M}C_{\kappa}=2^{M}$ solutions.

Lastly, consider the situation in which we fix the values of each $\omega_{p}$ while we allow $\mu\geq 0$ to increase. Observe from our above discussion of the homotopy method that the maximum value of $\mu$ (i.e., $t\mu$ with $t=1$) in Eq. (57) is simply an arbitrary real $\mu>0$. So, we expect that we may obtain not only the numerical solutions to $\vec{F}(\vec{z})=\vec{\mathcal{F}}(\vec{z},1)=\vec{0}$ but also the $n^{\mathrm{th}}$ ($n=0,1\ldots$) intermediate numerical solutions $\vec{z}^{(n)}$ to $\vec{\mathcal{F}}(\vec{z}^{(n)},t^{(n)})=\vec{0}$ for each $0\leq t^{(n)}\leq 1$, provided that our steps $t^{(n+1)}-t^{(n)}>0$ are sufficiently small for the discrete steps of our homotopy method to be appropriate. In fact, by applying this method to only $\vec{F}|_{\mu=\mu_{\max}}$, we may obtain the solutions to Eq. (39) for all $0\leq\mu\leq\mu_{\max}$ with an arbitrary maximum $\mu_{\max}$, as $\vec{z}^{(n)}$ may be interpreted as the solutions to $\vec{F}|_{\mu=t^{(n)}\mu_{\max}}$.

Appendix C Derivations of Eigenstates for $\kappa=2,3$

We can demonstrate the claim from Section  IV.2 that our general eigenstates can be written in terms of $\Lambda_{p}$ in lieu of $\zeta_{\alpha}$ for $\kappa=2,3$ when $j_{p}=1/2$ for all $p$. Here we present two ways to derive the eigenstate expressions for $\kappa=2,3$. First note that, for all $\kappa$, we have the identity:

[TABLE]

where $e_{\kappa}$ is the degree-$\kappa$ elementary symmetric polynomial in $M$ variables. Now, for $\kappa=2$ we have

[TABLE]

and so

[TABLE]

Thus, we have rewritten the eigenstates as well as the eigenvalues of the Hamiltonian in Eq. (11) completely in terms of our defined $\Lambda(\omega_{p})$ parameters, bypassing the need to directly compute the Bethe Ansatz variables $\zeta_{\alpha}$. This process can be carried out for $\kappa=3$ similarly:

[TABLE]

and

[TABLE]

where in the first line we are using $\beta\neq\gamma$ and $\beta,\gamma\neq\alpha$. Thus,

[TABLE]

There is a more direct approach to arrive at the same results, by decomposing terms of $S_{1}^{-}S_{2}^{-}$ and $S_{1}^{-}S_{2}^{-}S_{3}^{-}$ immediately without consideration of other products – as we will do in the next section. Also, note that the form of these results bears similarity to the characters of exterior powers of vector spaces (i.e., $\Lambda^{\kappa}V$ for $\kappa=2,3$), suggesting a possibility that knowledge of these characters may allow us to rewrite $\prod_{\alpha=1}^{\kappa}S_{\alpha}^{-}$ in terms of $\Lambda_{p}$ more swiftly.

Appendix D Alternative Derivation of Eigenstates for $\kappa=2,3$

Here, we present a more direct method of evaluating the product $\prod_{i=1}^{\kappa}S_{\alpha}^{-}\ket{j,+j}$ for the cases of $\kappa=2,3$. In $\kappa=2$, we can quickly compute:

[TABLE]

Analogously, we can extend this argument to $\kappa=3$:

[TABLE]

while

[TABLE]

[TABLE]

and

[TABLE]

Thus, we arrive at

[TABLE]