Slow and fast collective neutrino oscillations: Invariants and reciprocity

The flavor evolution of a neutrino gas can show ''slow'' or ''fast'' collective motion. In terms of the usual Bloch vectors to describe the mean-field density matrices of a homogeneous neutrino gas, the slow two-flavor equations of motion (EOMs) are $\dot{\mathbf{P}}_\omega=(\omega\mathbf{B}+\mu\mathbf{P})\times\mathbf{P}_\omega$, where $\omega=\Delta m^2/2E$, $\mu=\sqrt{2} G_{\mathrm{F}} (n_\nu+n_{\bar\nu})$, $\mathbf{B}$ is a unit vector in the mass direction in flavor space, and $\mathbf{P}=\int d\omega\,\mathbf{P}_\omega$. For an axisymmetric angle distribution, the fast EOMs are $\dot{\mathbf{D}}_v=\mu(\mathbf{D}_0-v\mathbf{D}_1)\times\mathbf{D}_v$, where $\mathbf{D}_v$ is the Bloch vector for lepton number, $v=\cos\theta$ is the velocity along the symmetry axis, $\mathbf{D}_0=\int dv\,\mathbf{D}_v$, and $\mathbf{D}_1=\int dv\,v\mathbf{D}_v$. We discuss similarities and differences between these generic cases. Both systems can have pendulum-like instabilities (soliton solutions), both have similar Gaudin invariants, and both are integrable in the classical and quantum case. Describing fast oscillations in a frame comoving with $\mathbf{D}_1$ (which itself may execute pendulum-like motions) leads to transformed EOMs that are equivalent to an abstract slow system. These conclusions carry over to three flavors.


I. INTRODUCTION
The refractive energy shift of a neutrino in a medium of neutrinos or antineutrinos of its own flavor is twice that for a different flavor. In a seminal paper thirty years ago, Pantaleone showed that therefore neutrino-neutrino refraction is a many-body phenomenon [1]. Shortly thereafter, in another foundational paper, Samuel found that neutrino flavor evolution feeding back on itself spawns intriguing forms of collective flavor evolution [2]. While the community was somewhat slow at fully catching on to these insights, the idea is now standard that neutrinoneutrino refraction can strongly affect flavor-dependent neutrino transport, notably in supernova cores or other (neutrino-)dense environments.
Some forty research papers on this topic have appeared in 2022 alone, making it impossible to do justice to this swelling body of literature in our brief introduction. So we merely refer to some early [3,4] and more recent reviews [5][6][7][8] that provide a glimpse of the richness of questions addressed in this quest. Much recent work is directed toward understanding space-time dependent solutions in the context of realistic astrophysical environments, including effects of collisions, or understanding possible limitations of the mean-field description.
We assume that the reader is largely familiar with this general subject and here return to far more elementary questions about the basic structure of the nonlinear flavor equations of motion (EOMs). In particular, we will juxtapose the basic reference case of "slow" (energydependent) flavor oscillations with the "fast" (angledependent) case and discuss a number of striking similarities and a certain reciprocity between them.
One key element of collective flavor evolution is the possible appearance of instabilities in the linearized equations, signifying soliton solutions (pendulum-like behavior) in the nonlinear regime. These were identified a long time ago in the simplest collective-oscillation case [9][10][11][12][13], which is a homogeneous and isotropic (or single-angle) neutrino gas with a nontrivial energy distribution. A similar behavior was recently identified in the simplest system of fast flavor conversion [14,15] that consists of a homogeneous and axisymmetric neutrino gas with a nontrivial angle distribution of the lepton-number flux. We here juxtapose these generic cases more systematically, discuss similarities and differences, and show a certain reciprocity between them.
To define the exact cases of comparison, we use the usual Bloch vectors to represent two-flavor density matrices. The simplest collective EOMs arė where ω = ∆m 2 /2E is the vacuum oscillation frequency, µ = √ 2G F (n ν +nν) is a measure of the neutrino-neutrino interaction strength, and P = dω P ω . Antineutrinos are included with negative ω and we use the flavor isospin convention: spin up means ν e orν µ , spin down means ν µ orν e . We usually assume |ω| ≪ µ, a condition that defines a "neutrino-dense" environment. Because flavor conversion is here driven by ω, this case has come to be called "slow flavor oscillations." If we interpret the Bloch vectors as classical angular momenta, these EOMs follow from the Hamiltonian 1 where P 0 = P = dω P ω and P 1 = dω ωP ω . The second case of comparison uses Bloch vectors for the density matrices of neutrinos minus that for antineutrinos, i.e., for the angle-dependent lepton-number occupation. Here spin up means ν e −ν µ occupation, spin down ν µ −ν e occupation. Assuming axial symmetry and using the fast flavor limit of vanishing neutrino mass splitting (ω = 0), the EOMs arė where v = cos θ is the mode velocity along the symmetry axis. In analogy to the slow case, we use the moments 2 The fast EOMs derive from the Hamiltonian once more using the canonical angular-momentum Poisson brackets for the D v . Our aim is to compare the properties of these minimal slow and fast cases. The slow Hamiltonian, or its quantum equivalent, arises in many other "spin problems" and in particular in the context of the BCS theory of superconductivity. Using Anderson's pseudo-spin formalism [17], a Cooper pair means "spin down" whereas unpaired electrons correspond to "spin up," implying a perfect analogy. For example, pendulum-like behavior was discovered in this system, collective oscillation between the paired and unpaired state, in a nonequilibrium situation when dissipation was small [18,19]. Many features of what we call collective neutrino oscillations were discussed in a long series of papers by Yuzbashyan and collaborators [20][21][22][23][24][25], where often the classical and quantum cases are clearly juxtaposed. We consider mostly the classical case (the mean-field description of neutrino oscillations) and the EOMs are central to our arguments rather than the Hamiltonians themselves. Indeed, we argue the reciprocity of the slow and fast system on the level of the EOMs, not on the basis of canonical transformations between the variables.
Returning to our discussion, as a final simplification we absorb the scale µ in the definition of dimensionless time and use dimensionless oscillation frequencies w = ω/µ. Then finally we consider the two sets of EOMṡ where −∞ < w < +∞ and −1 ≤ v ≤ +1. We will continue to denote the Bloch vectors of the slow system as P w and those of the fast system as D v . These equations look vaguely similar and it is the main point of our paper to show in which sense the two system are actually reciprocal to each other.
We begin our discussion in Sec. II with the mapping between the fast and an abstract slow system and we identify and interpret the constraint for the slow system required to enable the reverse mapping. In Sec. III we develop the machinery of Lax vectors to identify the invariants and prove the integrability of both systems in the mechanical sense. Moreover, we discuss the equivalence to a system with only few effective degrees of freedom based on the roots of the spectral polynomial. In Sec. IV we juxtapose the dispersion relations for the fast and slow cases and show the consistency of the original mapping on the level of the linearized system. We also explain the equivalence of the normal modes to the roots of the z-component of the Lax vectors if the neutrinos begin in flavor eigenstates. Next, in Sec. V, we extend the discussion to three flavors and derive the corresponding Gaudin invariants. We end with some general conclusions in Sec. VI. We address a number of additional topics in appendices. In App. A we explain the general transformation between moving frames in flavor space. In App. B we discuss the dispersion relation in the continuum (thermodynamic) limit following earlier work on the Vlasov equation. Finally, in App. C, we prove the integrability of the quantum version of the fast system.

A. Transforming a General Fast System to a Constrained Slow System
To compare Eqs. (6) and (7) we first identify a few conserved quantities. Integrating Eq. (7) over dv leads toḊ 0 = 0 and thus reveals that D 0 is a fixed vector.
Applying v dv on both sides provideṡ The evolution of D 1 is an instantaneous precession and so its length is conserved. For those conditions where the motion is analogous to a gyroscopic pendulum, D 0 plays the role of gravity, whereas D 1 that of the moving radius vector that executes precessions and nutations. There is an evident reciprocity between D 1 moving relative to the fixed D 0 or conversely D 0 moving in a frame comoving with D 1 . To express the fast oscillations of Eq. (7) in this comoving frame we observe that Eqs. (7) and (8) both engender differential rotations that can be added. Therefore, we can simply subtract the precession of the frame comoving with D 1 from the original precession of the individual D v . Denoting Bloch vectors in the moving frame with a tilde, we thus find We present a more formal derivation of this transformation in Appendix A. Applying v dv on both sides confirms that ∂ tD1 = 0. Moreover, applying dv on both sides reveals that ∂ tD0 = −D 2 ×D 0 and so it is indeed D 0 that now is a moving object with conserved length. All vectors coincide between the two frames at t = 0, so the spectra, stability conditions, and so forth are the same in both frames. In the fast EOMs (7), the common precession around the fixed vector D 0 does not affect the internal motion of the system. One usually assumes that initially all Bloch vectors are oriented in the flavor direction, apart from a small seed that can trigger possible instabilities, i.e., a small deviation of collinearity among the D v . We still assume that the conserved D 0 defines the z-direction. The flavor content of every mode is encoded in D z v which does not depend on the common precession. Therefore, the questions of physical relevance are the same in the co-rotating system so that one can simply remove the first term of Eq. (7) and finds the same physical answers. In this case, D 0 disappears from both Eqs. (7) and (8) and thus never appears in the transformed Eq. (9). Of course, it is physically clear that the internal motion of the D v relative to D 1 does not depend on the common precession of all D v and of So far we have treated w and v as continuous parameters. However, sometimes it is more practical to think of discrete sets of Bloch vectors to avoid worrying about integration measures. In this spirit integrals dv or dw correspond to v or w . With this understanding, the transformed EOMs (9) are identical in form to the slow EOMs (6) with the identifications where D 1 = |D 1 | = |D 1 |, −D 1 ≤ w ≤ +D 1 , and These identifications are possible because in the new frame, the first momentD 1 is constant and acts as a fixed direction around which all polarization vectors precess with different frequencies.
B. Constrained Slow System The slow system has a vector of conserved length that can play the role of a pendulum vector, 3 3 The integral in the continuous case would be interpreted in the principal value sense, provided Pw is integrable at w = 0. Notice, however, that w corresponds to infinite E and in a realistic medium would be exponentially suppressed. Moreover, Eq. (10) contains a factor w 2 times a slowly-varying function.
indeed R follows an instantaneous precession. Actually the Bloch vector R thus defined is a special case of a Lax vector, a topic that we will discuss later in Sec. III A. However, the transformed variables of Eq. (10) reveal This case of a slow flavor pendulum has not been discussed in the literature, but of course Eq. (12) always allows for the special case R = 0. Notice that, if R = 0 at the initial instant, it remains zero throughout the evolution. In our case, after transforming the fast EOMs, we have already shown that it isD 0 that is the moving quantity of conserved length.
In the slow case with R = 0, whether or not it derives from the fast-slow transformation, we should thus Because by assumption P −1 = B, the first term drops out and so indeedṖ −2 = P 0 ×P −2 is a precession equation that conserves the length of P −2 .
We can also reverse the transformation and start from the slow EOMs (6) and go to a frame that moves with the vector that stays of constant length. This is either R or P −2 , but in both cases the evolution is driven by P 0 . Following the same logic as in the earlier transformation, the moving EOMs are ∂ tPw = wB×P w . Next we can use Eq. (12) to eliminateB and find with the identification v = w and D v = −P w /w 2 the equivalent EOMṡ This has the form of a fast system, but only if R = 0 of the original slow system. Therefore, in both directions of the transformation one finds that the fast EOMs are equivalent to a specific slow system where R = 0.
The meaning of this constraint becomes more obvious for a system of N discrete modes D v in the fast system. After the transformation we obtain N slow modes P w and some combination of D v playing the role of the fixed vector B. The additional degree of freedom is absorbed by the constraint R = 0 so that the slow system also has N Bloch-vector degrees of freedom.
Therefore, we find that a slow system can be mapped to the same form as a fast system if the vector R = 0. However, more generally, we can transform the EOMs to a frame co-rotating around the z axis with frequency w c , where the EOMs arė and perform the same steps as before. In this frame, we need to require that This expression has the same form as a Lax vector Eq. (17) of this system. In other words, the slow-to-fast mapping can be done if there exists any vanishing Lax vector L u with real u.

C. Conserved Energy
The slow and fast EOMs derive respectively from the classical Hamiltonians Eqs. (2) and (5) if we interpret the Bloch vectors as classical angular momenta. The two terms in H fast are separately conserved and indeed, the same motion obtains if we leave out the first piece 1 2 D 2 0 , except for an overall precession around D 0 . In the slow case, on the other hand, the two pieces exchange energy in the course of the motion. For the pendulum solutions, the first term plays the role of potential energy, the second term the role of kinetic energy.
After mapping the fast system on a slow one, we recall that B ∼ D 1 , P 1 ∼ D 3 , and P 0 ∼ D 2 . Therefore one expects that the quantity is conserved as one can verify using the fast EOMs (7). This is one of many invariants of the fast system.

A. Lax Vectors
The slow EOMs (6) admit an infinite set of constants of the motion which are best expressed in terms of the so-called Lax vectors [22,24] where in general u is a complex parameter, whereas w is understood as a set of discrete frequencies. The time derivative is found by inserting Eq. (6) forṖ w on the right-hand side, leading tȯ and where P 0 = w P w . This is a formal precession equation with complex frequency u and implies that the length of L u is conserved, i.e., ∂ t L 2 u = 0, where L 2 u is a complex number. For the fast system, we find that the corresponding Lax vectors are In analogy, the time derivative is found by inserting Eq. (7) forḊ v on the right-hand side, leading tȯ with D 0 = v D v and D 1 = v vD v . So here as well, the complex length of L u is conserved for any complex u.

Slow System
The Lax vectors provide a continuous infinity of invariants of the motion, which however are not linearly independent. However, considering the Lax vectors only for the same discrete set w of frequencies, then there is a real Lax vector L w for every P w . It provides the Gaudin invariant or Gaudin Hamiltonian [26] This is a conserved quantity of our original system, but can also be interpreted as a Hamiltonian that defines a new system of a "central" spin P w coupled to an external B field and to a set of "environmental" spins, the "central spin problem" (see e.g. Ref. [20]). The Gaudin invariants generate all the other integrals of motion. One example is as can be seen immediately from the antisymmetry of the second term in Eq. (21) under exchange w ↔ w ′ . Likewise, the Hamiltonian of the original slow system can be expressed as 4 Notice that w P 2 w is separately conserved and thus could be left out without changing the EOMs. This term is included to make up for the piece w = w ′ that is missing in sums involving 1/(w − w ′ ). 4 This is seen if we consider where we have interchanged the summation variables in the last expression, leading to the minus sign. Taking this term to the left-hand side finally reveals that the original left-hand side is simply 1 A generalization of these invariants is given by the ex- where C 0 was already given in Eq. (22) and C 1 is H slow given in Eq. (23).

Fast System
To find the Gaudin invariants of the fast system, once more we consider the Lax vectors L v for the same set of discrete velocities v as the original fast system. Once more they follow the same EOMs as the original set D v and so we find the Gaudin invariants A trivial conserved quantity is The next invariant vanishes identically. Further nontrivial conserved quantities are obtained for n ≥ 2 In terms of this chain of invariants, the fast Hamiltonian Eq. (5) is On the other hand, the Hamiltonian in the frame comoving with D 1 that was given in Eq. (16) is the same as C 4 . The invariants of the system can be written in many ways, each providing a different glance on its properties. 5 For additional expressions see Appendix A of Ref. [16]. Their Eq. (A3) correspond to our Eq. (24a) with their µ → 1/2.

C. How Many Invariants per Mode?
We have identified two invariants for every mode, the length of the Bloch vector and the Gaudin invariant. The latter is the scalar product of the Bloch vector and the corresponding Lax vector. So for the fast system, we have D 2 v and D v · L v as invariants, but in addition L 2 v is also conserved because the Lax vector also follows a precession equation. While it is physically obvious that the system cannot have another independent invariant per mode, it is instructive to show this explicitly. Following the definition of the Lax vector in Eq. (19), its square is We modify the right-hand side with the identity so that Therefore, the set of invariants L 2 v is indeed a linear combination of the Gaudin invariants I v .

D. Integrability
The existence of two constants of motion for each of the Bloch vectors makes the system completely integrable in the mechanical sense. First of all, the nature of the motion as a precession (or in matrix form, the commutator structure) reveals that the length of each Bloch vector is conserved, reducing the cartesian variables to two polar coordinates. Therefore, the motion involves periodic variables. In addition, there is another conserved quantity for each P w given, for example, by the Gaudin invariants. (The nature of the fast case as a special case of the slow one makes redundant a separate discussion of the former.) Therefore, the system has only as many underlying variables as independent Bloch vectors, i.e., the number of different discrete w. Indeed, one can introduce canonical variables which separate the Hamilton-Jacobi equations as those frequencies ω * for which the vectors L ω * are along the z-axis, and their canonical conjugate variables, namely the z component of L ω * . An explicit solution for any initial condition can then be constructed following, e.g., Refs. [20,27].
When the system is quantized, the integrability of the equations shows up in the possibility of constructing all of the eigenstates of the Hamiltonian analytically by the Bethe ansatz. (See Appendix C for the integrability of the quantum system.)

E. Degenerate Solutions
The large number of conserved quantities implies that the motion in the 2N tot -dimensional phase space -N tot being the number of Bloch vectors -is restricted to an N tot -dimensional torus. However, depending on the set of frequencies ω i and initial P i , the motion of the system can be constrained to a reduced phase space instead of ergodically filling all of the torus. Such solutions were called "degenerate" in Ref. [20], whereas in a previous study by one of us [13] they were dubbed N-mode coherent solutions. Either way, the solutions P i (t) with i = 1, . . . , N tot would be linear combinations of a smaller set of functions J i (t) with i = 1, . . . , N < N tot and in this sense the full set P i (t) moves coherently rather than each of them independently. The J i (t) were termed "auxiliary spins" in Ref. [20] or "carrier modes" in Ref. [13].
To show that this phenomenon can indeed exist, we first reverse the problem and introduce a fictitious system of a smaller number N < N tot of Bloch vectors J i with frequencies χ i obeying an analogous slow EOM with the same unit vector Ḃ where We here return to dimensionful frequencies and coupling strength µ.) This is yet another slow system and as such integrable.
Therefore, we can construct Lax vectors on an arbitrary set of frequencies ω i with i = 1, . . . , N tot > N and suppose that none of the chosen ω i coincides with any of the χ i . This set of Lax vectors of the fictitious system is and as such obeys the EOMṡ Let us assume that the {L i , ω i } are fixed such that the matching condition L = J is obeyed. This can be achieved by another factor depending on i that adjusts the arbitrary length of L i and/or the choice of ω i . The new system L i is a higher-dimensional slow system with solutions L i (t) that depend on B and N < N tot independent functions J i (t) with i = 1, . . . , N . More difficult is the reverse problem where we begin with a given system {P i , ω i } with i = 1, . . . , N tot and ask for its true dimensionality. One numerical diagnostic requires to solve the EOMs for some chosen time interval and determine the components of the Gram matrix, [13]. Its rank (the number of independent eigenvalues) reveals the number of linearly independent underlying functions.
However, this dimensionality can also be directly determined from the initial conditions without solving the EOMs [20]. The key information is encoded in the complex length of the Lax vector L 2 u . If it vanishes for a real u at the initial time, it must vanish at all times, which means that one of the polarization vectors can be expressed in terms of the others, and therefore the number of degrees of freedom drops by one. Indeed, the root diagram of L 2 u , the generally complex solutions u i of L 2 u = 0, provides a lot of crucial information about a given system [20]. Formally one defines the "spectral polynomial" of our system by Here, p u is the product of the denominators of the Lax vectors and so, Q u is the numerator of L 2 u and as such a polynomial in u of order 2N tot . For all real u, we have Q u ≥ 0, which implies that the 2N tot zeroes are pairs of complex conjugate solutions. As a consequence, real zeros come in pairs of degenerate solutions.
In the fast flavor case, the constant vector B is absent from the Lax vector Eq. (19), but on the other hand there is a factor u in the numerator, so the spectral polynomial once more has degree 2N tot if there are N tot Bloch vectors D v . However, one trivial real double root is now u = 0, so there are at most N tot − 1 independent degrees of freedom, whereas one is eliminated, corresponding to the conserved D 0 .
In summary, the number N of pairs of complex conjugate roots of the spectral polynomial is the number of independent underlying degrees of freedom. For an arbitrary set {P i , ω i }, very special initial conditions are required to obtain a reduced number of degrees of freedom. In the continuum of possible initial conditions, only a discrete set of measure zero fulfills this condition [20]. However, in the neutrino context we consider very special initial conditions, i.e., all Bloch vectors are initially collinear with the flavor direction. In this case, the situation is reversed in that there are very few independent degrees of freedom, providing soliton solutions.

IV. NORMAL MODES AND LAX VECTORS
In the context of dense neutrino environments, one usually imagines that they were produced in flavor eigenstates and that the in-medium mixing angles are small. So for both slow and and fast oscillations, one typically studies initial conditions where all Bloch vectors are collinear with the flavor direction that we identify with the z-direction in flavor space. Analogous assumptions were made in the cited discussions of the BCS Hamiltonian where the initial state consists of paired or unpaired electrons, but not initially of coherent superpositions of paired with unpaired states [20][21][22][23][24][25].
Any evolution away from this initial state begins with small excursions from the z-directions, lending itself to a linearized analysis and to the identification of possible run-away modes that can lead to large flavor-conversion effects. Eventually these modes become nonlinear, representing what were called "pendulum solutions" in the neutrino literature or "solitons" in the BCS literature. We juxtapose the fast and slow cases in the context of a linear normal-mode analysis and we highlight the connection in both cases between the linear normal-mode analysis and the Lax-vector analysis.
A. Normal-Mode Analysis

Slow System
The fast and slow systems can both be analyzed in the linear regime of small off-diagonal components of the density matrices. This approach makes sense in the usual picture where all Bloch vectors start essentially collinear except for a small seed to trigger run-away modes. For both systems, the starting point are the "spectra" that define the systems. This includes both continuous functions or collections of discrete modes that are represented by collections of δ functions at different values of w or v. In this case, integrations over the spectra are effectively summations over discrete modes. The most general space-time dependent linearization and concomitant dispersion relation was developed in Ref. [28], but here it is somewhat more transparent to linearize our simple systems directly. For the slow case we may follow the first study of linear stability analysis [29] and note that initially P xy w ≡ P x w − iP y w is small, whereas P z w = g w is at its starting value. 6 Then the linearized version of Eq. (6) is where P z 0 = dw g w . We have explicitly included B z for later comparison with D z 1 . While B is defined as a unit vector collinear with the z-direction, both B z = ±1 are in principle possible.
A collective normal mode is of the form P xy where Ω is the eigenfrequency that can be complex. If P w with length |g w | is expressed in polar coordinates (ϑ w , ϕ w ), we have Q w = sin ϑ w e −iϕw . The linear EOMs imply the eigenfunction 6 From the relation to the density matrix ̺w = 1 2 Pw ·σ we observe that P xy w is the "upper right" entry of ̺w. We could also use the lower-left off-diagonal element, the complex conjugate of P xy w , leading however to the same information.
where A does not depend on w. Inserting this result back into the linear EOM reveals that the eigenfrequency is fixed by the condition Recall that while P 0 moves, its projection P z 0 = dw g w on B is conserved.
If the spectrum consists of n discrete modes, one finds n solutions that can be real or pairs of complex conjugate ones. If g w is continuous and nonzero for all w, except perhaps some crossings, the poles in the integral prevent any real solution of the "dispersion relation" Eq. (40). The singularity associated with the vanishing of the denominator in Eq. (40) requires special care in the interpretation, although it leads to no practical difficulties for the special problem we have at hand, namely identifying the instabilities of the system. We discuss this issue in Appendix B. For the discussion in the main text, we limit only to complex frequencies, for which there is no pole in the integration; this is enough to understand whether there are instabilities.
If Ω = ω P + iΓ is complex, the condition of Eq. (40) consists of a real and imaginary part that are explicitly The second equation implies that if there is any solution with Γ = 0, then there are two unstable solutions ±Γ, i.e., an exponentially growing and shrinking one, a point that is of course obvious in the discrete case anyway. Moreover, for Γ = 0, the second equation can only be true if g w has regions of positive and negative sign, an argument that also applies to both the discrete and continuum case. Unstable solutions require a "spectral crossing" as a necessary condition, which is also sufficient in the more general case of inhomogeneous solutions [30,31].

Fast System
For the fast case, we follow the same steps or use directly the explicit derivation in Ref. [15]. Here the linear EOMs are with the additional constraint v D xy v = 0 because the z-direction is defined to be the direction along D 0 . We will return to this question below. The eigenfunction is providing the condition where D z 0 = dvG v and D z 1 = dvvG v . Notice that D z 0 only shifts the real part of Ω and thus has no impact on the stability and was left out in Ref. [15]. In this sense one can use a co-rotating frame where D 0 drops out from the original EOMs. However, to avoid jumping between too many moving frames, we here keep all terms explicitly.
To compare Eqs. (40) and (44) we recall that D 1 = |D 1 | = |D z 1 | is positive, whereas D z 1 = dv G v can be both positive or negative. If we absorb D z 0 and P z 0 in the real part of Ω, Eq. (40) follows from Eq. (44) with the substitution D z , meaning g(w) = −(w/D 1 ) 2 G(w/D 1 )/D 1 , and the limits of integration are ±D 1 , or alternatively one can say that g(w) = 0 for |w| > D 1 . These identifications correspond to Eqs. (10) and (11), but now for a continuous system that involves a factor D 1 as a Jacobian from the v → w transformation.
Equation (44) is trivially solved by the real solution Ω = D z 0 . The corresponding eigenvector Q v from Eq. (43) is independent of v, and therefore seems to represent a uniform rotation of all polarization vectors, which however we had excluded in the beginning by the condition dvD xy v = 0. However, the Ω = D z 0 is not spurious, but actually can be constructed explicitly. Notice that for Ω = D z 0 , the eigenvector Q v can possess an arbitrary contribution from the Bloch vector with v = 0, which automatically satisfies the linearized EOM (42). Therefore, the eigenvector for the mode Ω = D z 0 is which correctly satisfies the constraint dvD xy v = 0. Physically, this mode corresponds to all Bloch vectors aligned along a direction slightly tilted from the z-axis, while the v = 0 vector is tilted in the opposite way so that D 0 is still along the z-axis. Since D v=0 does not contribute to D 1 (and in fact all D n with n > 0), the latter is still aligned with all Bloch vectors with v = 0, which therefore just precess around D 0 , whereas the polarization vector with v = 0 automatically satisfies the same pure precession motion sinceḊ v=0 = D 0 × D v=0 . Therefore, the fast system always admits a stable, uniformly rotating mode, in which all the moments D n with n > 0 are slightly tilted from the z axis. The existence of this mode is crucially connected with the conservation of D 0 .
By that same token, any function g w derived from G v as a fast-to-slow conversion fulfills where the sign depends on B z . On the other hand, starting with a generic slow system with any function g w , normally it will not fulfill this constraint which indeed corresponds exactly to the constraint derived in the previous section that the vector R vanishes, which is a necessary condition for a slow system to be mapped to a fast one. We now understand the physical origin of this issue: a fast system always admits a conserved vector D 0 because of its rotational invariance, and in turn it admits a uniformly precessing mode. A slow system can only be mapped to a fast system if it accidentally also admits a conserved vector. In the slow system case, this conservation is not protected by any fundamental symmetry, and therefore it only arises if the condition R = 0 is accidentally satisfied. For later reference, we finish by providing an alternative form of the fast eigenvalue equation. To this end, we first absorb D z 0 in the definition of Ω in Eq. (44). On the right-hand side, we use the manipulation providing the alternative form where we have restored the original meaning of Ω.

B. Sign of Spectral Crossing
A complex solution for Ω requires a spectral crossing as discussed earlier. In the slow case, it is well known that in addition the crossing must be positive, meaning that g w as a function of w must cross from negative to positive values of g w [32]. This condition applies for B oriented in the positive z-direction or rather, for B defining what is the positive z-direction. The reason for this condition is that the potential energy of the slow Hamiltonian in Eq. (2) is already minimal for a negative crossing and cannot be lowered by the P w moving away from their initial position. Of course, a system with a positive crossing need not be unstable-it could be stuck in the "sleeping top" position.
Therefore, also in the fast system the crossing should be positive measured relative to the direction defined by D z 1 = dv G v . For the explicit examples constructed in Ref. [15], this is indeed the case. While the crossings of G w shown in their Fig. 1 look negative, also D z 1 < 0. So a positive crossing in this sense is a necessary condition, but also in the fast case not a sufficient one, as seen in example A of that paper.

C. Normal Modes from Lax Vectors
The eigenvalue equation for normal modes presented earlier can be derived in a much simpler way using the Lax vectors of a given system. We focus on the example of the fast system, but the proof proceeds identically for the slow system as well. The Lax vector defined in Eq. (19) obeys Eq. (20) or explicitlẏ where we use the compact notation h u ≡ D 0 −uD 1 . This means in particular thaṫ where L ± u = L x u ± iL y u . Notice that in general both L x u and L y u are complex numbers. Let us assume that initially all polarization vectors are aligned with the z axis, except for a small perturbation δD v orthogonal to them. We denote by δ all perturbed quantities. If in the unperturbed situation we can find some u such that the condition holds, then it follows that the first term on the right-hand side of Eq. (50) drops out and the remaining motion is a rotation around the z axis with frequency h z u = D z 0 −uD z 1 , i.e., the eigenfrequency is Ω = D z 0 − uD z 1 that can be both real or complex. The condition Eq. (51) is recognized to be equivalent (in integral form) to the condition Eq. (48). It also follows that the corresponding eigenvector is aligned exactly with the Lax vector, namely that if the initial perturbation is chosen as where e i is the unit vector along the axis i, it will evolve in time proportionally to itself with a frequency Ω = D z 0 − uD z 1 . By replacing u = (D z 0 − Ω)/D z 1 , we recover the eigenstates obtained by linear dispersion analysis in Eq. (43).
Notice that this derivation closely parallels the derivation of the eigenfrequencies of the quantum Hamiltonian in Appendix C.

V. THREE-FLAVOR EXTENSION
The three-flavor case is notoriously more complicated, where the mean field of neutrino flavor is represented by 3 × 3 density matrices or by 8-dim Bloch vectors. Yet the technique of Lax vectors carries through in this case as well and one may identify similar Gaudin invariants and prove the integrability of the homogeneous EOMs.

A. Bloch Vectors and Equations of Motion
In the three-flavor case, the homogeneous flavor mean field is represented by 3 × 3 density matrices. For the axially symmetric fast flavor case each mode v is represented by ̺ v . It can be associated to an 8-dimensional Bloch vector using the decomposition [33,34] where λ i are the 8 Gell-Mann matrices. Here ̺ v is already written as a traceless matrix, because the trace of the matrix is the total number of neutrinos which is conserved in the absence of collisions. Therefore, one can still describe the dynamics in terms of vector equations of motion, although in 8-dimensional space. We recollect the main properties of the Gell-Mann matrices which will be useful in this discussion, namely where f ijk and d ijk are the structure constants. Therefore, the commutator of two matrices ̺ 1 and ̺ 2 can be compactly expressed as a vector product of their Bloch vectors [̺ 1 , Notice that f ijk is antisymmetric under the exchange of any two indices.
We can now write the EOMs for the Bloch vectors of lepton number in an identical form as for two flavorṡ The antisymmetry of f ijk implies that this is once more a "precession equation" in the sense that D v moves in a direction orthogonal to D v so that its length is conserved. This formal analogy implies that all of the properties we derived in the previous sections are still valid in the three-flavor case. For example, we can perform the exact same mapping that leads us to the equations of slow flavor conversions. Furthermore, it is still true that we can associate to any set of Bloch vectors obeying the EOMs (56) a new set of Lax vectors obeying the same equations of motion. This means that a continuous set of Bloch vectors can still exhibit collective motions corresponding to a small number of degrees of freedom, e.g. a pendulum motion.

B. How Many Degrees of Freedom?
To investigate the integrability of the system, we need to identify both the number of dynamical degrees of freedom and the number of independent invariants. In the two-flavor case, the number of Lagrangian degrees of freedom of N Bloch vectors D v is exactly N , since the length of each D v is conserved, and two angles are needed to describe the motion, one of which is the canonical conjugate of the other. In addition, the length of the Lax vector L v is also conserved, providing N constants of the motion that are equivalent to the N Gaudin invariants (Sec. III C). Therefore, we had N linearly independent invariants and the system was therefore integrable in the mechanical sense.
What is the number of degrees of freedom in the threeflavor system? To answer this question, we first need to identify the analogue of the length conservation for the three-dimensional Bloch vectors. The density matrices obey a commutator equation of the form It is easily shown that the trace of ̺ v is conserved. As a consequence, the traces of the powers ρ n v are conserved as well for any n. We usually consider ̺ v to denote the traceless part of the density matrix which can be expressed in terms of the Bloch vector D v and in this case for n = 1, Tr ρ v = 0 identically. The non-trivial conservation laws come from n = 2 and n = 3, giving and where we define the triple product as A, B, C = d ijk A i B j C k . For higher values of n, we find no new conservation laws. These two invariants of motion are the analogue of the fixed-length property of the two-flavor Bloch vectors. Therefore, for N Bloch vectors with 8 components each, only 6 angles are needed to parameterize each of them. In analogy with the two-flavor case, we may say that the system admits 3N Lagrangian variables and 3N canonically conjugate variables. Another way of looking at the number of variables comes from observing that the diagonal elements of the original density matrix of one mode of the neutrino radiation field are the usual occupation numbers, whereas the complex off-diagonal elements encode flavor coherence. In the two-flavor case, there are two diagonal elements and one independent off-diagonal one (one complex number or two real parameters). Altogether there is the trace, the length of the Bloch vector, and two angle variables. In the three-flavor case, there are three real diagonal elements and three independent complex off-diagonal ones (six real parameters). Altogether we have the trace, the length of the Bloch vector, and its triple product, and in addition six angle variables encoding flavor coherence information.
The physical interpretation of the conserved-length parameters follows if we imagine that initially all neutrinos begin in flavor eigenstates so that ̺ = diag(f νe , f νµ , f ντ ) before making it traceless. We have already observed that the trace conservation implies that flavor evolution without collisions conserves ℓ f n ν ℓ for any n. In the two-flavor case, the two conserved quantities are particle number (f νe + f νµ )/2 and the coefficient of σ 3 in the Pauli-matrix expansion, which is (f νe − f νµ )/2. Its absolute value is the length of the Bloch vector. In general, (f νe −f νµ )/2 is the z component of the Bloch vector which is not conserved as flavor oscillations proceed. The initial quantity (f νe − f νµ )/2 is the conserved length of the Bloch vector if neutrinos start in flavor eigenstates.
In the three-flavor case, the relevant parameters are the coefficients in the Gell-Mann matrix expansion, i.e., the eight components of the Bloch vector P, and in addition the coefficient P 0 = 1 3 Tr̺ of the unit matrix, which is not part of the Bloch vector. To represent the initial ̺ that has only diagonal components, we only need P 0 and P 3,8 , the coefficients of the diagonal matrices λ 3,8 , and find The squared length of the Bloch vector P·P = P 2 is conserved under flavor oscillations. The third conserved quantity is the triple product P, P, P = √ 3P 2

C. Gaudin Invariants and Integrability
To identify the Gaudin invariants, it is more straightforward to express the Lax vectors in matrix form. The original EOMs are The Lax vectors are in vector and matrix form obeying the original EOMs From this follows in particular that the entire family of operators D n v L m v obey the EOMs and therefore their traces are conserved. This implies a new set of conserved quantities These are the generalizations of the two-flavor Gaudin invariants to three flavors. Notice that these invariants are all linearly independent of one another. Indeed, J v can obviously not be expressed in terms of I v , because it involves the tensor d ijk which does not appear in the definition of I v .
The third family of invariants, K v , can be written using the strategy that we adopted to derive Eq. (32) as but it cannot be further simplified, and therefore is not expressible in terms of J v . Just as in the two-flavor case, one can introduce additional sets of invariants which are linearly dependent on I v , J v , and K v , namely L 2 v and L v , L v , L v . We have already proved in the two-flavor case (Sec. III C) that L 2 v can be expressed in terms of I v alone. One can analogously prove that Therefore, we have 3N Lagrangian degrees of freedom, and 3N associated conservation laws for I v , J v , and K v . This guarantees integrability in the mechanical sense also for the three-flavor motion.

VI. CONCLUSIONS
The homogeneous nonlinear equations of motion for collective neutrino flavor oscillations have two oftenstudied limits, the single-angle case of slow oscillations and fast oscillations, corresponding to the multi-angle case in the limit of vanishing neutrino masses. In the classical limit (mean-field approach), both cases correspond to ensembles of classical interacting spins, but with different interaction structure.
Whereas the "slow Hamiltonian" corresponds to wellstudied cases in other fields, in particular the BCS Hamiltonian in the theory of superconductivity, the "fast Hamiltonian" appears to be less prominent. On the other hand, both cases have many similarities. It has been previously observed that both cases have soliton (pendulumlike) solutions, the nonlinear manifestation of the runaway solutions of a linear mode analysis.
We have systematically juxtaposed the two cases and developed the similarities and differences between them. Our starting point is the observation that the EOMs of the fast system can always be mapped on an abstract slow system. For the special case when the fast system has a pendulum solution, our transformation has the intuitive interpretation that, from the perspective of the moving pendulum, it is the "direction of gravity" the moves like a "reciprocal" pendulum with EOMs formally equivalent to a slow system. While transforming the EOMs to co-rotating frames is a standard procedure in this field, probably this is the first time one has used a more general transformation.
On the other hand, a slow system can be mapped on an equivalent fast system only if its system of spins (or Bloch vectors) obeys an additional constraint. This fundamental difference between the two cases is traced to the presence of an external vector in the slow system, the mass direction B in flavor space, whereas the fast system has no preferred direction except for the conserved total angular momentum, which is determined by initial conditions, not the Hamiltonian itself.
We have adapted the powerful tool of Lax vectors, introduced in this field in studies of the BCS Hamiltonian, to the fast system. In this way, it is straightforward to identify the invariants of the systems and especially the Gaudin invariants, but also offers another perspective on the normal mode analysis, allowing one to derive the dispersion relation in a single line.
In the continuum (thermodynamic) limit, this dispersion relation covers only the collective modes, whereas the non-collective ones require a different interpretation. Using the equivalence of our linearized system to the Vlasov equation, we highlight the equivalence of the noncollective modes to the well-known Case-Van Kampen modes. In this way, decoherence of a generic initial condition is a natural feature of our system. The Gaudin invariants immediately prove the integrability of the system. We explicitly show that the integrability also translates to the quantum case.
These techniques and many of our conclusions carry over to the three-flavor case, where the Bloch vectors and Lax vectors live in an eight-dimensional space. While part of this extension is straightforward, it involves a number of nontrivial steps and new invariants, owing to the SU(3) structure of this problem.
Relative to the question of neutrino flavor evolution in real astrophysical environments, our homogeneous system is probably too constrained to serve as a realistic proxy. On the other hand, it provides surprising insights on the foundational properties of the nonlinear EOMs and their possible solutions.

Appendix A: Comoving frames
Precession equations of the type used throughout this paper can become more transparent in a moving frame in flavor space. Frames rotating around the z-direction [10] have become standard. However, one has to be careful to change everything consistently. In the fast EOM (7) we may simply leave out D 0 by going to a frame rotating around the fixed D 0 with unit frequency. On the other hand, going to a frame co-rotating with frequency w c around B in the slow EOM (6) changes w → w − w c .
To derive the formal transformation between two frames we consider the precession equationṡ where P w is a Bloch vector that depends on some attribute w and precesses around some other vector J w that may itself depend on time. Moreover, we consider another precession equation of the forṁ and we want to express the former in a frame that comoves with R.
The formal transformation is more easily derived using instead the commutator form of the EOMs. We recall that any 2 × 2 Hermitean matrix P can be expressed as P = 1 2 (Tr P + P · σ) in terms of the Pauli matrices and a Bloch vector P. The commutation relations of the Pauli matrices imply that an EOM of the form is represented by Tr P = 0 andṖ = H × P. The EOM (A3) is formally solved by where T is the time-ordering operator. This transformation provides the compound effect of infinitesimal precessions around H t that is itself moving. The content of Eq. (A4) can also be expressed as where it was assumed that H is a Hermitean matrix. One can verify that indeed We now denote matrices in the new frame with a tilde, i.e.,P w = U † P w U. In particular,R = U † RU = U † UR 0 U † U = R 0 . So indeed R = R 0 at all times, which was the original goal of the transformation. The equivalent of Eq. (A1) is i∂ t P w = [J w , P w ]. In the new frame it is explicitly i∂ tPw = i∂ t (U † P w U) = (i∂ t U † )P w U + U † P w (i∂ t U) + U † (i∂ t P w )U = −U † HP w U + U † P w HU + U † [J w , P w ]U = −[H,P w ] + [J w ,P w ] and finally In this form, our result also applies to three flavors. Returning to Bloch vectors, the result is equivalent to This is the same prescription that was argued in the main text: To take out the motion engendered by H we must subtractH in the transformed frame from the vectorJ w that generates the original motion.
Appendix B: Linear stability analysis in the thermodynamic limit The dispersion relation Eq (40) is unsatisfactory from a mathematical standpoint, because it involves an integral which is singular for real frequencies. A prescription on how to deal with this singularity should be given. This issue has historically first been solved in the context of collisionless plasma oscillations. Indeed, Eq. (38) is mathematically identical to the linearized Vlasov equation describing the evolution of an initial spatially monochromatic small perturbation δf (r, v) = δf k (v)e ik·r in the distribution function of electrons in a plasma (here v is the electron speed) 7 where f 0 (v) is the unperturbed distribution, e and m are the electron charge and mass; see, e.g., Ref. [35]. Let k be aligned with the z axis. We now define and similarly Integrating Eq. (B1) over v x and v y , we obtain This equation is indeed identical to Eq. (38), with the identification kv z = w. The singularities in the dispersion relation appear in the same way in this analogous problem; below, we summarize the main results from the literature concerning this issue. The eigenmodes of the linearized problem can be separated into collective and non-collective modes, as is done in Ref. [36]. The collective modes p xy w,i appear as discrete complex solutions Ω i of the dispersion relation (40), and, as shown in the text, they always appear in pairs of complex conjugate solutions. The number of such collective modes depend on the spectrum g w ; for example, for noncrossed spectra there are no collective modes. The non-collective modes are not predicted by the dispersion relation (40). Rather, they correspond to a continuum of modes with real frequencies Ω v = vB z − P z 0 , and do not show up in our dispersion relation because their eigenfunctions are singular where The pole in the eigenmode (B5) is meant to be integrated in principal value (PV). By explicit substitution, one can verify that these eigenmodes indeed satisfy Eq. (38). These modes were first identified as the so-called Case-Van Kampen modes [37] in the case of collisionless plasma oscillations, where it was shown that the collective and non-collective modes together form a complete set of eigenmodes for the linearized Vlasov equation. By the same token, in our problem of neutrino oscillation, any initial perturbation P xy w (t = 0) can be decomposed in collective and noncollective eigenmodes which evolve independently, so that the full evolution can be written (see, e.g., Ref. [38]) where the coefficients a i and b(v) are determined by the initial conditions. However, the integrated lepton number P xy (t) = dwP xy w (t) has an entirely different dynamics in the continuum limit. The reason is that the superposition of the Case-Van Kampen modes, evolving with a continuum of frequencies in time, translates into a damped behavior. This behavior is the exact analog of Landau damping in collective plasma oscillations; even though the individual polarization vectors have non-decaying terms in their perturbation, because of their fast precession with incommensurate frequencies they lead to a perturbation in the integrated lepton number that decays in time.
The dynamics of P xy (t) could be obtained by explicitly integrating Eq. (B7). However, it is easier to circumvent the expansion in eigenstates and directly look at the evolution in time of the initial condition using the Laplace method, following the strategy in Ref. [35]. By this method, one finds that P xy (t) evolves in time as where A α are determined by the initial condition andΩ α are the solutions of the modified dispersion relation where the contour of integration C is chosen as the real axis if Im(Ω α ) > 0, therefore identical to Eq. (40), whereas for Im(Ω α ) < 0 the contour of integration is deformed so as to pass below the pole w = (Ω α + P z 0 )/B z . The difference in the contour definition for growing (ImΩ α > 0) and shrinking modes (ImΩ α < 0) originates from causality requirements. Because of this difference, the frequenciesΩ α do not appear in complex conjugate pairs. The modified dispersion relation admits new shrinking modes without the corresponding growing ones. Physically, these shrinking modes come from Landau damping of the initial perturbation, due to the phase mixing of the many polarization vectors with incommensurate frequencies. The symmetry between growing and shrinking modes is therefore broken by the second principle of thermodynamics, which implies the decoherence of the motion of individual polarization vectors into a damping of the integrated lepton number.
A first, fundamental difference to the slow flavor system is the following. The Hamiltonian (C1) is invariant under a simultaneous rotation of all polarization vectors. This is not the case for the slow system, where there is a privileged direction fixed by the B vector. The rotational invariance of the fast flavor system implies the conservation of the generator of collective rotation which indeed commutes with the Hamiltonian. Therefore, we reach again the conclusion stated multiple times in the main text that the fast flavor system is analogous to a slow flavor system endowed with an exactly conserved vector. For the fast flavor system, such a vector is D 0 , which is conserved by rotational invariance.
Let us now determine the eigenstates. The Hamiltonian (C1) admits two trivial eigenstates, in which all the spins are aligned along, say, the z direction. We will use the notation |S, S z to denote the eigenstates of the total spin with eigenvalues S and z component S z . Then the two eigenstates identified above are | N 2 , ± N 2 . However, because of rotational invariance, there is nothing special about the z direction, and these two states rotated in an arbitrary way are still eigenstates of the Hamiltonian with the same eigenvalues. Therefore, the entire multiplet of 2N +1 states | N 2 , S z , with S z ranging from − N 2 to N 2 , are degenerate eigenstates of the Hamiltonian. This stands in sharp contrast to a generic slow flavor system, for which only the two eigenstates | N 2 , ± N 2 with spin totally aligned to the z axis defined by B are eigenstates.
We now need to identify the remaining eigenstates. We will classify them in terms of their z component of spin S z , and we will start looking for eigenstates containing all spins aligned, for example negatively, along the z axis, except for a single spin flipped (therefore, S z = − N 2 + 1). If there were no interaction, all N states of the form | ↑↓ ... ↓ , | ↓↑ ... ↓ ,..., | ↓↓ ... ↑ would be degenerate eigenstates of the Hamiltonian. The presence of the interaction splits the degeneracy and selects N combinations of these states, which we now identify; these states are essentially the analog of the magnon excitations in the one-dimensional Ising model, with the difference that our Hamiltonian (C1) breaks translational invariance, so that the form of the eigenstates is not as simple as in the Ising model.
To identify the form of the eigenstates with a single spin flip, we use the Bethe ansatz and assume that they can be written as where Q u is an operator we are now going to identify. In order to do so, let us recall that the Lax vectors introduced in the main text follow a pure precession motioṅ Defining the vector h u = D 0 − uD 1 implies that where L ± u = L x u ± iL y u and so forth. We now make the ansatz that the operator Q u introduced above coincides with L + u , and determine how the Hamiltonian acts on the state |ψ u The last term is H| N 2 , − N 2 = E 0 | N 2 , − N 2 , where E 0 is the ground state energy. In the first term, we notice that L + u h u z | N 2 , − N 2 = e u |ψ u , where e u is the eigenvalue of the operator h u z acting on the ground state | N 2 , − N 2 , namely Therefore, we obtain Therefore, if u is chosen as one of the roots of L z u | N 2 , − N 2 , namely one of the roots of the Bethe equation the last term vanishes and the state |ψ u is indeed an eigenstate of the Hamiltonian. The Bethe equation admits therefore N −1 eigenstates with a single spin flip; an additional such eigenstate is contained among the multiplet of 2N + 1 eigenstates of the total spin D z 0 , namely the state | N 2 , − N 2 + 1 . This procedure can be generalized to states with more than one spin flip. For example, one can prove that the commutator Therefore, the state L + u L + v | N 2 , − N 2 is an eigenstate of the Hamiltonian provided that the last two terms in parenthesis vanish identically, namely if the Bethe equations are satisfied with the compact notation L z u | N 2 , − N 2 = J z u | N 2 , − N 2 . The eigenvalue of the Hamiltonian on this new eigenstate is E 0 + e u + e v + 1 − uv.
Therefore, the repeated application of the operator L + u1 ...L + un on the ground state generates a novel eigenstate with eigenvalue E 0 + i e ui + 1 2 i,j =i (1 − u i u j ), provided that the Bethe equations are satisfied which completes our proof.