Eigenvalues and eigenstates of the many-body collective neutrino oscillation problem
EEigenvalues and eigenstates of the many-body collective neutrino oscillation problem
Amol V. Patwardhan,
1, 2, ∗ Michael J. Cervia, † and A. Baha Balantekin ‡ Department of Physics, University of California, Berkeley, California 94720-7300, USA Department of Physics, University of Wisconsin, Madison, Wisconsin 53706, USA (Dated: August 19, 2019)We demonstrate a method to systematically obtain eigenvalues and eigenstates of a many-bodyHamiltonian describing collective neutrino oscillations. The method is derived from the Richardson-Gaudin framework, which involves casting the eigenproblem as a set of coupled nonlinear “BetheAnsatz equations”, the solutions of which can then be used to parametrize the eigenvalues andeigenvectors. The specific approach outlined in this paper consists of defining auxiliary variablesthat are related to the Bethe-Ansatz parameters, thereby transforming the Bethe-Ansatz equationsinto a different set of equations that are numerically better behaved and more tractable. We showthat it is possible to express not only the eigenvalues, but also the eigenstates, directly in terms ofthese auxiliary variables without involving the Bethe Ansatz parameters themselves. In this paper,we limit ourselves to a two-flavor, single-angle neutrino system.
PACS numbers: 14.60.Pq, 97.60.Bw, 13.15.+g, 26.30.-k, 26.50.+x
I. INTRODUCTION
Experimental evidence has ascertained that neutri-nos can undergo flavor oscillations, as a result of mass-differences between propagation eigenstates, which aredistinct from the eigenstates of the weak interaction [1–3]. Additionally, it has been shown that the mixingbetween the two sets of eigenstates can be modified inthe presence of coherent forward scattering of neutrinoson charged leptons [4–7], as well as coherent neutrino-neutrino forward scattering [8–13]. The latter contribu-tion is of particular interest on account of its nonlinearnature, and becomes relevant in environments with highneutrino fluxes, such as the hot and dense early uni-verse [14–23], as well as in compact-object systems, suchas the neutrino emission accompanying a core-collapsesupernova explosion or black hole-neutron star or binaryneutron star merger. The interplay between the linearand nonlinear terms can result in various forms of in-teresting collective flavor oscillation phenomena [24–27],including spectral splits/swaps [28–50], matter-neutrinoresonances [23, 51–60], and fast flavor oscillations arisingfrom spatial or temporal instabilities [61–73].Collective oscillations of neutrinos have been investi-gated in literature predominantly using the “mean field”approach, wherein each neutrino is considered to be in-teracting with a background mean field composed of allother neutrinos. However, it has been pointed out thatthis problem also lends itself to a full many-body descrip-tion [27, 49, 74–79]. Such a description is much morecomplete than the mean field approach, in the sense thatit operates within a larger Hilbert space, and captures ex-clusive many-body effects such as the formation and theevolution of entangled neutrino states. The many-body ∗ [email protected] † [email protected] ‡ [email protected] Hamiltonian describing neutrino oscillations exhibits amathematical analogy to that of a one-dimensional spinchain with one-particle and “long-range” two-particle in-teractions (in momentum space), as well as the separablepairing Hamiltonians describing nucleon pairs present indifferent Shell Model orbitals [80, 81]. Eigenvalues andeigenvectors of the latter class of Hamiltonians were al-ready constructed in the 1960s by Richardson [82]. Thissolution was cast into an algebraic form by Gaudin [83],whose formalism we shall use extensively in this article.The article is organized as follows. In Sec. II, wepresent the many-body neutrino oscillation Hamiltonianwith vacuum and self-interaction terms, and simplify itusing the single-angle approximation. Sec. III introducesthe Bethe Ansatz method, also known as the Richardson-Gaudin diagonalization technique, and demonstrate howthe eigenvalues and eigenvectors of the Hamiltonianmay be systematically expressed using this method. InSec. IV, we re-cast the Bethe Ansatz equations into a dif-ferent set of simpler, more tractable equations, and alsooutline a procedure for calculating the eigenvalues andeigenstates of the Hamiltonian in terms of the solutionsto these simpler equations. In Sec. V, we describe theanalytic and numerical solutions for specific cases. Weconclude in Sec. VI. Many of the mathematical detailsand derivations are presented in the Appendices.
II. THE MANY-BODY NEUTRINOHAMILTONIAN
In scenarios where the neutrino scattering rates aresufficiently low for the neutrinos to be essentially free-streaming across the relevant physical scales, the flavorevolution of neutrinos is dominated by coherent forward-scattering processes. In that case, the interacting neu-trinos may be described as a many-body Hamiltoniansystem. Generally speaking, such a Hamiltonian shall a r X i v : . [ nu c l - t h ] A ug consist of terms that represent neutrino oscillations invacuum, as well as interactions of neutrinos with ordi-nary background matter and with other neutrinos, alongwith the corresponding terms for anti-neutrinos.In this article, for ease of discussion and to reduce nu-merical complexity, we make a number of simplifying as-sumptions. First and foremost, we restrict the discus-sion to two flavor/mass states of neutrinos, rather thanthe full three-flavor picture. Additionally, we pick a sys-tem consisting only of neutrinos, and no antineutrinos.Furthermore, keeping the environments where neutrino-neutrino interactions are dominant in mind, we ignorethe interactions between neutrinos and ordinary (non-neutrino) background matter. With these assumptions,the Hamiltonian may be written as a sum of vacuum andself-interaction terms, H vac and H νν , given by H = (cid:88) p ω p (cid:126)B · (cid:126)J p + √ G F V (cid:88) p , q (1 − cos ϑ pq ) (cid:126)J p · (cid:126)J q , (1)where (cid:126)B = (0 , , − mass = (sin 2 θ, , − cos 2 θ ) flavor is aunit vector indicating the direction of the mass basis inisospin space, with θ being the vacuum mixing angle, and ω p = δm / (2 | p | ) are the vacuum oscillation frequencies.Here, ϑ pq is the intersection angle between the trajec-tories of neutrinos with 3-momenta p and q , V is thequantization volume, and G F is the Fermi coupling con-stant. Here, we have defined the neutrino mass-basisisospin operators (cid:126)J p in terms of the Fermionic creationand annihilation operators [77] J + p = a † ( p ) a ( p ) , (2) J − p = a † ( p ) a ( p ) , (3) J z p = 12 (cid:16) a † ( p ) a ( p ) − a † ( p ) a ( p ) (cid:17) . (4)An analogous set of weak isospin operators may bedefined in the flavor basis, using the corresponding cre-ation/annihilation operators, which are related to theirmass basis counterparts via the following unitary trans-formation a e ( p ) = cos θ a ( p ) + sin θ a ( p ) (5) a x ( p ) = − sin θ a ( p ) + cos θ a ( p ) . (6)The isospin operators form an SU(2) algebra, obeyingthe usual commutation relations[ J + p , J − q ] = 2 δ pq J z p , [ J z p , J ± q ] = ± δ pq J ± p . (7)At this point, it is important to note that the strengthof the interaction between any two neutrinos is depen-dent on the intersection angle between their momenta,as can seen from the second term of the Hamiltonian inEq. (1). This geometric dependence makes the collectiveneutrino oscillation problem extremely complex, even inthe mean field limit. To avoid these complexities, and to facilitate a qualitative understanding of various collectiveflavor evolution phenomena, the so-called “single-angle”approximation has been frequently adopted in the litera-ture, wherein the angle-dependent interaction is insteadreplaced by an overall, angle-averaged coupling strength.The implication of this approximation is that it removesany trajectory-dependence from the flavor evolution ofneutrinos. It has been demonstrated that the single-anglecalculations are able to capture many of the qualita-tive behaviors observed in the more sophisticated multi-angle calculations. To this end, we define a direction-independent weak isospin operator (cid:126)J ω as follows: (cid:126)J ω = (cid:88) | p | = δm ω (cid:126)J p . (8)The self-interaction term may then be approximatedas H νν ≈ √ G F V (cid:104) − cos ϑ pq (cid:105) (cid:126)J · (cid:126)J ≡ µ ( r ) (cid:126)J · (cid:126)J, (9)where we have defined (cid:126)J = (cid:80) ω J ω = (cid:80) p J p . The radialdependence of the coupling strength µ arises due to thegeometric dilution of the neutrino fluxes and the narrow-ing of the intersection angles as one moves further fromthe source. For example in a “neutrino bulb” model,where neutrinos are assumed to be emitted isotropicallyfrom a single spherical emission surface, the dependenceis given by µ ( r ) ∝ (cid:32) − (cid:114) − R ν r (cid:33) , (10)where R ν is the radius of the neutrinosphere. III. DIAGONALIZATION AND THE BETHEANSATZ EQUATIONS
In the single-angle picture, the Hamiltonian fromEq. (1) may be re-written as H = (cid:88) p ω p (cid:126)B · (cid:126)J p + µ ( r ) (cid:126)J · (cid:126)J, (11)where p is now just an index denoting the oscillation fre-quencies present in the system. It has been shown thatthis particular Hamiltonian is amenable to diagonaliza-tion using the Richardson-Gaudin procedure. To beginwith, one may observe that, among the common eigen-states | j, m (cid:105) of the total weak isospin operators (cid:126)J and J z , some can be written as direct products of eigenstates | j p , ± j p (cid:105) of (cid:126)J p and J zp . For instance, the states | j, + j (cid:105) ≡ (cid:79) p | j p , + j p (cid:105)| j, − j (cid:105) ≡ (cid:79) p | j p , − j p (cid:105) , (12)where j = (cid:80) p j p , can be shown to be simultaneous eigen-states of (cid:126)J and each J zp , and are therefore eigenstatesof the Hamiltonian . In particular, the highest and low-est weight states | N , N (cid:105) and | N , − N (cid:105) , where N is thetotal number of neutrinos in the system, are eigenstates.These represent the states where all the neutrinos areeither isospin-up or isospin-down, i.e., | ν , . . . , ν (cid:105) and | ν , . . . , ν (cid:105) , respectively. The corresponding eigenvaluesfor these two states can be easily shown to be E ± N/ = ∓ (cid:88) p ω p N p µ N (cid:18) N (cid:19) . (13)where N p is the number of neutrinos at the oscillationfrequency ω p , and where we have suppressed the radialdependence in the notation of µ for convenience. Theremaining eigenstates of the Hamiltonian may then besystematically constructed by first defining the Gaudinalgebra with the operators (cid:126)S ( ζ α ) ≡ (cid:88) p (cid:126)J p ω p − ζ α (14)where we have introduced a sequence of Bethe Ansatzvariables ζ α , which are yet to be determined, followingthe formalism introduced in [78].Our Bethe Ansatz is the claim that the eigenstates ofour system are states of the form | ζ , . . . , ζ κ (cid:105) = N ( ζ , . . . , ζ κ ) (cid:18) κ (cid:89) α =1 S − α (cid:19) | j, + j (cid:105) (15)where | j, + j (cid:105) is an eigenstate, S − α = S − ( ζ α ), and N ( ζ , . . . , ζ κ ) is a normalization factor. An eigenstateof the Hamiltonian, such as the one defined in Eq. (15),is typically a linear superposition of several mass basisstates (i.e., Kronecker-products of mass states of individ-ual neutrinos), each consisting of N/ j − κ neutrinos inthe ν state and N/ − j + κ neutrinos in the ν state. Thefact that each eigenstate has a well-defined number of ν and ν is a consequence of the total J z = (cid:80) p J zp commut-ing with the Hamiltonian H from Eq. (11). Eigenstatesof the Hamiltonian are therefore also eigenstates of J z : J z | ζ , . . . , ζ κ (cid:105) = ( j − κ ) | ζ , . . . , ζ κ (cid:105) . (16) Note that, for j < N/
2, not every state of the form | j, ± j (cid:105) is aneigenstate. j = (cid:80) p j p is a necessary condition. A similar procedure may be carried out in which we in-stead apply S + α operators to | j, − j (cid:105) to obtain such eigen-states. In fact, these two procedures yield identical re-sults. This fact can be shown by defining a rotation op-erator T = e i π ( J + + J − ) , for which T J z T − = − J z and T J ± T − = J ∓ , and so T | j, − j (cid:105) = | j, + j (cid:105) . Moreover,one can see that, if | ψ (cid:105) is an eigenstate obtained fromEq. (15), then T | ψ (cid:105) is an eigenstate of a Hamiltonian T HT − = − (cid:80) p ω p (cid:126)B · (cid:126)J p + µ ( r ) (cid:126)J · (cid:126)J that may be ob-tained from a “raising formalism” using S + α instead.For a particular eigenstate consisting of κ neutrinosin the isospin-down ( ν ) configuration, we will have κ different ansatz variables ζ α to determine, and one canshow [27, 49] using the commutation relations betweenthe Gaudin operators that this requirement is equivalentto the condition − µ − M (cid:88) p =1 j p ω p − ζ α = κ (cid:88) β =1 β (cid:54) = α ζ α − ζ β , (17)for α = 1 , . . . , κ . Here, M is the total number of energy(or equivalently, ω ) values in the system, and where j p is related to the eigenvalue of the Casimir operator (cid:126)J p (each j p can take values from 0 or 1 / N p / κ , and j p s, Eq. (17) can admitmultiple solutions of the form { ζ , . . . , ζ κ } . Obtainingthe complete set of eigenstates and eigenvalues for an N -neutrino system therefore involves solving multiple setsof equations of the form Eq. (17), for κ = 1 , . . . , N andeach j p = 0 or 1 / , . . . , N p /
2. For the specific case where j p = 1 / p , the number of solutions of Eq. (17)for a given κ is equal to N C κ , and the total number ofsolutions is 2 N , across all values of κ .One can always choose j = N/ E = E + N/ + κ (cid:88) α =1 ζ α − µκ ( N − κ + 1) . (18) IV. THE LAMBDA METHOD
The Bethe Ansatz equations, Eq. (17), constitute a setof coupled algebraic equations in κ variables. In princi-ple, numerical solutions to this set of equations may besought—however, in their original form the equations areunwieldy. For starters, the variables ζ α admit complexvalues, and the equations contain singularities for certainvalues of parameters ω p and µ where the different ζ α ap-proach each other. Moreover, if each equation were tobe converted into a coupled polynomial form (by cross-multiplying all the denominators), then the order of eachpolynomial would be M + κ −
2. Therefore, it is worth-while to explore the possibilities of recasting the BetheAnsatz equations into a different, more tractable form.In order to accomplish this, one may introduce certainauxiliary functions which depend on the Bethe Ansatzvariables. For instance, following Refs. [84, 85] we maydefine the functionΛ( λ ) = κ (cid:88) α =1 λ − ζ α , (19)and transform Eq. (17) into a first-order ordinary differ-ential equationΛ( λ ) + Λ (cid:48) ( λ ) + 1 µ Λ( λ ) = M (cid:88) q =1 j q Λ( λ ) − Λ( ω q ) λ − ω q , (20)where the prime represents derivative with respect to λ .Eq. (20) is not straightforward to integrate because of thepresence of the parameters Λ( ω q ), whose values are notknown a priori , and are in fact dependent on the equationitself. These parameters can be determined by taking thelimit of Eq. (20) as λ → ω p for each p = 1 , . . . , M . Doingso yields the following system of equations:Λ p + (1 − j p )Λ (cid:48) p + 1 µ Λ p = M (cid:88) q =1 q (cid:54) = p j q Λ p − Λ q ω p − ω q , (21)where Λ p = Λ( ω p ) and Λ (cid:48) p = Λ (cid:48) ( ω p ), for p = 1 , . . . , M .In particular, if j p = 1 / p , then our equations forΛ p reduce simply to the formΛ p + 1 µ Λ p = M (cid:88) q =1 q (cid:54) = p Λ p − Λ q ω p − ω q , (22)yielding a system of coupled algebraic equations ofquadratic order in the parameters Λ p . Physically, thisrepresents the case where the system is composed of neu-trinos that all have pairwise distinct momenta, allowingus to choose a discrete set of ω p bins in which each binincludes exactly one neutrino. Details of the derivationof Eqs. (20) and (22) are given in Appendix A.Eq. (22) is manifestly much simpler than the origi-nal Bethe Ansatz equations Eq. (17). Moreover, unlikeEq. (17) where each value of κ = 1 , . . . , N requires solvinga separate set of Bethe Ansatz equations, Eq. (22) repre-sents just a single set of equations that can be solved toyield all the solutions corresponding to different valuesof κ . In Sec. IV A, we show that the following relationholds between the variables Λ p and ζ α : M (cid:88) p =1 j p ω p Λ p = κ M (cid:88) p =1 j p − µ κ (cid:88) α =1 ζ α − κ ( κ − , (23)which may be used to express our energy eigenvalues fromEq. (18) in terms of Λ p instead of the Bethe Ansatz vari-ables ζ α : E = E N/ − µ M (cid:88) p =1 j p ω p Λ p . (24) Thus, we seek to determine the parameters Λ , . . . , Λ M from Eq. (22). Following that, one may use Eq. (19),which may then be inverted to obtain the Bethe Ansatzvariables ζ α . The ζ α may then be used to reconstructthe states and their energies that solve our model usingEqs. (14) and (18). Alternatively, we show in this paperthat for the system in which j p = 1 / p , one mayinstead directly reconstruct the eigenstates and their en-ergies using the variables Λ p , without involving the BetheAnsatz variables. These two procedures are detailed inSecs. IV A and IV B, respectively. A. Obtaining the Bethe Ansatz roots ζ α from thevariables Λ p After solving Eq. (22) for Λ p , it is possible to reducethe problem of obtaining the Bethe Ansatz variables ζ α to that of solving a single polynomial equation of order κ .This process involves two steps, the first being deriving aset of constraint relations between the Λ p s and the powersums of ζ α s. Using the definition Eq. (19), one has M (cid:88) p =1 j p ω kp Λ p = κ (cid:88) α =1 M (cid:88) p =1 j p ω kp ω p − ζ α = (cid:88) α (cid:88) p j p ω kp − ζ kα ω p − ζ α + (cid:88) α ζ kα (cid:88) p j p ω p − ζ α = (cid:88) α (cid:88) p j p (cid:32) k (cid:88) l =1 ω k − lp ζ l − α (cid:33) − (cid:88) α ζ kα (cid:88) β (cid:54) = α ζ α − ζ β − (cid:88) α ζ kα µ . (25)Here, in the final step, we have used the Bethe Ansatzequation Eq. (17) to replace the second inner sum. Us-ing symmetry arguments and changing the order of thesummations, the above equation becomes M (cid:88) p =1 j p ω kp Λ p = k (cid:88) l =1 (cid:32) κ (cid:88) α =1 ζ l − α (cid:33) (cid:32) M (cid:88) p =1 j p ω k − lp (cid:33) − µ κ (cid:88) α =1 ζ kα − κ (cid:88) α,β =1 α (cid:54) = β ζ kα − ζ kβ ζ α − ζ β = k (cid:88) l =1 (cid:32)(cid:88) α ζ l − α (cid:33) (cid:32)(cid:88) p j p ω k − lp (cid:33) − µ (cid:88) α ζ kα − (cid:88) α,βα (cid:54) = β k (cid:88) l =1 ζ k − lα ζ l − β . (26)Here, it is useful to define the power sums of the BetheAnsatz variables, Q k ≡ (cid:80) κα =1 ζ kα . This allows us to writethe above expression as M (cid:88) p =1 j p ω kp Λ p = − µ Q k + k (cid:88) l =1 Q l − (cid:32) M (cid:88) p =1 j p ω k − lp (cid:33) − (cid:34) k (cid:88) l =1 Q k − l Q l − − k Q k − (cid:35) . (27)For the first few values of k , the above equation takesthe following forms: (cid:88) p j p Λ p = − Q µ = − κ µ , (28) (cid:88) p j p ω p Λ p = − Q µ + κ (cid:88) p j p − κ ( κ − , (29) (cid:88) p j p ω p Λ p = − Q µ + κ (cid:88) p j p ω p + Q (cid:88) p j p − ( κ − Q , (30)for k = 0 , , and 2, respectively. In particular, the equa-tion for k = 0 may be treated as a constraint relationfor the solutions Λ p obtained by solving Eq. (22), andcan also be used to classify those solutions according to κ . Also note that the equation for k = 1 is identical toEq. (23), which can be used to express the energy eigen-values in terms of the Λ p , as shown earlier.Using Eq. (27), one can successively calculate thepower sums Q k to any desired value of k , once all theΛ p are known. And the first κ power sums, Q , . . . , Q κ ,can be used to obtain all the roots ζ α . This involves firstcalculating the elementary symmetric polynomials of the ζ α s from the power sums. The elementary symmetricpolynomials are e = 1, e = (cid:80) α ζ α , e = (cid:80) α,β<α ζ α ζ β ,and so on. For k ≤ κ , these can be calculated recursivelyfrom the power sums using Newton’s identities: k e k ( ζ , . . . , ζ κ ) = k (cid:88) i =1 ( − i − e k − i Q i . (31)Therefore, one has e = e Q = Q , (32) e = 12 ( e Q − e Q ) , (33)and so on. Once the κ elementary symmetric poly-nomials e , . . . , e κ are evaluated, then the polynomial P ( λ ) ≡ (cid:81) κα =1 ( λ − ζ α ), whose roots are ζ α : α = 1 , . . . , κ may be constructed as P ( λ ) = κ (cid:88) k =0 ( − k e k λ κ − k . (34)Any one-dimensional polynomial root-finding algo-rithm can be employed to numerically obtain the roots ζ α of this polynomial. Once the ζ α are determined, thenthe eigenstates of the Hamiltonian may be explicitly con-structed using Eq. (15). As an aside, it is also interest-ing to note that P ( λ ) is related to the function Λ( λ ) asΛ( λ ) = P (cid:48) ( λ ) /P ( λ ). B. Constructing the eigenstates directly using Λ p Alternatively, we show that it is possible to directlycompute the eigenstates of the Hamiltonian in terms ofthe auxiliary variables, without needing to calculate theBethe Ansatz variables first. Unlike the procedure de-scribed in Sec. IV A, this method does not involve anyadditional numerical root-finding after obtaining the Λ p ,and therefore eliminates a potential source of numericalerror. We can use the fact [ J − p , J − q ] = 0 to find that theright hand side of Eq. (15) can be rewritten directly interms of Λ p (without involving the ζ α ) for any κ . Explicitderivations of these identities for κ = 2 , κ . Let us define the κ × κ matrix A with thematrix elements A ij = S − i δ ij , where δ ij is the Kroneckerdelta. Disregarding the normalization for the time being,one can rewrite Eq. (15) as | ζ , . . . , ζ κ (cid:105) = e κ ( S − , . . . , S − κ ) | j, + j (cid:105) = det( A ) | j, + j (cid:105) , (35)where e κ is the κ -th elementary symmetric polynomialand where we have used the fact that e κ ( x , . . . , x κ ) = (cid:81) κj =1 x j for any variables x j . Since A is a square matrixover a commutative ring, we can use the Cayley-HamiltonTheorem to infer that | ζ , . . . , ζ κ (cid:105) = 1 κ ! (cid:88) σ ∈ S( κ ) sgn( σ )tr σ ( A ) | j, + j (cid:105) (36)where S ( κ ) is the symmetry group of κ letters andtr σ ( A ) = tr( A f ) · · · tr( A f n ) for a permutation σ of cycletype ( f , . . . , f n ) [86]. The traces may be expressed interms of the Λ p s via the following steps:tr( A f ) = κ (cid:88) α =1 ( S − α ) f = κ (cid:88) i =1 M (cid:88) p =1 · · · M (cid:88) p f =1 J − p · · · J − p f ( ω p − ζ α ) . . . ( ω p f − ζ α )= κ (cid:88) α =1 M (cid:88) p =1 · · · M (cid:88) p f =1 J − p · · · J − p f f (cid:88) m =1 ω p m − ζ α f (cid:89) l =1 l (cid:54) = m ω p l − ω p m = M (cid:88) p =1 · · · M (cid:88) p f =1 J − p · · · J − p f f (cid:88) m =1 Λ p m f (cid:89) l =1 l (cid:54) = m ω p l − ω p m . (37)Here we have used the Heaviside cover-up rule betweenthe second and third equalities, and the definition fromEq. (19) between the third and fourth equalities (after ex-changing the order of summation). An alternative deriva-tion for the individual overlaps of energy eigenstates withthe mass basis states is given in Ref. [87, 88].Since the form of sgn( σ )tr σ ( A ) depends on the cycletype—which have multiplicities c σ —but not the particu-lar σ , we can reduce Eq. (36) to the following: | ζ , . . . , ζ κ (cid:105) = 1 κ ! (cid:88) cycle types sgn( σ ) c σ tr σ ( A ) | j, + j (cid:105) . (38)Lastly, we note that reversing the order of l and m in the denominator of Eq. (37) for each factor tr( A f )in tr σ ( A ) produces a factor of sgn( σ ), which cancelswith the same factor already written in each term ofEq. (36). In fact, using a decomposition into cycle-types, σ = ( f , . . . , f n ), we see that sgn( σ ) = (cid:81) ni =1 ( − f i +1 =( − κ + n , while a factor of (cid:81) ni =1 ( − f i − = ( − κ − n isproduced from reversing the differences mentioned in theprevious sentence. Alternatively, since tr( A f ) is simply the power sum (cid:80) κi =1 ( S − i ) f , one may use Newton’s identities to system-atically construct the elementary symmetric polynomial e κ ( S − , . . . , S − κ ) using the traces tr( A ) , . . . , tr( A κ ). Thiswas the basis of our numerical approach to calculatingthe eigenstates, the results of which are shown in Sec. V. V. SOLUTIONS OF THE BETHE ANSATZEQUATIONS
Having obtained a set of algebraic equations in theauxiliary variables Λ p , we are now in a position to dis-cuss the solutions to these equations. To begin with, onecan examine the solutions in the limit µ →
0, which we Interestingly, the resulting form of Eq. (36) with this cancella-tion suggests that the module spanned by S − , . . . , S − κ forms a κ -th-power symmetric (as opposed to exterior) algebra, whosecharacter can be described with the same formula. henceforth shall also refer to as the “asymptotic limit.”At this point, it is instructive to define the quantities (cid:101) Λ p ≡ µ Λ p to rewrite Eq. (22) as (cid:101) Λ p + (cid:101) Λ p = µ M (cid:88) q =1 q (cid:54) = p (cid:101) Λ p − (cid:101) Λ q η pq . (39)For convenience, we have defined η pq = ω p − ω q . Tak-ing the limit µ → (cid:101) Λ p = 0 or −
1, indepen-dently for each p . These asymptotic solutions can serveas a starting point for efficiently calculating numerical so-lutions at a generic µ >
0. This also makes it easy to seethat for the case where j p = 1 / p , 2 M solutionsexist in total. In fact, one can infer that the number ofsolutions is 2 M even when µ >
0, by noting that Eq. (39)is a set of M mutually co-prime coupled quadratic equa-tions in M variables, and therefore has a finite solutionset. One can then invoke Homotopy continuation to ar-gue that each solution for a generic µ > µ → A. Algebraic Solutions for M = 2 For the purposes of gaining some mathematical andphysical intuition, we first show the analytic solutionsfor a simple system consisting only of two interactingneutrinos, at frequencies ω and ω . This correspondsto solving Eq. (22) for M = 2 and can serve as a testbed for the numerical technique that we later implementfor dealing with the M > p for the model with M = N = 2 with j = j = 1 /
2. These Λ p can easily befound analytically from Eq. (22). The four solutions arelisted in Table I, and categorized by their corresponding κ values. For a particular solution { Λ p : p = 1 , . . . , M } , TABLE I. Algebraic solutions for Λ p from Eq. (22) with M = 2.Λ Λ κ = 0 0 0 κ = 1 − µ + 1 η − sgn( η ) (cid:115) η + 14 µ − µ − η + sgn( η ) (cid:115) η + 14 µ − µ + 1 η + sgn( η ) (cid:115) η + 14 µ − µ − η − sgn( η ) (cid:115) η + 14 µ κ = 2 − µ − µ TABLE II. Energy eigenvalues and eigenstates determined from Eq. (22) M = 2. E | E (cid:105) κ = 0 2 µ −
12 ( ω + ω ) | ν , ν (cid:105) κ = 1 − ( ω + ω ) − µ (cid:20) (cid:115) ω − ω ) µ (cid:21) N , + (cid:0) Λ , + | ν , ν (cid:105) + Λ , + | ν , ν (cid:105) (cid:1) − ( ω + ω ) − µ (cid:20) − (cid:115) ω − ω ) µ (cid:21) N , − (cid:0) Λ , − | ν , ν (cid:105) + Λ , − | ν , ν (cid:105) (cid:1) κ = 2 2 µ + 12 ( ω + ω ) | ν , ν (cid:105) κ can be determined using the identity M (cid:88) p =1 Λ p = − κµ , (40)which was derived in Sec. IV A (Eq. (28)). Note that,across a complete set of solutions, κ takes all values from0 , , . . . , N , and Eq. (40) holds for arbitrary M = N .Furthermore, using the algebraic Λ p solutions fromTable I, we may compute the energy eigenvalues usingEqs. (13) and (24) and eigenstates using Eqs. (36) and(37) (or more specifically, Eq. (C3) for M = 2). Theseeigenvalues and eigenstates are aggregated in Table II.Note that the Λ , ± and Λ , ± referenced in the latter ta-ble for κ = 1 are respectively from the two solutions inTable I for κ = 1. Additionally, normalization coefficientscalculated for states with κ = 1 , |N , ± | = 4 (cid:20)(cid:18) µ + 1 η (cid:19) ∓ | η | (cid:115) µ + 1 η (cid:21) (41) |N | = 1 µ (42)respectively. Note that even though these normalizationcoefficients are singular as µ →
0, the eigenstates them-selves are not.
B. Numerical Solutions for
M > Extending from our argument of the previous sectionthat demonstrated there must be finitely many solutionsto Eq. (22) (or, equivalently, to Eq. (39)), we can further observe that the problem with
M > M ≥ µ → µ → M solutions in thelimit µ → µ >
0. Infact, we build upon this existing homotopy method forour purposes: if 0 ≤ t ≤ t = 1, but also all solutions determined for intermediatevalues 0 < t < µ . We apply Newton’s method tosolve our system of equations at the n th homotopy step: (cid:126) F ( (cid:101) Λ ( n )1 , . . . , (cid:101) Λ ( n ) M , t ( n ) ) = (cid:126) , (43) − − − − ω Λ p µ/ω FIG. 1. One out of the 210 solutions with κ = 6 to Eq. (22) for a system with M = N = 10 and ω p = pω , calculatednumerically using the modified Newton-Raphson method with homotopy continuation. Observe that, as µ →
0, six out of theten Λ p values approach −∞ (as − /µ ), as expected for a κ = 6 solution. for each t ( n ) >
0, where (cid:126) F is defined in Eq. (B3). Atthe n -th Homotopy step, the initial guess for Newton’smethod is taken to be the numerical solution (cid:101) Λ ( n − p ob-tained by solving the system of equations at the ( n − th step, i.e., (cid:126) F ( (cid:101) Λ ( n − , . . . , (cid:101) Λ ( n − M , t ( n − ) = (cid:126) . (44)Given the form of (cid:126) F in Eq. (B3), we can easily computethe analytic form of the Jacobian ∂ (cid:126)z (cid:126) F ( (cid:126)z , t ), to used ineach step of Newton’s method. Finally, all 2 M solutionsto our system can be obtained by repeating the proce-dure, starting from each (cid:101) Λ p independently taking values0 or − µ = 0.To test this method, we computed the numerical solu-tions to Eq. (22) for M = 2, and compared them with theanalytic results described in Sec. V A. The root-mean-squared relative error between the analytic and numeri-cal values in comparing each component of each solutionwas (cid:46) − , suggesting that the method was able toattain a high level of numerical accuracy.A potential issue with the utilization of homotopy con-tinuation with large- M -dimensional algebraic systems isthat the close proximity of solutions in R M may causenumerical 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 dur-ing each Newton-Raphson step helps in obtaining theexpected continuity of solutions as µ is varied. We applied this method to systems of interacting neu-trinos with equally spaced oscillation frequencies givenby when ω p = pω , just that only a single neutrino re-sides at each oscillation frequency (therefore, j p = 1 / p and so M = N ). Such a system of N neutrinosadmits 2 N solutions. Solutions of the system for up to N ∼
15 can be computed within a reasonable time frameusing a personal computer. In Fig. 1, we show the evolu-tion with µ of one of the solutions to such a system with N = 10 neutrinos. Shown in the figure are the quantitiesΛ p = (cid:101) Λ p /µ , for p = 1 , . . . , κ = 0 , . . . ,
5, for the same system as in Fig. 1. We alsonote that results for κ > S + α operators from our Gaudin algebra to the state | j, − j (cid:105) . In analogy to (cid:101) Λ p , one can define the variables (cid:101) Λ ( ↑ ) p in the raising formalism, which obey the coupledquadratic equations (cid:101) Λ ( ↑ ) p − (cid:101) Λ ( ↑ ) p = µ M (cid:88) q =1 q (cid:54) = p (cid:101) Λ ( ↑ ) p − (cid:101) Λ ( ↑ ) q η pq . (45)Similar to Eq. (40), the solutions (cid:101) Λ ( ↑ ) p obey the constraintrelations M (cid:88) p =1 j p Λ ( ↑ ) p = + κ ( ↑ ) µ (46)where Λ ( ↑ ) p = (cid:101) Λ ( ↑ ) p /µ . In particular, a solution (cid:101) Λ ( ↑ ) p ofEq. (45) can be shown to correspond to a solution (cid:101) Λ ( ↓ ) p of Eq. (39), via the identity (cid:101) Λ ( ↑ ) p = − (cid:101) Λ ( ↓ ) p . (47)Using Eq. (40), it can be shown that the corresponding κ ( ↑ ) of a solution in the raising formalism can be relatedto a given κ ( ↓ ) via κ ( ↑ ) = M − κ ( ↓ ) , further reinforcingthis duality. This correspondence of solutions for a given κ and M − κ is also discussed in Refs. [49, 78, 87]. Inparticular, Ref. [87] notes that, in the limit µ → ∞ ,Λ ( ↑ ) p → Λ ( ↓ ) p .For each value of κ , one can observe that the en-ergy eigenvalues form b + 1 distinct branches, where b ≡ min { κ, M − κ } , in the limit of large µ , as previ-ously noticed in [49]. One may also observe that theself-interaction term in the Hamiltonian from Eq. (11) be-comes dominant as µ → ∞ , and therefore the eigenstatesof H approximately align with the eigenstates | j, m (cid:105) of (cid:126)J · (cid:126)J , with eigenvalues E ≈ µj ( j +1). For an energy eigen-state | ψ (cid:105) with a given κ , we can see that J z | ψ (cid:105) = m | ψ (cid:105) where m = M/ − κ ; therefore as µ → ∞ , the quantumnumber j of all the eigenstates for that κ may take valuesin the range M/ , . . . , M/ − b . So, as µ → ∞ , we expectthat the energy eigenvalues of states with a given κ willsplit into b + 1 branches with the approximate energiesgiven above.In Fig. 3a, we zoom in on the small- µ regions of Fig. 2f(all eigenvalues of the κ = 5 eigenstates) to better illus-trate the various energy-level crossings that are presentin this region. To examine the nature of these level cross-ings in better detail, we show in Fig. 3b the energy eigen-values of all the states (spanning all permitted κ values)in a system with N = M = 3. Eigenstates correspondingto each κ are assigned a particular color, as described inthe figure caption. Ref. [49] observes that the highest-energy states from among the solutions of each κ (high-lighted here as dashed lines) do not have energy levelcrossings with other states of the same κ , which is con-sistent with our observation. Here, we would additionallylike to point out that the highest energy eigenstate of aparticular κ can cross with other eigenstates of a different κ . However, since J z is a symmetry of the Hamiltonian,and each κ corresponds to a unique J z eigenvalue (seeEq. (16)), these crossings do not result in mixing betweenthese eigenstates.In addition to our energy eigenvalues, we may calcu-late energy eigenstates for generic µ using Eqs. (35)–(37).These eigenstates may be encoded as a sequence of their2 M coefficients in the mass basis. Particular examplesof eigenstates for systems with for M = 3 , µ ) of an arbitrary initialstate, in the adiabatic limit. This can be used for thepurposes of comparing the results of a many-body cal-culation to the corresponding results in the mean-fieldlimit, as well as for studying exclusive many-body effectslike the emergence of quantum entanglement between thevarious modes as the system evolves. VI. CONCLUSIONS
Adiabatic evolution of a many-neutrino system in thesingle-angle approximation is exactly solvable in thesense that such an evolution can be completely charac-terized by the solutions of the appropriate Bethe ansatzequations. However, solving those nonlinear Betheansatz equations is a highly nontrivial problem. In thispaper we presented a technique to evaluate the exactadiabatic eigenvalues and eigenstates of the collectiveneutrino oscillation Hamiltonian by casting Bethe ansatzequations first into a differential form, then into a set ofalgebraic equations which are numerically more tractablethan the original Bethe ansatz equations. With our out-lined procedure, to determine solutions for up to N ∼ µ ≤ N would requiremore computational power. In the future we will explorethe behavior of the solutions as N gets larger.An immediate benefit of obtaining such exact solutionsis the ability to explore the limits of applicability of thecommonly used mean-field solutions. Notably, our ex-act problem (with j p = 1 /
2) has 2 N solutions, whilethe mean-field problem has only 2 N solutions in total.Clearly, in the mean-field case we are either losing or com-bining many states. Earlier studies of the question of en-tangled neutrinos also explored modifications of the col-lective oscillations due to the many-body effects presentwhen one goes beyond the one-body description inher-ent to the mean-field approximation [74–76]. These pa-pers primarily investigate part of the Hamiltonian whichsurvives in the large µ limit and explore if the oscilla-tions speed up due to many-particle entanglements. Thegrowth rates they obtain differ by a factor of √ N rel-ative to each other, depending on the setup. In futurepublications we plan to explore many-body entanglementeffects using our approach, for generic values of µ whereboth the one-body and two-body terms in the Hamilto-nian play a role. Clearly exploring these issues is a stepwhich needs to be taken before moving on to astrophysicsapplications.Many of the elements heavier than iron were formed byrapidly capturing neutrons on seed nuclei (r-process nu-cleosynthesis). In the astrophysical site of the r-processnucleosynthesis many reactions take place during a rathershort duration, a feature which is typically associatedwith explosive phenomena. Currently the leading candi-dates for the sites of r-process nucleosynthesis are core-0 − . . . .
53 0 2 4 6 8 10 − E / ω µ/ω (a) κ = 0 − . . . .
53 0 2 4 6 8 10 − E / ω µ/ω (b) κ = 1 − . . . .
53 0 2 4 6 8 10 − E / ω µ/ω (c) κ = 2 − . . . .
53 0 2 4 6 8 10 − E / ω µ/ω (d) κ = 3 − . . . .
53 0 2 4 6 8 10 − E / ω µ/ω (e) κ = 4 − . . . .
53 0 2 4 6 8 10 − E / ω µ/ω (f) κ = 5 FIG. 2. Energy eigenvalues calculated from Eq. (24) after obtaining the values of each Λ p . System parameters are identical tothose in Figure 1. The solutions for κ > M − κ counterparts, up toa constant vertical offset. − . − . − . . . . . . . . . . . . . . . . − E / ω µ/ω (a) − − . . E / ω µ/ω (b) FIG. 3. In the region of small µ , we observe numerous level crossings between different energy states of the Hamiltonian inEq. (11). Left: energy eigenvalues of all the κ = 5 solutions for the case N = M = 10, with the same conditions as in Figure2f. Right: energy eigenvalues for all the solutions of a system with N = M = 3, with eigenvalues corresponding to different κ coded as follows: κ = 0 (dotted line), 1 (solid lines), 2 (dashed lines), and 3 (dot-dashed line). . . . .
81 0 2 4 6 8 10 | h ν i ,..., ν i | ψ i | µ/ω (a) The mass-basis decomposition for the fourth excited state | ψ (cid:105) in the N = M = 3 system, as a function of µ . Here, κ = 1. . . . .
81 0 2 4 6 8 10 | h ν i ,..., ν i | ψ i | µ/ω (b) The mass-basis decomposition of the third excited state | ψ (cid:105) in the N = M = 4 system, as a function of µ . Here κ = 2. FIG. 4. Overlaps, | (cid:104) ν i , . . . , ν i N | ψ n (cid:105) | , of an excited state | ψ n (cid:105) with the mass basis states (with i , . . . , i N = 1 , N . Here, the n th energy eigenstate | ψ n (cid:105) ( n = 0 , . . . , N −
1) is connected to the mass eigenstate | ν j , . . . , ν j N (cid:105) in the limit µ →
0, where j N +1 − k = 1 + ( k th digit of n in binary representation). Observe that a state with a given κ has N C κ nontrivial components for µ >
0. For Figure 4b, two of the six nontrivial overlaps are numerically indistinguishable.
ACKNOWLEDGMENTS
We would like to thank E. Armstrong, S. Birol, S. Cop-persmith, G. Fuller, E. Grohs, A. Hashimoto, C. John-son, E. Rrapaj, J. Schmidt, M. Sen, and S. Shalgar forvaluable conversations. This work was supported in partby the U.S. National Science Foundation Grants PHY-1630782 and PHY-1806368.
Appendix A: Deriving the Lambda equations fromthe Bethe Ansatz
In this appendix, we shall present the derivation of thecoupled quadratic equations in Lambda, Eq. (22), fromthe Bethe Ansatz equations Eq. (17). These results aredispersed through the condensed-matter physics litera-ture. Here for the convenience of the reader we gatherthem in one place. Using the function Λ( λ ) defined inEq. (19), one can writeΛ ( λ ) + Λ (cid:48) ( λ ) = κ (cid:88) α,β =1 α (cid:54) = β λ − ζ α )( λ − ζ β ) . (A1)Using partial fraction decomposition and symmetry ar-guments, this may be rewritten asΛ ( λ ) + Λ (cid:48) ( λ ) = 2 κ (cid:88) α =1 λ − ζ α κ (cid:88) β (cid:54) = α ζ α − ζ β ≡ κ (cid:88) α =1 W ( ζ α ) λ − ζ α . (A2)Replacing W ( ζ α ) with the left hand side of the BetheAnsatz equations Eq. (17), one obtainsΛ ( λ ) + Λ (cid:48) ( λ ) = 2 κ (cid:88) α =1 λ − ζ α (cid:34) − µ − (cid:88) p j p ω p − ζ α (cid:35) . (A3)Using the definition of Λ( λ ) from Eq. (19) after chang-ing the order of summation and using partial fractiondecomposition, this equation reduces toΛ ( λ ) + Λ (cid:48) ( λ ) + 1 µ Λ( λ ) = 2 (cid:88) p j p Λ( λ ) − Λ( ω p ) λ − ω p . (A4) It can be shown that the ordinary differential equa-tion, Eq. (A4), is exactly equivalent to the Bethe Ansatzequations, i.e., every solution of Eq. (A4) correspondsto a unique solution of the Bethe Ansatz equations, andvice versa. This equivalence can be proven using the factthat every step of the above derivation is reversible—onecould just as easily start with Eq. (A4) and derive theBethe Ansatz equations, Eq. (17).The ordinary differential equation in Eq. (A4) may beconverted to a set of algebraic equations in Λ( ω p ), foreach ω p in the system. In order to do so, one can Taylor-expand Λ( λ ) around λ = ω q . Following this Taylor ex-pansion, the right-hand side of Eq. (A4) may be writtenas 2 (cid:88) p j p Λ( λ ) − Λ( ω p ) λ − ω p = 2 (cid:88) p (cid:54) = q j p Λ( λ ) − Λ( ω p ) λ − ω p +2 j q (cid:20) Λ (cid:48) ( ω q ) + 12 Λ (cid:48)(cid:48) ( ω q )( λ − ω q ) + . . . (cid:21) (A5)where Λ (cid:48) ( ω q ) , Λ (cid:48)(cid:48) ( ω q ) , . . . are the successive derivatives ofΛ( λ ) with respect to λ , as evaluated at λ = ω q . Usingthis expansion and then taking the limit λ → ω q , oneobtains the following set of equations (one for each ω q )Λ ( ω q ) + (1 − j q )Λ (cid:48) ( ω q ) + 1 µ Λ( ω q )= 2 (cid:88) p (cid:54) = q j p Λ( ω q ) − Λ( ω p ) ω q − ω p . (A6)If j q = 1 / q (corresponding to there beingjust one neutrino per bin), then the Λ (cid:48) term vanishesand the equations become purely algebraic (Eq. (22)),and can in principle be solved to obtain the solutions { Λ( ω q ) : q = 1 , . . . , M } . If j q > /
2, then one can takesuccessive derivatives of Eq. (A4) with respect to λ , asnoted in Refs. [84, 85, 90, 91]. For example, taking thefirst derivative of Eq. (A4) followed by the λ → ω q limitgives us2Λ( ω q )Λ (cid:48) ( ω q ) + (1 − j q )Λ (cid:48)(cid:48) ( ω q ) + 1 µ Λ (cid:48) ( ω q )= 2 (cid:88) p (cid:54) = q j p (cid:20) Λ (cid:48) ( ω q ) λ − ω p − Λ( ω q ) − Λ( ω p )( ω q − ω p ) (cid:21) . (A7)Now, if j q = 1, then the Λ (cid:48)(cid:48) term vanishes, and one mayuse Eqs. (A6) and (A7) to eliminate Λ (cid:48) ( ω q ) from the sys-tem of equations. Likewise, if j q = 3 /
2, one may obtainan additional equation by taking the second derivativeof Eq. (A4), and use it along with Eqs. (A6) and (A7)to eliminate Λ (cid:48) ( ω q ) and Λ (cid:48)(cid:48) ( ω q ). And so on, for highervalues of j q .3 Appendix B: Proof of the applicability of thehomotopy continuation method
To carry out this homotopy method, let us first rewriteEq. (39) as the system (cid:126)F ( (cid:126)z ) = (cid:126) F p ( (cid:126)z ) = z p + z p − µ (cid:88) q (cid:54) = p z p − z q η pq (B1)for p = 1 , . . . , M . In particular, since we know each valueof (cid:101) Λ p must be real , we restrict (cid:126)z ∈ R . Additionally, wecan define a set of functions (cid:126)G with simpler solutions by G p ( (cid:126)z ) = z p + z p (B2)and the family of functions (cid:126) F ( (cid:126)z, t ) = (1 − t ) (cid:126)G ( (cid:126)z ) + t (cid:126)F ( (cid:126)z )for 0 ≤ t ≤
1. Then, from the explicit form of F p ( (cid:126)z, t ) = z p + z p − tµ (cid:88) q (cid:54) = p z p − z q η pq (B3)we see that (cid:126) F is an analytic (and therefore a C ) mappingof R M × [0 , → R M . By Sard’s Theorem, we find that (cid:126) F is a regular function. Recall that regular, in the con-text of algebraic geometry, means that: (i) (cid:126) F ( (cid:126)z , t ) = (cid:126) ∂ (cid:126)z (cid:126) F ( (cid:126)z , t ) has rank M for fixed t = 0 , (cid:126)z ∈ R M and (ii) (cid:126) F ( (cid:126)z , t ) = (cid:126) ∂ (cid:126)z,t (cid:126) F ( (cid:126)z , t ) hasrank M for fixed 0 < t < (cid:126)z ∈ R M . (Toclarify, ∂ (cid:126)z,t (cid:126) F is the total Jacobian of (cid:126)H while ∂ (cid:126)z (cid:126) F is theJacobian of (cid:126) F stripped of its final column, ∂ t (cid:126) F .) Fur-ther, a theorem by Garcia and Zangwill [89] states thatif, in addition, for all sequences with || (cid:126)z ( n ) || → ∞ thereis a p such that F p , G p satisfylim m →∞ F p ( (cid:126)z ( n m ) ) G p ( (cid:126)z ( n m ) ) ≮ α (B4)for some α > (cid:126)z ( n m ) atwhich G p ( (cid:126)z ( n m ) ) (cid:54) = 0, then all roots of (cid:126)F can be foundone-to-one from the known roots of (cid:126)G (i.e., the asymp-totic solutions as µ →
0) by homotopy curves as we vary t = 0 to t = 1. The symbol ≮ implies that the real part ofthe limit (which generically may be complex if our chosen (cid:126)F were not strictly real) must be greater than or equalto a positive α that we are free to choose. We can seethat this additional hypothesis is true for our system byfirst calculating for an arbitrary p that F p ( (cid:126)z ( n ) ) G p ( (cid:126)z ( n ) ) = 1 − µ (cid:80) q (cid:54) = p ( z ( n ) p − z ( n ) q ) /η pq z ( n )2 p + z ( n ) p (B5) Since the eigenvalues of a Hermitian operator are real (cf. Eqs.(23) and (24)), any complex-valued solutions ζ α of Eq. (17) mustcome in complex conjugate pairs [49, 78], implying by Eq. (19)that Λ( ω p ) and therefore (cid:101) Λ( ω p ) must be real. for any sequence of points (cid:126)z ( n ) after eliminating any ofthe points in the sequence solving G p ( (cid:126)z ) = 0 for some p . (We know we may choose such a p since we knowthere are finitely many points solving G p ( (cid:126)z ) = 0 for all p .) Moreover, we may observe that this sequence musthave at least one r for which | z ( n ) r | → ∞ at the fastestrate amongst the (cid:126)z ( n ) . Selecting p = r we find that F p ( (cid:126)z ( n ) ) /G p ( (cid:126)z ( n ) ) → (cid:126)F is homotopic to (cid:126)G as wevary tµ from 0 to µ , in the sense described by Ref. [89],and so all our solutions to (cid:126)F ( (cid:126)z ) = (cid:126) (cid:126)G ( (cid:126)z ) = (cid:126) (cid:101) Λ( ω p ) solutions.Additionally, we confirm from this process that theremust be exactly 2 M solutions to Eq. (39) for µ (cid:54) = 0 aswell. Moreover, for a given κ = − (cid:80) Mp =1 (cid:101) Λ p , there are M C κ solutions. We can observe that there are equally asmany ways for us to permute the values 0 , − (cid:101) Λ p for a given κ . Thus, we expect that byconsidering all possible permutations of 0 , − (cid:101) Λ p at µ = 0 we may find all (cid:80) Mκ =0 M C κ = 2 M solutions.Lastly, consider the situation in which we fix the val-ues of each ω p while we allow µ ≥ µ (i.e., tµ with t = 1)in Eq. (B3) is simply an arbitrary real µ >
0. So,we expect that we may obtain not only the numeri-cal solutions to (cid:126)F ( (cid:126)z ) = (cid:126) F ( (cid:126)z,
1) = (cid:126) n th ( n = 0 , . . . ) intermediate numerical solutions (cid:126)z ( n ) to (cid:126) F ( (cid:126)z ( n ) , t ( n ) ) = (cid:126) ≤ t ( n ) ≤
1, provided that oursteps t ( n +1) − t ( n ) > (cid:126)F | µ = µ max , we mayobtain the solutions to Eq. (39) for all 0 ≤ µ ≤ µ max withan arbitrary maximum µ max , as (cid:126)z ( n ) may be interpretedas the solutions to (cid:126)F | µ = t ( n ) µ max . Appendix C: Derivations of Eigenstates for κ = 2 , We can demonstrate the claim from Section IV B thatour general eigenstates can be written in terms of Λ p inlieu of ζ α for κ = 2 , j p = 1 / p . Here wepresent two ways to derive the eigenstate expressions for κ = 2 ,
3. First note that, for all κ , we have the identity:1 κ ! (cid:18) κ (cid:88) α =1 S − α (cid:19) κ | j, + j (cid:105) = 1 κ ! (cid:18) M (cid:88) p =1 Λ p J − p (cid:19) κ | j, + j (cid:105) = e κ (cid:0) Λ J − , . . . , Λ M J − M (cid:1) | j, + j (cid:105) (C1)4where e κ is the degree- κ elementary symmetric polyno-mial in M variables. Now, for κ = 2 we have (cid:88) α =1 ( S − α ) | j, + j (cid:105) = (cid:88) α =1 (cid:18) M (cid:88) p =1 J − p ω p − ζ α (cid:19)(cid:18) M (cid:88) p =1 J − p ω p − ζ α (cid:19) | j, + j (cid:105) = (cid:88) α =1 M (cid:88) p =1 M (cid:88) q =1 q (cid:54) = p ω p − ζ α )( ω − ζ α ) J − p J − q | j, + j (cid:105) = M (cid:88) p =1 M (cid:88) q =1 q (cid:54) = p (cid:88) α =1 (cid:18) ω p − ζ α − ω q − ζ α (cid:19) − ω p − ω q J − p J − q | j, + j (cid:105) = − M (cid:88) p,q =1 q (cid:54) = p Λ p − Λ q ω p − ω q J − p J − q | j, + j (cid:105) (C2) and so S − S − | j, + j (cid:105) = 12 (cid:20)(cid:18) (cid:88) α =1 S − α (cid:19) − (cid:88) α =1 ( S − α ) (cid:21) | j, + j (cid:105) = 12 M (cid:88) p,q =1 p (cid:54) = q (cid:18) Λ p Λ q + Λ p − Λ q ω p − ω q (cid:19) J − p J − q | j, + j (cid:105) . (C3)Thus, we have rewritten the eigenstates as well as theeigenvalues of the Hamiltonian in Eq. (11) completely in terms of our defined Λ( ω p ) parameters, bypassing theneed to directly compute the Bethe Ansatz variables ζ α .This process can be carried out for κ = 3 similarly: (cid:88) α =1 ( S − α ) | j, + j (cid:105) = (cid:88) α =1 M (cid:88) p,q,r =1 p,q,r distinct J − p J − q J − r ( ω p − ζ α )( ω q − ζ α )( ω r − ζ α ) | j, + j (cid:105) = M (cid:88) p,q,r =1 p,q,r distinct 3 (cid:88) α =1 (cid:20) − ω p − ω q (cid:18) ω p − ζ α − ω q − ζ α (cid:19) ω r − ζ α (cid:21) J − p J − q J − r | j, + j (cid:105) = M (cid:88) p,q,r =1 p,q,r distinct − ω p − ω q (cid:88) α =1 (cid:20) − ω p − ω r (cid:18) ω p − ζ α − ω r − ζ α (cid:19) − − ω q − ω r (cid:18) ω q − ζ α − ω r − ζ α (cid:19)(cid:21) × J − p J − q J − r | j, + j (cid:105) = M (cid:88) p,q,r =1 p,q,r distinct ω p − ω q (cid:18) Λ p − Λ r ω p − ω r − Λ q − Λ r ω q − ω r (cid:19) J − p J − q J − r | j, + j (cid:105) (C4)5and (cid:88) α =1 ( S − α ) (cid:88) β (cid:54) = α S − β | j, + j (cid:105) = (cid:88) α =1 M (cid:88) p,q,r =1 p,q,r distinct J − p J − q ( ω p − ζ β )( ω q − ζ γ ) (cid:18) J − r ω r − ζ k − J − r ω r − ζ l (cid:19) | j, + j (cid:105) = M (cid:88) p,q,r =1 p,q,r distinct 3 (cid:88) i =1 − ω p − ω q (cid:18) ω p − ζ α − ω q − ζ α (cid:19)(cid:20) Λ r − ω r − ζ α (cid:21) J − p J − q J − r | j, + j (cid:105) = − M (cid:88) p,q,r =1 p,q,r distinct (cid:20) Λ r Λ p − Λ q ω p − ω q − ω p − ω q (cid:88) α =1 (cid:18) ω p − ζ α − ω q − ζ α (cid:19) ω r − ζ α (cid:21) J − p J − q J − r | j, + j (cid:105) = − M (cid:88) p,q,r =1 p,q,r distinct (cid:20) Λ r Λ p − Λ q ω p − ω q + 1 ω p − ω q (cid:18) Λ p − Λ r ω p − ω r − Λ q − Λ r ω q − ω r (cid:19)(cid:21) J − p J − q J − r | j, + j (cid:105) (C5)where in the first line we are using β (cid:54) = γ and β, γ (cid:54) = α . Thus, S − S − S − | j, + j (cid:105) = 13! (cid:20)(cid:18) (cid:88) α =1 S − α (cid:19) − (cid:88) α =1 ( S − α ) − (cid:88) α =1 ( S − α ) (cid:88) β (cid:54) = α S − β (cid:21) | j, + j (cid:105) = 13! M (cid:88) p,q,r =1 p,q,r distinct (cid:20) Λ p Λ q Λ r + 3Λ r Λ p − Λ q ω p − ω q + 2 ω p − ω q (cid:18) Λ p − Λ r ω p − ω r − Λ q − Λ r ω q − ω r (cid:19)(cid:21) J − p J − q J − r | j, + j (cid:105) . (C6)There is a more direct approach to arrive at the same results, by decomposing terms of S − S − and S − S − S − imme-diately without consideration of other products – as we will do in the next section. Also, note that the form of theseresults bears similarity to the characters of exterior powers of vector spaces (i.e., Λ κ V for κ = 2 , (cid:81) κα =1 S − α in terms of Λ p more swiftly. Appendix D: Alternative Derivation of Eigenstates for κ = 2 , Here, we present a more direct method of evaluating the product (cid:81) κi =1 S − α | j, + j (cid:105) for the cases of κ = 2 ,
3. In κ = 2,we can quickly compute: S − S − | j, + j (cid:105) = 12 ( S − S − + S − S − ) | j, + j (cid:105) = 12 M (cid:88) p,q =1 (cid:20) ω p − ζ )( ω q − ζ ) + 1( ω p − ζ )( ω q − ζ ) (cid:21) J − p J − q | j, + j (cid:105) = 12 M (cid:88) p,q =1 p (cid:54) = q (cid:20) ω p − ζ ) (cid:18) Λ q − ω q − ζ (cid:19) + 1( ω p − ζ ) (cid:18) Λ q − ω q − ζ (cid:19)(cid:21) J − p J − q | j, + j (cid:105) = 12 M (cid:88) p,q =1 p (cid:54) = q (cid:18) Λ p Λ q + Λ p − Λ q ω p − ω q (cid:19) J − p J − q | j, + j (cid:105) . (D1)6Analogously, we can extend this argument to κ = 3: S − S − S − | j, + j (cid:105) = M (cid:88) p,q,r =1 p,q,r distinct J − p J − q J − r ( ω p − ζ )( ω q − ζ )( ω r − ζ ) | j, + j (cid:105) = M (cid:88) p,q,r =1 p,q,r distinct ω p − ζ )( ω q − ζ ) (cid:18) Λ r − ω r − ζ − ω r − ζ (cid:19) J − p J − q J − r | j, + j (cid:105) = 13! (cid:88) σ ∈ Sym(3) M (cid:88) p,q,r =1 p,q,r distinct ω p − ζ σ (1) )( ω q − ζ σ (2) ) (cid:18) Λ r − ω r − ζ σ (1) − ω r − ζ σ (2) (cid:19) J − p J − q J − r | j, + j (cid:105) (D2)while M (cid:88) p,q,r =1 p,q,r distinct Λ r ( ω p − ζ )( ω q − ζ ) J − p J − q J − r | j, + j (cid:105) → M (cid:88) p,q,r =1 p,q,r distinct Λ r (cid:20) ω p − ζ )( ω q − ζ ) + 1( ω p − ζ )( ω q − ζ )+ 1( ω p − ζ )( ω q − ζ ) (cid:21) J − p J − q J − r | j, + j (cid:105) = ⇒ M (cid:88) p,q,r =1 p,q,r distinct Λ r ( ω p − ζ )( ω q − ζ ) J − p J − q J − r | j, + j (cid:105) → M (cid:88) p,q,r =1 p,q,r distinct Λ r ω p − ζ (cid:18) Λ q − ω q − ζ (cid:19) J − p J − q J − r | j, + j (cid:105) = ⇒ (cid:88) σ ∈ Sym(3) M (cid:88) p,q,r =1 p,q,r distinct Λ r ( ω p − ζ σ (1) )( ω q − ζ σ (2) ) J − p J − q J − r | j, + j (cid:105) = M (cid:88) p,q,r =1 p,q,r distinct Λ r (cid:18) Λ p Λ q + Λ p − Λ q ω p − ω q (cid:19) J − p J − q J − r | j, + j (cid:105) (D3)and M (cid:88) p,q,r =1 p,q,r distinct ω p − ζ )( ω q − ζ )( ω r − ζ ) J − p J − q J − r | j, + j (cid:105) → (cid:88) σ ∈ Sym(3) M (cid:88) p,q,r =1 p,q,r distinct ω p − ζ σ (1) )( ω r − ζ σ (1) ) × (cid:18) ω q − ζ σ (2) + 1 ω q − ζ σ (3) (cid:19) J − p J − q J − r | j, + j (cid:105) = 112 (cid:88) σ ∈ Sym(3) M (cid:88) p,q,r =1 p,q,r distinct ω p − ζ σ (1) )( ω r − ζ σ (1) ) × (cid:18) Λ q − ω q − ζ σ (1) (cid:19) J − p J − q J − r | j, + j (cid:105) = 16 M (cid:88) p,q,r =1 p,q,r distinct (cid:18) Λ q Λ p − Λ r ω p − ω r + Λ p − Λ q ω p − ω r − Λ q − Λ q ω q − ω r ω p − ω r (cid:19) J − p J − q J − r | j, + j (cid:105) . (D4)Thus, we arrive at S − S − S − | j, j (cid:105) = 13! (cid:88) p,q,r =1 p,q,r distinct (cid:20) Λ p Λ q Λ r + 3Λ r Λ p − Λ q ω p − ω q + 2 ω p − ω q (cid:18) Λ p − Λ r ω p − ω r − Λ q − Λ r ω q − ω r (cid:19)(cid:21) J − p J − q J − r | j, j (cid:105) . (D5) [1] Y. Fukuda et al. (Super-Kamiokande), Phys. Rev. Lett. , 1562 (1998), arXiv:hep-ex/9807003 [hep-ex]. [2] Q. R. Ahmad et al. (SNO), Phys. Rev. Lett. , 071301 (2001), arXiv:nucl-ex/0106015 [nucl-ex].[3] F. P. An et al. (Daya Bay), Phys. Rev. Lett. , 171803(2012), arXiv:1203.1669 [hep-ex].[4] L. Wolfenstein, Phys. Rev. D , 2369 (1978).[5] S. P. Mikheyev and A. Y. Smirnov, Yad. Fiz. , 1441(1985).[6] W. C. Haxton, Physical Review Letters , 1271 (1986).[7] H. A. Bethe, Phys. Rev. Lett. , 1305 (1986),[,324(1986)].[8] G. M. Fuller, R. W. Mayle, J. R. Wilson, and D. N.Schramm, Astrophys. J. , 795 (1987).[9] D. N¨otzold and G. Raffelt, Nuclear Physics B , 924(1988).[10] J. T. Pantaleone, Phys. Lett. B287 , 128 (1992).[11] G. Sigl and G. Raffelt, Nuclear Physics B , 423(1993).[12] G. Raffelt, G. Sigl, and L. Stodolsky, Physical ReviewLetters , 2363 (1993), arXiv:hep-ph/9209276.[13] S. Samuel, Phys. Rev. D , 1462 (1993).[14] V. A. Kostelecky and S. Samuel, Phys. Rev. D49 , 1740(1994).[15] V. A. Kosteleck´y and S. Samuel, Phys. Rev. D , 621(1995), arXiv:hep-ph/9506262.[16] A. D. Dolgov, S. H. Hansen, S. Pastor, S. T. Petcov,G. G. Raffelt, and D. V. Semikoz, Nucl. Phys. B632 ,363 (2002), arXiv:hep-ph/0201287 [hep-ph].[17] C. M. Ho, D. Boyanovsky, and H. J. de Vega, Phys. Rev.
D72 , 085016 (2005), arXiv:hep-ph/0508294 [hep-ph].[18] B. H. J. McKellar and M. J. Thomson, Phys. Rev. D ,2710 (1994).[19] K. Enqvist, K. Kainulainen, and J. Maalampi, NuclearPhysics B , 754 (1991).[20] S. Pastor, G. Raffelt, and D. V. Semikoz, Phys. Rev. D , 053011 (2002), arXiv:hep-ph/0109035.[21] K. N. Abazajian, J. F. Beacom, and N. F. Bell, Phys.Rev. D , 013008 (2002), arXiv:astro-ph/0203442.[22] S. Hannestad, G. G. Raffelt, G. Sigl, and Y. Y. Y. Wong,Phys. Rev. D , 105010 (2006), arXiv:astro-ph/0608695.[23] L. Johns, M. Mina, V. Cirigliano, M. W. Paris,and G. M. Fuller, Phys. Rev. D94 , 083505 (2016),arXiv:1608.01336 [hep-ph].[24] H. Duan and J. P. Kneller, J. Phys.
G36 , 113201 (2009),arXiv:0904.0974 [astro-ph.HE].[25] H. Duan, G. M. Fuller, and Y.-Z. Qian, Annual Re-view of Nuclear and Particle Science , 569 (2010),arXiv:1001.2799 [hep-ph].[26] S. Chakraborty, R. Hansen, I. Izaguirre, and G. Raffelt,Nuclear Physics B , 366 (2016), arXiv:1602.02766[hep-ph].[27] A. B. Balantekin, J. Phys. G45 , 113001 (2018),arXiv:1809.02539 [hep-ph].[28] H. Duan, G. M. Fuller, J. Carlson, and Y.-Z. Qian, Phys.Rev. D , 105014 (2006), arXiv:astro-ph/0606616.[29] H. Duan, G. M. Fuller, J. Carlson, and Y.-Z. Qian,Physical Review Letters , 241101 (2006), arXiv:astro-ph/0608050.[30] H. Duan, G. M. Fuller, and Y.-Z. Qian, Phys. Rev. D , 123004 (2006), arXiv:astro-ph/0511275.[31] H. Duan, G. M. Fuller, and Y.-Z. Qian, Phys. Rev. D , 085013 (2007), arXiv:0706.4293.[32] H. Duan, G. M. Fuller, J. Carlson, and Y.-Z. Qian, Phys.Rev. D , 125005 (2007), arXiv:astro-ph/0703776.[33] H. Duan, G. M. Fuller, J. Carlson, and Y.-Z. Qian, Phys-ical Review Letters , 241802 (2007), arXiv:0707.0290. [34] H. Duan, G. M. Fuller, and Y.-Z. Qian, Phys. Rev. D , 085016 (2008), arXiv:0801.1363.[35] G. G. Raffelt and A. Y. Smirnov, Phys. Rev. D ,081301 (2007), arXiv:0705.1830.[36] G. G. Raffelt and A. Y. Smirnov, Phys. Rev. D ,125008 (2007).[37] G. Fogli, E. Lisi, A. Marrone, and A. Mirizzi, Journalof Cosmology and Astro-Particle Physics , 10 (2007),arXiv:0707.1998.[38] B. Dasgupta and A. Dighe, Phys. Rev. D77 , 113002(2008), arXiv:0712.3798 [hep-ph].[39] B. Dasgupta, A. Dighe, A. Mirizzi, and G. G. Raffelt,Phys. Rev. D , 113007 (2008), arXiv:0801.1660.[40] B. Dasgupta, A. Dighe, G. G. Raffelt, and A. Y.Smirnov, Physical Review Letters , 051105 (2009),arXiv:0904.3542.[41] A. Friedland, Physical Review Letters , 191102(2010), arXiv:1001.0996 [hep-ph].[42] B. Dasgupta, A. Mirizzi, I. Tamborra, and R. Tomas,Phys. Rev. D81 , 093008 (2010), arXiv:1002.2943 [hep-ph].[43] S. Galais and C. Volpe, Phys. Rev.
D84 , 085005 (2011),arXiv:1103.5302 [astro-ph.SR].[44] A. Mirizzi, Phys. Rev.
D88 , 073004 (2013),arXiv:1308.1402 [hep-ph].[45] S. Chakraborty and A. Mirizzi, Phys. Rev.
D90 , 033004(2014), arXiv:1308.5255 [hep-ph].[46] Y. Pehlivan, A. L. Suba¸sı, N. Ghazanfari, S. Birol,and H. Y¨uksel, Phys. Rev.
D95 , 063022 (2017),arXiv:1603.06360 [astro-ph.HE].[47] J. Y. Tian, A. V. Patwardhan, and G. M. Fuller,Phys. Rev.
D96 , 043001 (2017), arXiv:1703.03039 [astro-ph.HE].[48] A. Das, A. Dighe, and M. Sen, JCAP , 051 (2017),arXiv:1705.00468 [hep-ph].[49] S. Birol, Y. Pehlivan, A. B. Balantekin, and T. Kajino,Phys. Rev.
D98 , 083002 (2018), arXiv:1805.11767 [astro-ph.HE].[50] V. Cirigliano, M. Paris, and S. Shalgar, JCAP ,019 (2018), arXiv:1807.07070 [hep-ph].[51] A. Malkus, J. P. Kneller, G. C. McLaughlin, and R. Sur-man, Phys. Rev. D , 085015 (2012), arXiv:1207.6648[hep-ph].[52] A. Malkus, A. Friedland, and G. C. McLaughlin, ArXive-prints (2014), arXiv:1403.5797 [hep-ph].[53] M.-R. Wu, H. Duan, and Y.-Z. Qian, Physics Letters B , 89 (2016), arXiv:1509.08975 [hep-ph].[54] C. J. Stapleford, D. J. V¨a¨an¨anen, J. P. Kneller, G. C.McLaughlin, and B. T. Shapiro, Phys. Rev. D94 , 093007(2016), arXiv:1605.04903 [hep-ph].[55] A. Malkus, G. C. McLaughlin, and R. Surman, Phys.Rev. D , 045021 (2016), arXiv:1507.00946 [hep-ph].[56] M. Frensel, M.-R. Wu, C. Volpe, and A. Perego,Phys. Rev. D95 , 023011 (2017), arXiv:1607.05938 [astro-ph.HE].[57] D. V¨a¨an¨anen and G. C. McLaughlin, Phys. Rev. D ,105044 (2016), arXiv:1510.00751 [hep-ph].[58] Y.-L. Zhu, A. Perego, and G. C. McLaughlin, Phys. Rev. D94 , 105006 (2016), arXiv:1607.04671 [hep-ph].[59] S. Shalgar, JCAP , 010 (2018), arXiv:1707.07692[hep-ph].[60] A. Vlasenko and G. C. McLaughlin, Phys. Rev.
D97 ,083011 (2018), arXiv:1801.07813 [astro-ph.HE].[61] R. F. Sawyer, Phys. Rev. D , 045003 (2005), arXiv:hep- ph/0503013.[62] S. Chakraborty, R. S. Hansen, I. Izaguirre, and G. Raf-felt, JCAP , 042 (2016), arXiv:1602.00698 [hep-ph].[63] F. Capozzi, B. Dasgupta, and A. Mirizzi, JCAP ,043 (2016), arXiv:1603.03288 [hep-ph].[64] M.-R. Wu and I. Tamborra, Phys. Rev. D95 , 103007(2017), arXiv:1701.06580 [astro-ph.HE].[65] M. Sen,
Proceedings, XXII DAE High Energy PhysicsSymposium: Delhi, India, December 12 -16, 2016 ,Springer Proc. Phys. , 533 (2018), arXiv:1702.06836[hep-ph].[66] A. Dighe and M. Sen, Phys. Rev.
D97 , 043011 (2018),arXiv:1709.06858 [hep-ph].[67] S. Abbar and H. Duan, Phys. Rev.
D98 , 043014 (2018),arXiv:1712.07013 [hep-ph].[68] M.-R. Wu, I. Tamborra, O. Just, and H.-T. Janka,Phys. Rev.
D96 , 123015 (2017), arXiv:1711.00477 [astro-ph.HE].[69] B. Dasgupta and M. Sen, Phys. Rev.
D97 , 023017 (2018),arXiv:1709.08671 [hep-ph].[70] B. Dasgupta, A. Mirizzi, and M. Sen, Phys. Rev.
D98 ,103001 (2018), arXiv:1807.03322 [hep-ph].[71] F. Capozzi, B. Dasgupta, A. Mirizzi, M. Sen, and G. Sigl,Phys. Rev. Lett. , 091101 (2019), arXiv:1808.06618[hep-ph].[72] S. Abbar and M. C. Volpe, Phys. Lett.
B790 , 545 (2019),arXiv:1811.04215 [astro-ph.HE].[73] S. Shalgar and I. Tamborra, (2019), arXiv:1904.07236[astro-ph.HE].[74] N. F. Bell, A. A. Rawlinson, and R. F. Sawyer, Phys.Lett.
B573 , 86 (2003), arXiv:hep-ph/0304082 [hep-ph].[75] A. Friedland and C. Lunardini, Phys. Rev.
D68 , 013007(2003), arXiv:hep-ph/0304055 [hep-ph].[76] A. Friedland and C. Lunardini, JHEP , 043 (2003), arXiv:hep-ph/0307140 [hep-ph].[77] A. B. Balantekin and Y. Pehlivan, Journal of Physics GNuclear Physics , 47 (2007), arXiv:astro-ph/0607527.[78] Y. Pehlivan, A. B. Balantekin, T. Kajino, andT. Yoshida, Phys. Rev. D84 , 065008 (2011),arXiv:1105.1182 [astro-ph.CO].[79] M. J. Cervia, A. V. Patwardhan, and A. B. Balantekin,(2019), arXiv:1905.00082 [nucl-th].[80] F. Pan, J. P. Draayer, and W. E. Ormand, Phys. Lett.
B422 , 1 (1998), arXiv:nucl-th/9709036 [nucl-th].[81] A. B. Balantekin, J. H. de Jesus, and Y. Pehlivan, Phys.Rev.
C75 , 064304 (2007), arXiv:nucl-th/0702059 [nucl-th].[82] R. W. Richardson, Phys. Rev. , 949 (1966).[83] M. Gaudin, J. Physique , 1087 (1976).[84] A. Faribault, O. E. Araby, C. Strater, and V. Gritsev,Phys. Rev. B83 , 235124 (2011), arXiv:1103.0472 [cond-mat.mes-hall].[85] O. Babelon and D. Talalaev, Journal of Statistical Me-chanics: Theory and Experiment , P06013 (2007).[86] N. Bourbaki,
Algebra I: Chapters 1-3 , Actualit´es scien-tifiques et industrielles (Springer, 1998).[87] P. W. Claeys, D. V. Neck, and S. D. Baerdemacker,SciPost Phys. , 028 (2017).[88] P. W. Claeys, Richardson-Gaudin models and brokenintegrability , Ph.D. thesis, Ghent University (2018),arXiv:1809.04447 [math-ph].[89] C. B. Garcia and W. I. Zangwill, Mathematical Program-ming , 159 (1979).[90] O. E. Araby, V. Gritsev, and A. Faribault, Phys. Rev. B85 , 115130 (2012), arXiv:1201.5542 [cond-mat.mes-hall].[91] P. W. Claeys, S. De Baerdemacker, M. Van Raemdonck,and D. Van Neck, Phys. Rev. B91