Collisions and collective flavor conversion: Integrating out the fast dynamics

In dense astrophysical environments, notably core-collapse supernovae and neutron star mergers, neutrino-neutrino forward scattering can spawn flavor conversion on very short scales. Scattering with the background medium can impact collective flavor conversion in various ways, either damping oscillations or possibly setting off novel collisional flavor instabilities (CFIs). A key feature in this process is the slowness of collisions compared to the much faster dynamics of neutrino-neutrino refraction. Assuming spatial homogeneity, we leverage this hierarchy of scales to simplify the description accounting only for the slow dynamics driven by collisions. We illustrate our new approach both in the case of CFIs and in the case of fast instabilities damped by collisions. In both cases, our strategy provides new equations, the slow-dynamics equations, that simplify the description of flavor conversion and allow us to qualitatively understand the final state of the system after the instability, either collisional or fast, has saturated.


I. INTRODUCTION
Collective flavor oscillations in neutrino-dense environments are driven by neutrino-neutrino refraction that also possesses off-diagonal flavor components [1].The magnitude of the induced refractive energy shift is set by the neutrino number density n ν , and can roughly be measured by the scale µ = √ 2G F n ν , with G F the Fermi constant; everywhere in this work, we use natural units such that c = ℏ = k B = 1.For typical conditions of a supernova (SN) and the remnant of a compact object binary, the neutrino-neutrino interaction strength can be as large as µ = 10 5 km −1 in the vicinity of the neutrino decoupling regions and can drive neutrino fast flavor conversions (FFC) that operate on timescales of a few nanoseconds.The timescale µ −1 is much smaller than the one of neutrino slow oscillations √ µω −1 , where ω = ∆m 2 /2E ν ; for a typical neutrino energy of 15 MeV and the largest mass-squared difference one finds ω ≃ 0.4 km −1 .Since the potentials span over a wide range, collective flavor oscillation can also span over a wide range of time and length scales, for recent reviews see Refs.[2][3][4][5][6].
Recently, there has been a surge of interest in the interplay between the coherent forward scattering (∝ G F ) and the incoherent scattering of neutrinos off matter (∝ G 2 F ). Roughly, the ratio of coherent to incoherent rates is given by Γ/µ ≃ G F E 2 ν ∼ 10 −9 for a 15 MeV average neutrino energy; the regime where Γ ≪ µ is referred to as the weak-damping regime.Otherwise, if Γ ∼ µ oscillations are generally damped, similarly to the quantum Zeno effect [7], but such conditions never seem to occur in realistic astrophysical or cosmological environments.
Conventional wisdom suggests that incoherent collisions should damp the off-diagonal coherence required to exhibit collective motion, and therefore destroy collective oscillations.However, as realized recently, this is not always the case if the collisional rates depend on energy or differ for neutrinos and antineutrinos [8][9][10][11][12][13][14], which might trigger novel collisional flavor instabilities (CFIs).Even if CFIs have a subdominant role, collisions might still play a role in damping the fast oscillations which likely occur in the inner regions of SNe.
Investigating the impact of collisions on the final state of flavor conversions requires an investigation of the nonlinear evolution of the latter.To tackle this issue, hierarchies among different scales are one key to simplifying the problem.As mentioned before, one such hierarchy is Γ ≪ µ.Therefore, whether collisions drive instabilities or simply damp the collective fast oscillations driven by forward scattering, their dynamics is slow, and the separation of these scales can be used to simplify the problem.
A similar separation is key to understanding the traditional MSW effect [15][16][17] where fast in-matter oscillations separate from the slow change of density along the beam.Likewise, in a dense medium, a large matter effect is often lumped into a small effective in-medium mixing angle and otherwise ignored because the fast oscillations can be averaged out.On the other hand, for oscillations driven by neutrino-neutrino refraction, such separationof-scale techniques have not been explored.This is a more challenging scenario because even the rapid oscillations are nontrivial, being driven by nonlinear dynamics.
Here we address precisely this problem.We show that due to the hierarchy Γ ≪ µ, the dynamics of oscillations impacted by collisions can be separated into a fast dynamics, on timescales µ −1 , and a slow dynamics 1 , on timescales Γ −1 .In practice, the separation is performed by averaging over the fast dynamics, and identifying a set of physically relevant quantities that only change slowly.This allows us to obtain a new set of equations, henceforth referred to as slow-dynamics equations (SDEs), which do not contain terms of order µ.Therefore, solving the SDEs requires to follow the dynamics only over timescales of order Γ −1 , without keeping track of the fast motion which is not relevant to the final flavor outcome of the evolution.Admittedly, these new equations are still nonlinear and coupled; yet, they can be solved much faster than the original equations of motion (EOMs), while offering significant advantages in terms of intuitive understanding.
We apply this approach to two complementary problems.First, we investigate a homogeneous neutrino gas that exhibits a CFI without a fast instability.Representing the neutrino flavor coherence in terms of polarization vectors, the fast dynamics is a simple precession of the polarization vectors around the common flavor axis, so the amplitudes and phases of this precession are left invariant by the fast dynamics.The slow change of these quantities induced by collisions is described by a set of SDEs.We verify that the SDEs lead to the same evolution as the standard EOMs, and we use them to conclude that after the instability has saturated, flavor equipartition is reached for low energies while no conversion is attained for high energies.
Second, we consider a homogeneous neutrino gas that possesses a fast instability, subject to energy-independent collisions.In this case, the fast dynamics, without collisions, is pendular in nature [18,19].Therefore, the slowly-varying quantities can be chosen as the pendulum parameters.In this case, the SDEs allow one to understand immediately several features that had previously only been recovered empirically, such as the observation that the final outcome of the oscillation only depends on the depth of the pendulum and not on any other parameter or the collision rate [18].
The structure of the paper follows three main subjects.In Sect.II, we consider a toy system consisting of an active and a sterile neutrino species subject to oscillations and collisions.This system has been thoroughly studied in the literature and we use it as a simple introductory example for fast dynamics, driven by oscillations, that can be cleanly separated from a slow dynamics, driven by collisions.In Sect.III, we discuss the SDEs for a homogeneous system that exhibits CFIs.In Sect.IV, we discuss the SDEs for a homogeneous system subject to a fast instability.Finally, in Sec.V, we give a general discussion of our method and its limitations.

II. ACTIVE-STERILE NEUTRINO CONVERSION
The simultaneous appearance of fast and slow dynamics is more general than phenomena in collective neutrino oscillations.A first example for such a separation of time scales is the production of sterile neutrinos ν s from a population of active ones ν a in the early Universe [20][21][22][23][24][25][26][27] or in SN cores [28][29][30][31][32][33][34][35][36].Even in the standard cosmological model, averaging over the rapid oscillations has proven a good strategy to simplify the numerical solution of the flavor evolution of the neutrino plasma [37,38].For illus-tration, we use the simplest toy model, a homogeneous gas of ν a and ν s , no cosmic expansion, and only ν a interacting with the ambient plasma.The ν s production rate is considerably simplified by the oscillation rate being much faster than the ν a collision rate.This simplification is analogous to the slow-dynamics approximation that we will use in the later collisional and fast instabilities examples.
The neutrino radiation field of this simple two-flavor system is described, in the mean-field approximation, by a 2 × 2 density matrix ρ p (t) for every momentum mode p.Assuming also isotropy, it depends only on p = |p| and we write it in the form where the time dependence of ρ p is not explicitly shown.
The entries are the usual occupation numbers, defined as f αβ (p) = ⟨a † β,p a α,p ⟩ with α, β = a or s, i.e., the expectation values of number operators, where the mixed off-diagonal elements encode flavor coherence.For the diagonal elements, we do not show the repeated flavor index and we note that the matrix is Hermitian with f sa (p, t) = f * as (p, t).In the absence of collision, ρ p evolves by flavor oscillations caused by masses and flavor mixing and modified by refractive effects of the ambient medium.Without spelling out the details of the latter, we simply assume for every p a Hamiltonian 2 × 2 matrix H 0 p for every mode that drives flavor oscillations according to (We use sans-serif letters such as H to denote matrices in flavor space.)We assume that H 0 p does not depend on the neutrino medium itself so that this EOM is indeed linear in ρ p .In this case, it is explicitly solved by which can be spelled out analytically because H 0 p itself does not depend on time.
The eigenstates 1 and 2 of H 0 p are the propagation eigenstates, which for simplicity we call mass eigenstates, with the energies , where the effective masses depend on p, the relativistic approximation was taken, and we have removed the common p from the definition of H 0 p .The diagonalization of H 0 p is enacted by the in-medium mixing angle θ through 1 2p where it is understood that θ and m 1,2 depend on p.
In addition, the active states interact with the medium, breaking ν a -ν s flavor coherence.We model this effect by a scattering rate Γ p→p ′ between ν a states and denote with Γ p = p ′ Γ p→p ′ the total loss rate.With the matrix the modified EOMs are [39] where we neglect Pauli blocking factors.The separation of fast and slow timescales naturally happens if the scattering rates are much slower than the oscillations induced by H 0 p ; more precisely, we require that Γ p→p ′ is much smaller than the difference between the eigenvalues of H 0 p , i.e., ∆m 2 /2p ≫ Γ p .In this case, neutrinos oscillate many times between collisions and we can significantly simplify the equations.Based on the nocollision solution of Eq. ( 3) we first pass to the interaction representation by where the density matrix ϱ p in the interaction representation would be constant in the absence of collisions.Including collisions, the EOM becomes ρp = We can now perform an average over many periods of the rapidly oscillating exponentials.In the mass basis, this means that the off-diagonal entries of any matrix average to zero.For the density matrices, this means explicitly where f 1 (p) and f 2 (p) are the distribution function of the two mass eigenstates.After inserting this parameterization in the EOM (8), we obtain finally These equations refer to the mass-basis distribution functions, which coincide with the flavor-basis distributions if the mixing angle is small.In this case, if the ν a are in equilibrium, one finds that the ν s production rate equals the ν a scattering rate times a factor sin 2 θ cos 2 θ = sin 2 (2θ)/4.True equilibrium, however, can only obtain in the mass basis.
The main feature of this formulation is that the righthand side of Eqs.(10) only contains the scale set by the collision rate.The terms depending on the rapid oscillations have disappeared after we have averaged them out over many periods of oscillation.In the remainder of this paper, we will apply the same strategy to the much more complicated case in which the rapid oscillations are induced by neutrino-neutrino refraction.

III. COLLISIONAL INSTABILITY
A. Setup of the problem As a first nontrivial example of the slow/fast dynamics separation, we consider a homogeneous neutrino gas without a fast instability which exhibits CFIs with energy-dependent collisions.We consider the limit of vanishing neutrino masses, i.e., ω = ∆m 2 /2p → 0, which are not required for driving the CFI dynamics.We limit ourselves to an isotropic, two-flavor setting with neutrinos subject to mutual forward scattering and numberchanging processes with an external medium.In the limit of vanishing neutrino masses and for a homogeneous and isotropic system, the refractive matter term is eliminated by a uniform rotation around the flavor axis.Moreover, the ∆m 2 → 0 limit also implies the absence of vacuum mixing.Of course, we still imagine that the mass term provides an initial seed for the CFI.
As in Sect.II, the system is described by density matrices ρ E (t) for the flavors e and x, where the latter stands for µ or τ or an admixture of both.We here use E ≃ |p| to label the modes and in addition, we use a negative sign for E to denote the density matrices for antineutrinos.On the other hand, we do not use the flavor isospin convention, i.e., the upper left entry of ρ E is the ν e or ν e occupation number.In the language of polarization vectors to be used later, this convention means that a given P E pointing in the positive z-direction means ν e or ν e , depending on the sign of E, and pointing down the other flavor.The evolution of the system in the absence of collision is now ruled by the equation where s E ′ = sign E ′ needs to be explicitly included.
In order to avoid cluttering notation, and to connect with the definition of a typical energy scale, we define µ = √ 2G F n 0 ν , where n 0 ν = 3T 3 ζ 3 /4π2 is the equilibrium number density of a neutrino population at the temperature T of the thermal bath.We then define a sum over energies with a modified measure E = E 2 dE/(2π 2 n 0 ν ).Therefore, the density matrix for the full ensemble is With these definitions, in the absence of collisions, the EOM for every mode has the usual form ρE = iµ[ρ E , ρ]; µ is a typical energy scale associated with the selfinteraction, while the sum over energies, as well as ρ E , is now dimensionless.Notice that, as defined here, µ does not fully coincide with the usual definition in terms of the total neutrino number density n ν (see Sec. I); in the context of CFIs such a definition cannot be adopted, since the number of neutrinos is not conserved.
To model the collisional effect, we envision a SN medium where ν e and ν e are efficiently absorbed or emitted by beta processes, whereas the x flavor does not interact at all.In this sense, our setup is similar to the activesterile system of Sect.II.The EOM with beta processes included can be written 2 in the form [40] where in the flavor basis the production and absorption matrices P E and A E have only upper-left entries like B in Eq. ( 5), relevant for ν e and ν e .
Here, Q E is the spontaneous emission rate of ν e by the medium, Γ E the reduced absorption rate, and analogous for ν e .Ignoring temporarily the matrix structure and considering only the e-flavor, the collision equa- rate, although it is actually enhanced, the terminology deriving from photon radiative transport, where the Bose stimulation factor has the opposite sign from the Pauli blocking factor.Moreover, if the medium is in thermal equilibrium, detailed balance implies for ν e (positive E) Q E = Γ 0 E e −(E−µν e )/T and thus where µ νe is the ν e chemical potential implied by the properties of the medium.Using Γ E as our primary parameter for the damping strength, the production rate is where ± refers to ν e and ν e .In our convention, the latter are denoted by negative E and have the opposite chemical potential.Also, for negative E, Γ E is the reduced absorption rate for antineutrinos.The equations are most conveniently expressed in terms of the total number P 0 E = Tr ρ E = f e (E) + f x (E) of neutrinos in mode E and the polarization vectors P E , such that In terms of these variables, the EOMs become Here e z identifies the direction of the flavor axis, and It will prove convenient to express them in terms of a modified neutrino density N E = P 0 E − 2Q E /Γ E , a quantity that need not be positive.In thermal and chemical equilibrium, both e and x flavors are occupied with the Fermi-Dirac distribution f e,x (E) = 1/[e (E−µν e )/T +1] and P 0 E = f e (E) + f x (E) with Eq. ( 14) reveals that N E is the deviation from equilibrium.Finally, is the set of EOMs that we plan to solve.

B. Slow-dynamics equations (SDEs)
We now proceed to integrate out the fast dynamics in the spirit of our earlier discussion.When Γ E = 0, the exact solution consists of a fast precession of each P E around P. This motion is simple to understand and has many conserved quantities.In the absence of collisions, the total P = E s E P E is conserved, reflecting total lepton-number conservation.Furthermore, each P E moves uniformly around P on a cone with a fixed angle.To write the general motion explicitly, we introduce a right-handed coordinate system {e 1 , e 2 , e 3 }, where e 3 is along P, whereas e 1 = e z × P/|e z × P| is orthogonal to both P and the flavor direction e z , and finally e 2 = e 3 × e 1 .In other words, e 1 and e 2 span the plane transverse to P. The projection of P E on the direction of P is expressed through α E = P E • e 3 /P , whereas the length in the transverse plane is b E = |P E × e 3 |.In this way, the explicit solution is where the ϕ E are fixed individual phases and Φ t = µP t the common oscillation phase.When we later assume that P slowly changes, it will be Φ t = µ t 0 dt ′ P (t ′ ).The solution is entirely specified by the conserved quantities P, the projection α E of each P E on the precession axis, the amplitude of precession b E , and the phase ϕ E .Because P = E s E P E (t), these quantities are constrained by and Collisions cause these quantities to slowly evolve over timescales given by the inverse damping rate.By the assumption of weak damping, each P E performs many precession cycles over this timescale and therefore the precession itself can be averaged out.In analogy to Sect.II, we replace the form of the solution Eq. ( 18) in the full EOMs with damping, where we now consider the parameters P, α E , b E , e 1 , e 2 , and ϕ E as dynamical quantities.Notice that the directions e 1 , e 2 and P themselves depend on time and therefore must be differentiated as well.To average over the rapidly varying phase Φ t , we use ⟨sin , and the mixed ⟨cos . Performing these averages directly provides us with the SDEs: These equations are admittedly not much simpler than the original EOMs; they are still nonlinear and couple neutrinos of all energies.However, they offer some practical and conceptual advantages.There is no term of order µ and thus a numerical solution requires only a relatively coarser grid, with a step of order Γ −1 that mirrors the scale over which the quantities change.More importantly, the new equations provide some intuitive understanding of the role of damping, by cleaning the dynamics of the complicated precession terms.As a first step, we can eliminate the phase variables ϕ E .From Eq. (19c) we see that the phases follow a Kuramotolike dynamics [41][42][43], with the damping term tending to synchronize the phases of precession and make all the polarization vectors coplanar.If two vectors P E and P E ′ have identical phases ϕ E = ϕ E ′ , one can see from Eq. ( 18) that at a given time they lie in the plane spanned by P and cos(Φ t + ϕ E )e 1 + sin(Φ t + ϕ E )e 2 .(Notice that for a two-beam case the two vectors are necessarily co-planar, so one can exactly set ϕ = 0.) One can always choose the initial condition such that all P E are exactly coplanar, and the subsequent dynamics likely does not depend on this choice, since anyway they would soon become coplanar because of phase locking.In the following, we simply choose ϕ E = 0 for every E.
We can now understand the existence of two separate branches of CFIs from the SDEs directly.Starting from Eq. (19d), let us first consider the component transverse to the flavor direction, Thus, the evolution of the transverse component is driven by the effective damping rate E s E Γ E α E .If initially this expression is negative, it becomes an effective growth rate, leading to an instability and coincides with one of the two branches of instability that have been identified in the literature [8][9][10].
Even when this branch corresponds to damping, an instability can still arise from the growth of the individual precession amplitudes b E .To see how this can happen, we take the simplest case of a two-bin model with b = b and α = α−1, where we denote by α, b, and Γ the quantities for monochromatic neutrinos, and by α, b, and Γ the ones for antineutrinos.Equation (19b) then simplifies to ḃ = α(Γ − Γ)b − Γb that can be written as If the expression in brackets is positive, b grows, corresponding to the second branch of instability identified in the previous literature.In the two-bin model, the first branch corresponds to (αΓ − Γα) being positive (growth) or negative (damping).Thus, for this simple two-bin system, we recover the linear stability results for the isotropic, non-resonant collisional instability [11,12].The identification of the instability allows one to provide a clearer criterion for the applicability of the SDE approach.So far, we have schematically required a hierarchy of the form µ ≫ Γ, but since the interaction rate Γ E is actually energy-dependent, this requirement should be made more specific.The typical timescale over which the CFI evolves is indicated by the collective growth rate, which we denote by γ, determined by linear stability analysis.On the other hand, the precession of the polarization vectors happens on typical timescales of the order of µP z (t = 0).Therefore, a more concrete criterion for the applicability of the SDEs is that µP z (t = 0) ≫ γ.
C. Understanding the nonlinear outcome

Simplified equations of motion
The second instability, in the two-bin case represented by Eq. ( 21), corresponds to the presence of more ν e than ν e , and ν e more strongly damped, which is the situation in a SN core with regard to beta processes.Therefore, we focus on this case and gather some physical intuition about its impact on the flavor evolution.This case corresponds to damping for the first type of mode, i.e., damping in Eq. (19d) so that P never develops a large transverse component if we begin with all P E aligned with the z-direction except for a small seed.Therefore, we can simplify our equations by taking P = P z .
We can further simplify the EOMs by writing them in terms of dimensionless parameters through N E = ν E P and b E = β E P .The interpretation of the new parameters is that for every E, α E is the dimensionless z component of P E , β E the dimensionless amplitude of precession, and ν E represents the dimensionless number of neutrinos in that mode.We further define The EOMs of the dimensionless variables are then whereas regulates the dynamics of P .For the absorption rates we assume a quadratic energy dependence Γ 0 E = Γ(E/3T ) 2 for ν e (positive E) and Γ(E/3T ) 2 for ν e (negative E), where Γ and Γ are positive parameters.As explained in the text above Eq.( 14), the reduced damping parameters then have the form 1 + e −(E−µν e )/T for E > 0, (25a) 1 + e +(E−µν e )/T for E < 0, (25b) where µ νe is the ν e chemical potential defined by the medium properties.Here Γ and Γ are the damping rates for typical ν e or ν e energies.
Our simplified EOMs immediately reveal that the final state only depends on the ratio Γ/Γ because all damping parameters are linear in Γ E so that Γ can be pulled out and absorbed in the definition of the units of time.

Numerical example
To build some intuition, we now turn to a specific example, for which we solve numerically the newly derived SDEs 23.As initial distributions, we choose thermal ones for ν e and ν e corresponding to the distribution of beta equilibrium with the medium.Here we use ϵ = |E| as a positive energy variable.For ν x and ν x , we assume that their initial population essentially vanishes, so they can only be produced by flavor conversion.The numerical values chosen for the parameters are listed in Table I.We also show the initial number densities that follow from Eq. ( 26).The number density of neutrinos with vanishing chemical potential for T = 4 MeV is n 0 ν = 5.864 MeV 3 .With our definition of the sum over energies, the initial value of the total polarization vector is P z (0) = 1 2 (n νe − n νe )/n 0 ν = 0.468.As discussed earlier, the absolute values of Γ and Γ are irrelevant, and only their relative values to one another matter, as long as we are in the weak-damping regime Γ, Γ ≪ µ.We choose the ratio Γ/Γ = 0.1 and measure time in units of Γ −1 .We initialize the system with seeds of small randomly chosen off-diagonal components P x E and P y E (|P x E , P y E | ≲ 10 −3 P z E for all energies) .As a reference comparison, we solve the SDEs and compare them with the solution of the full EOMs for µ = 50 Γ.Linear stability analysis shows that the growth rate of the CFI is γ ≃ 0.39Γ, so that µP z (0) ≃ 122 γ, justifying the applicability of the SDE approach.
As an example of the simplification induced by the SDEs, we show in Fig. 1 the evolution of the squared transverse component (P x E ) 2 of a specific polarization vector, corresponding to an energy E = 9.5 MeV, where the ν e spectrum peaks, both solving the full EOMs and the SDEs.For the latter, from Eq. ( 18) we see that, since P remains aligned with the z axis, the mean squared component is simply This has a simple interpretation; since b E is the amplitude of the precession around the z axis, the squared amplitude along a specific direction is simply one-half of the transverse squared amplitude.The comparison in Fig. 1 shows that the solution of the SDEs tracks the average of the exact solution of the EOMs over the very rapid oscillations.
Figure 2 shows the evolution of this system, as predicted from the SDEs, as well as the final state reached after the CFI saturates.The final ν e and ν e distributions are identical to the initial ones, both corresponding to beta equilibrium with the medium, but during the evolution they deviate.By assumption, there is no initial population of ν x and ν x .The dashed lines show the distributions at the intermediate stage chosen at the instant when P z is maximal as indicated by the vertical dashed line in the right panels.
The CFI is manifested in the growth of the in-plane components of the polarization vectors, namely of the parameters β E introduced above.Correspondingly, the collective amplitude B = ϵ s E Γ E β E grows exponentially in time.We easily check that the instability for these conditions belongs to the second branch identified above, since the in-plane components P x and P y never grow, so the total polarization vector remains aligned with the z axis.Its magnitude initially grows exponentially.After the instability grows nonlinear, both P z and B return to their original value; in the case of B, this means that the polarization vectors return to align with the z axis.We also show the evolution of the function A defined in Eq. (22a); the complementary function N defined in Eq. (22c) remains always essentially identical to N ≃ −A.As the instability grows, the ν e and ν e population are temporarily depleted since they convert into ν x and ν x , as visible from the dashed lines in Fig. 2.
After the instability has saturated, ν e and ν e return to their initial spectrum.This is only to be expected since the equilibrium state of the e flavor is determined by thermal and chemical equilibrium with the medium and is enforced by the collisional term.On the other hand, ν x and ν x do not interact with the medium, and therefore their final amount is not so trivially understood.The simulation reveals that their final distributions are characterized by an equipartition of ν e and ν x for low energies.At high energies, the ν x remain unpopulated and are not efficiently produced by flavor conversion.

Low-energy behavior
To show that these main features identified in the numerical solution descend naturally from the SDEs, we begin with the low-energy behavior.Beginning with Eq. ( 23), we first perform the transformation and a similar transformation for β E and ν E .In terms of the new variables, the equations become We now notice that the behavior of the system is critically different according to whether the ratio Γ E /B is very large or very small.Since the function B depends on time and is initially very small, clearly for all polarization vectors we have initially B ≪ Γ E .However, as time progresses, B grows exponentially until a maximum value B max ; for the numerical example analyzed earlier, B max ≃ 0.35 Γ, as seen from the right panel of Fig. 2.
The coefficients Γ E have a quadratic energy dependence, while the Fermi-Dirac factor does not have a dramatic impact, so we can approximately identify a range of energies E ≲ 3T B max /Γ for neutrinos and E ≲ 3T B max /Γ for antineutrinos.For these energies, we may get a qualitative insight by considering the limit Γ E → 0, where Eqs.(28) result in ( αE − βE ) = −B(α E − βE ) which implies that αE − βE decreases exponentially in time as e − t 0 B(t ′ )dt ′ .Therefore, in this lowenergy range, α E and β E are maintained approximately equal.Since β E finally reaches the value of zero, after the transverse motion has exhausted, we conclude that in this low-energy range, we must finally have α E = 0, i.e., equipartition among the e and x flavor.Our numerical simulation completely verifies this prediction, as seen in Fig. 2. Notice that this allows us to conclude that the maximum energy at which equipartition is reached for neutrinos and antineutrinos, Ẽν and Ẽν respectively, are related by Ẽν / Ẽν ≃ Γ/Γ, a conclusion that is verified by our numerical simulation.

High-energy behavior
At high energies, such that Γ E ≫ B max , the qualitative behavior of the solution can be understood for the limit B → 0, where we find (α 2 E − ν2 E ) to be constant at all times with αE , νE ∝ e −iΓ E t .Moreover, β E never FIG. 2. Numerical solution of the SDEs for the model system described in the text.Left: Final (solid) and intermediate (dashed) spectrum of νe and νx (νe and νx are shown on the negative energy axis).For νe and νe, the final spectrum is identical to the initial one, corresponding to beta equilibrium with the medium.The intermediate spectrum corresponds to the instant in time when P z is maximum, identified by the dashed line in the right panels.The vertical dotted line identifies the energy of the polarization vector whose dynamics is illustrated in Fig. 1.Right: Time evolution of the integrated quantities P z , B, and A, as defined in the main text.
becomes large, and in turn, no large flavor conversion occurs; neutrinos remain pinned to their initial configuration.This prediction is again confirmed by Fig. 2, where we see a clear break to a region where flavor equipartition is not reached, and at sufficiently large energies (E ≳ 20 MeV), the ν x spectrum is hardly populated.

IV. DAMPED FAST FLAVOR PENDULUM
A. Setup of the problem A somewhat more complex situation arises when the fast dynamics, unperturbed by collisions, is not simply a precession.One obvious example is the fast flavor pendulum [18,44,45], which in its simplest manifestation arises in a homogeneous, anisotropic, azimuthally symmetric, and monoenergetic neutrino system.In the absence of collisions, if the angular distribution of the lepton number has crossings obeying a certain criterion [46], the system possesses a fast instability and, due to the presence of an infinite number of integrals of motion [19,44], exhibits a pendulum-like dynamics.
We now add a collisional damping term which is taken to be equal for neutrinos and antineutrinos, a case that was previously examined numerically [14] (see also Ref. [47] for a discussion of the impact of collisions on fast instabilities).We will see that our separation of fast/slow dynamics applied to this system directly recovers the earlier features.In this case, the polarization vectors P v and P v for neutrinos and antineutrinos depend on the velocity −1 < v < 1 along the axis of azimuthal symmetry.The dynamics is entirely driven by the lepton number where D 0 = v D v and D 1 = v vD v and e z the flavor direction.The EOM for D 0 is corresponding to damping of its component transverse to the flavor direction.
In the absence of collisions, no direction in flavor space is singled out and the conserved D 0 , given by initial conditions, defines the z-direction.Collisions with the background medium, for example beta processes, introduce a flavor direction and define e z .If the small seeds chosen to start the motion imply that D 0 is not exactly parallel to e z , this deviation is quickly damped but otherwise, this small misalignment has no further impact and we may assume that D 0 is conserved, simplifying the EOMs.In this case we may remove the term D 0 × D v in the EOMs by a uniform corotation.After these simplifications, are the EOMs to be solved.

B. Pendulum dynamics
In the absence of collisions, for a fast-unstable system with a single crossing, the motion of the polarization vectors is periodic and can be expressed in terms of three fundamental vectors only.The dynamics of these three vectors is identical to that of a spherical pendulum [18,19].In other words, we can introduce a fictitious system of three polarization vectors P α , with three velocities v α , obeying the EOMs with P 1 = α v α P α .The original polarization vectors can be now expressed in terms of this fictitious system using the connection [19] The two systems have identical dynamics, provided they obey the matching condition this condition needs only be satisfied at the initial time, and it will be automatically true at all later times.What is surprising is that, for the special choice of the isotropic collisional term we consider for this work, the equivalence between the full system D i and the threebeam system P α is rigorously true also in the presence of collisions.Indeed, by explicit substitution, we find that if P α obey the EOMs the associate polarization vectors D i automatically satisfy the correct EOMs.Therefore, without loss of generality we restrict our discussion to a system of three velocity beams only.The dynamics of three beams in the absence of collisions can always be reduced to a pendulum dynamics.We do this by introducing the variables with p = α v α .For Γ = 0, these indeed coincide with the EOMs of a spherical pendulum.

C. Slow dynamics evolution of flavor pendulum
The unperturbed pendulum motion, for Γ = 0, is characterized by the invariants of motion D 1 , J z = J • e z , α = J • D 1 /D 1 , and E = 1 2 J 2 + pD 0 • D 1 .Collisions lead to a slow evolution of these quantities.Therefore, just as in the case of CFIs, we can now obtain the SDEs describing only the slow change of these quantities over timescale ∼ Γ −1 , rather than the fast pendular motion over timescales ∼ µ −1 .As in the case of CFI, we need to average the right-hand side of the corresponding equations over many periods of the fast motion induced by the self-interaction term, which here is the pendular motion performed by D 1 .
From the EOMs, we derive the following evolution equations Jz = 0, (37a) where X = cos θ = D 1 •e z /D 1 is the deviation of the pendulum from the flavor direction and the overline denotes the average over a pendulum period.This set of equations can be simplified by introducing the quantity Φ = αD 0 D 1 p − J z E, whose derivative is obtained from the previous set as Since J z is conserved, this tells us that Φ relaxes to its asymptotic value exponentially.However, if the pendulum is initially in its upright position, as we always assume, and E = 1 2 J 2 z + pD 0 D 1 , and so Φ(0) = − 1 2 J 3 z .Therefore, Φ is conserved throughout the whole motion and is an adiabatic invariant, which means that at every instant.Finally we have only two variables, α and D 1 , because J z is constant and E is determined at each instant by the conserved Φ.To simplify the equations, we express everything in dimensionless variables a = α/J z and d 1 = D 1 /D 1 (0).Initially, both a and d 1 are equal to 1.We redefine time units such that Γ = 1.Now a and d 1 obey the equations From the energy conservation equations, we find Here cos θ min is the maximum depth of flavor conversion, namely the cosine of the minimum zenith angle reached by the pendulum in its unperturbed state, defined as Therefore, we define where the integral is taken only over the region where the square root argument is positive and −1 < X < 1; one can check that this implies the integration range where ξ = (1+cos θ min ) −1 .With this definition, we write This expresses the averages in Eqs.(40) explicitly in terms of a and d 1 .The dynamics of a and d 1 is now expressed in the form of a closed set of equations, where at every next instant the new pendulum parameters depend only on the pendulum parameters at the previous instant, and on the control parameter cos θ min .Once more, the SDEs cannot be solved analytically, but offer considerable insight about the solutions.First, no terms of order µ appear in the equations.Second, only two quantities remain, namely the length of D 1 and the projection of the angular momentum over the pendulum length α, which completely characterize the motion.Most of the original parameters of the problem have disappeared, including the collision rate Γ.This confirms the empirical findings of Ref. [14] that, provided Γ ≪ µ, the precise value of Γ only determines the evolutionary speed towards a stationary state, but not the final state itself.Similarly, the SDEs depend only on the single parameter cos θ min .This proves the empirical finding that the final state only depends on cos θ min .
The SDEs determine directly the evolution of the pendulum parameters, without need to solve the EOMs for all polarization vectors.The evolution from integrating the SDEs is shown in Fig. 3.The dynamics for a few selected values of cos θ min is shown in the left panel.Initially, all curves start from a(t = 0) = 1 (pendulum in the upright position) and d 1 (t = 0) = 1.For θ min very small, the pendulum is always damped to a vertical position with reduced length; in the final state, a → 1, implying that the pendulum ends up vertical.As θ min increases, the final length d 1 eventually reaches 0. As θ min increases further, the final outcome is to invert the direction, which is vertical but points downward compared to the original direction.This is signaled by a → −1 in the final state.The entire evolution can be compactly represented in the a-d 1 plane, as in the right panel of Fig. 3.All curves start from the upper right corner (a = 1, d 1 = 1), and their subsequent path depends on θ min .The final length d 1 (t → ∞) can be directly obtained, for each value of θ min , by identifying the intersection of each path with the axis a = 1 or a = −1, corresponding to the pendulum finally damped to the vertical position, either upwards or downwards respectively.By integrating the SDEs, we recover the linear relation between the final height of the pendulum and cos θ min that was found in Ref. [14], as we show in the inset of Fig. 3.
At each instant, the pendulum parameters can be used to recover the angular distribution of the (anti)neutrino density matrices.The practical procedure is outlined in Appendix D of Ref. [46].Because the damped pendulum is always pushed toward its vertical position during flavor evolution, the off-diagonal terms of the final-state density matrices always vanish, as expected.Nevertheless, the specific shape of the ELN angular distribution is in most cases non-trivial and specific scenarios need to be systematically explored in a full multi-angle framework; particular examples with ELN angular crossings are shown in Fig. 2 of Ref. [14].

V. CONCLUSIONS
We have presented a novel, approximate approach to solving the collective evolution of a homogeneous neutrino system that is subject to collisions with a back-ground medium.The key element is the separation of the fast dynamics implied by the large scale µ of neutrinoneutrino refraction and the slow dynamics caused by the collision rate Γ.The justification for the separation of scales is the weak-damping limit Γ ≪ µ that applies in all practical cases and means that between collisions, the system performs many oscillations, allowing one to use oscillation-averaged density matrices (or polarization vectors) to follow the slow evolution caused by collisions.One traditional example is the production of sterile neutrinos in the early universe or in SN cores, where the fast dynamics can be analytically integrated out.We have presented this case as a first warm-up exercise.
The novelty of our work consists in applying this method to much more complicated systems, where integrating out the fast dynamics leads to new slow-dynamics equations that we call SDEs, which themselves are nonlinear equations that require numerical solution, yet provide insights about the evolution and final state that are otherwise obscured by the fast dynamics.
The approach of separating fast from slow dynamics in a controlled way by averaging over the cycle of a fast periodic motion, in the mechanical domain similar to Kapitza's pendulum [48], differs from other attempts to remove unnecessary detail from the EOMs.In statistical mechanics, one removes the motion of the microscopic degrees of freedom and studies macroscopic quantities such as temperature or pressure.Coarse graining over the microscopic degrees of freedom of a fast flavor system may lead to analogous simplifications [49], but the underlying assumption of the system relaxing to thermal equilibrium remains for the moment a conjecture.
We have studied two nontrivial applications, the collisional flavor instability (CFI) and the fast flavor pendulum with an energy-independent collision term.There are three main advantages compared to a full solution: (i) A much coarser numerical step size.(ii) A smaller number of variables in the slow dynamics; this shows up most clearly in the fast flavor pendulum, where the large set of polarization vectors reduce to two coupled differential equations for the pendulum length and spin.(iii) Most importantly, the SDEs reveal certain aspects of the final state without actually solving them.
Notably in the case of CFIs, this approach allowed us to deduce the main features of the final distributions, namely equilibrium between the e and x flavor at low energies, and no flavor conversion at high energies.For the fast flavor pendulum, we proved that the final flavor composition depends only on the depth of the unperturbed pendulum, and not on Γ or any other property, a conclusion that had been previously reached empirically.
For the CFI, there is a subtle point about the scale of fast dynamics.Usually, the scale of the neutrino-neutrino interaction energy is estimated as µ = √ 2G F n ν .On the other hand, the relevant scale for the CFI derives from the length of the initial overall polarization vector and thus initially, µ CFI = √ 2G F (n νe − n νe − n νx + n νx ) is the relevant scale.Therefore, our method applies only when Γ ≪ µ CFI , whereas the opposite case µ CFI ≪ Γ is sometimes termed the resonance-like regime [9] in which a complete flavor swap may occur [50].Conceivably, even in this case, our method could be applied in a later phase when the z component of the global polarization vector has acquired a significant length.
In practice, our approach requires one to understand exactly the fast dynamics in the absence of collisions, so that the slow dynamics only touches those quantities that would be left invariant by the fast dynamics.This is why we had to restrict ourselves to homogeneous, azimuthally symmetric setups, for which the dynamics without collisions can be expressed analytically.However, if one were able to develop an understanding of the fast dynamics in inhomogeneous systems, our method could be applied to this more practically interesting case as well.
The impact of collisional flavor instabilities cannot be gauged by the limited information provided by linear stability analysis.A self-consistent solution should carefully incorporate matter-neutrino collisions and a multidimensional treatment of neutrino advection.Given the present difficulties in estimating flavor evolution in such complex environments, one possible way forward could be the approach presented in this work, where slowlychanging quantities can leverage some of the technical difficulties in the nonlinear regime.This new framework could pave the road toward understanding CFI in more complex systems that go beyond the common assumptions of isotropy and homogeneity.

FIG. 1 .
FIG.1.Comparison of the solution of the full EOMs (solid thin) and the solution of the SDEs (dashed thick) for the squared transverse component of the polarization vector corresponding to an energy E = 9.5 MeV.

1 FIG. 3 .
FIG. 3. Evolution of the pendulum parameters obtained from the SDEs.Left: Dynamical evolution of the pendulum spin along the axis a = J • D1/D1Jz and the dimensionless pendulum length d1 = D1/D1(0).Different colors correspond to different values of θmin, as identified in the right panel.Right: Evolution of the pendulum parameters in the a − d1 plane for various values of θmin; we show in gray a sampling of the paths for uniformly sampled values of θmin.In the inset, we show the final value of the pendulum height, ad1, as a function of cos θmin.

TABLE I .
Parameters used for the numerical simulation.