Fast Neutrino Flavor Conversion at Late Time

We study the fully nonlinear fast flavor evolution of neutrinos in 1+1 dimensions. Our numerical analysis shows that at late time the system reaches an approximately steady state. Using the steady state approximation we analytically show that the spatial variation of the polarization vectors is given by their precession around a common axis, which itself has a motion reminiscent of a gyroscopic pendulum. We then show that the steady state solution to the equations of motion cannot be separated in position and velocity, that is the motion is not collective in the usual sense. However, the fast evolution allows spectral-swap-like dynamics leading to partial decoherence over a range of velocities, constrained by conservation of lepton number(s). Finally, we numerically show that at late time the transverse components of the polarization vectors become randomly oriented at different spatial locations for any velocity mode and lepton asymmetry.


I. INTRODUCTION
Neutrinos emitted by stars present valuable opportunities to study neutrino properties [1]. While solar neutrinos have famously helped zero in on the large mixing angle scenario, neutrinos from supernovae may yet provide a unique opportunity to study neutrino-neutrino interactions -a crucial piece of the standard model of particle physics that has not been tested directly.
The rate of neutrino oscillations is typically dictated by the vacuum oscillation frequency, ω, and the matter potential, λ [2][3][4]. Until the early 2000s, it was believed that this paradigm was sufficient to describe neutrino oscillations inside supernovae as well [5]. At that time, the outstanding problem of the field appeared to be to understand the effect of fluctuations in the background matter density [6][7][8].
Following the pioneering papers by Pantaleone [9,10], however, it became clear that the issue is more subtle [11,12]. Owing to the large neutrino density, even free-streaming neutrinos experience significant forwardscattering off other neutrinos. Such scattering leads to a self-interaction potential, µ ω, that is proportional to the neutrino density and can dominate over the vacuum term. As a result, a gamut of new collective flavor transformations can occur inside supernovae.
The flavor evolution of a dense neutrino gas is governed by a large number of coupled nonlinear partial differential equations. These are almost always very difficult to solve. Although linearized stability analysis is useful to ascertain if and when fast conversion takes place, it cannot directly answer the question -what is the impact of fast flavor conversion on observable neutrino fluxes or the explosion mechanism? This is a significantly harder problem that requires understanding the nature of the solution in the nonlinear regime. A step in this direction was taken by Sen and one of the present authors [51], where the flavor evolution of a 4-beam model in 0+1 dimension was understood in the fully nonlinear regime.
In this work, we take another step in the same direction. We consider a dense neutrino gas in 1+1 dimensions, with a spectrum of velocity modes, and analyze the coupled flavor evolution of the neutrino system into the nonlinear regime. Our numerical analysis suggests that the system reaches an approximately steady state at late time. In the steady state approximation, we analytically show that the spatial variation of the polarization vectors is given by their precession around a gyrating flavor pendulum with a fixed length, spin, and energy, and the solution is not collective. The polarization vectors, when averaged over space, however, exhibit complete (partial) decoherence for zero (nonzero) lepton asymmetry. For partial decoherence, the non-vanishing range of velocity modes is dictated by conservation of lepton numbers. This kinematic decoherence stems from randomization of the transverse components. Numerical examples confirm these analytical insights.
The paper is structured as follows: Sec. II recollects the equations of motion, followed by our analytical results on the nature of the solution, conserved quantities, partial decoherence, and its dependence on lepton asymmetry. Sec. III has our numerical results and their comparison with analytical claims. We conclude in Sec. IV.

A. Equations of Motion
Neglecting momentum-changing collisions, the space and time evolution of two flavors of neutrinos with velocity v and vacuum oscillation frequency ω = |∆m 2 |/(2E) is given by [46,47], encodes the selfinteractions, with µ being the ν − ν potential. Note that we work in a basis {ê 1 ,ê 2 ,ê 3 }, where the Bloch vector for a ν e points along the directionê 3 .
In the fast flavor limit the vacuum and matter term in Eq.(1) are negligible compared to the neutrino potential term, so the solution for P ω,v becomes ω-independent. The self-term then enters the Hamiltonian only through the electron lepton number (ELN) distribution, i.e., the difference of occupation number densities integrated over energy, defined by G v = dω g ω,v . Thereby, allowing us to rewrite Eq.(1) in 1+1 D as We choose our units in Eq.(2) such that neutrino selfinteraction potential µ is 1, with length and time expressed in units of µ −1 . G v encodes the amount of lepton asymmetry of the system as a function of velocity modes.
One can obtain some more insight into the equations of motion, by expanding P v (x, t) for each velocity mode in terms of Legendre polynomials, L n (v). Using the expan- , and using the orthogonal property of Legendre polynomials, i.e., 1 −1 L r (v)L n (v)dv = 2 δ rn /(2r+1), following Ref. [29] one can rewrite Eq. (2) as

B. Steady State
We conjecture that at late time the system becomes approximately stationary in time, which says that in Eq. (2) we can drop the t-dependence from all the quantities. We will verify this conjecture in our numerical survey described in Sec. III.
In the steady state one can write Eq. (2) as and the equation for each multipole moment M r (x) as where From Eq.(5), for r = 1 and using 1n = (2n + 1)δ 0n , one gets that M 1 (x) is constant in space. This means that all the polarization vectors of different velocity modes precess about a fixed axis M 1 with the same frequency. This common motion can be hidden away by considering a rotating frame, using the transformation where we introduce the notation that all the quantities in the rotating frame are denoted by tilde, spatial derivatives in the rotating frame as d x , etc. Using the rotation formula described in Eq.(7), the equations of motion for the system at late time, described by Eqs.(4) and (5) in the rotating frame, look like Considering r = 0 and summing both sides of Eq.(8b) with weightage factor of 0r , from r = 0 to r = ∞, one can rewrite Eq.(8b) as two coupled equations: where, Similarly summing both sides of Eq.(8b), with a weightage factor of ∞ p=0 0p pr , from r = 0 to r = ∞, one obtains an equation for B(x): 0p pr rn M 0 (x) × M n (x) . (12) Eqs.(9a) and (9b) are nominally the same set of coupled equations which describe the "gyroscopic pendulum" (see Ref. [15]), with the difference that B(x), instead of being a constant, has an equation of motion in space described by Eq. (12). Nevertheless, one can follow Ref. [15] to derive an equation for a gyroscopic pendulum with a position-independent length, spin, and energy. Eq.(9a) clearly indicates M 0 , i.e., the length of M 0 (x) which is the same in rotating and non-rotating frames, is a constant in space. Taking dot products with D(x) or B(x), respectively, on both sides of Eq.(9a), and with M 0 (x) or D(x), respectively, on both sides of Eq.(9b), and with M 0 (x) for Eq.(12), one gets the following set of equations: Adding Eqs.(13a) and (13c) one obtains which implies that the "spin" i.e., position-independent and owing to the steady-state approximation also time-independent. Here, m 0 (x) = M 0 (x)/M 0 is the unit vector along M 0 (x). Similarly addition of Eqs.(13b), (13d), and (13e) reveals energy conservation: which implies Taking a cross product of Eq.(9a) with M 0 (x) one gets Dividing both sides of Eq.(18) with M 2 0 and then using Eq.(15), one can rewrite Eq.(18) as Differentiating Eq. (19) once, and using the spatial conservation of σ along with Eq.(9b), gives The vector M 0 (x) in flavor space plays the role of a gyroscopic pendulum. It has a fixed length at all spatial locations, so that it is restricted to move on a sphere of a fixed radius M 0 . According to Eq.(17) the energy of the pendulum is spatially invariant at late time, where the first term B(x) · M 0 (x) is equivalent to the potential energy of M 0 (x) in an inhomogeneous magnetic field B(x), and 1 2 D(x) · D(x) is the rotational energy of the system, with D(x) playing the role of the orbital angular momentum of the system. Eq.(15) describes one more conserved quantity σ, which says that the component of the angular momentum D(x), that is parallel to M 0 (x), is equivalent to the pendulum's spin and is spatially constant.
The above analysis shows that at late time the spatial structure of the solution is very simple: every P v has a precession about an axis m 0 (x) with a frequency M0 v , where m 0 (x) itself has a motion equivalent to that of a gyroscopic pendulum in an inhomogeneous magnetic field, keeping a fixed length, spin and energy. This similarity to the gyroscopic pendulum solution, in a corotating frame, is however limited to the conserved quantities and a formal similarity of the equations of motion. The actual motion is different, because B in this case has a nontrivial motion that is not necessarily much slower than that of the M n . These conserved quantities are formally identical to those in Ref. [65], but they have a slightly different interpretation. Unlike in the above, where exactly 0+1D or 1+0D was considered, ours are obtained under a steady state approximation in 1+1D and need not be exactly conserved.

Nature of the Solution
We now prove that the solution at late time, described by Eq.(8a), cannot be separable in position and velocity coordinates, i.e., it cannot be written in the form: In our notation, the components of S v in the rotated frame look like is not possible. Otherwise, the normalization of S v for a fixed velocity mode will be different at different points in space, which is obviously unphysical. By a similar argu- cannot be possible either, as that will lead to S v having different normalizations for different velocity modes at a fixed point in space. Plugging Eq.(21) into Eq.(8a), we get two separate sets of equations -one governing the spatial dependence and the other with the velocity dependence of the full solution. The equations governing the velocity dependence look like where, This already suggests that, again, no meaningful solutions exist. However, we should confirm that the spatial solutions do not diverge. The equations for the spatial dependence look like To solve the above set of equations, we multiply Eqs.(24a), (24b), and (24c) by f 1 (x), f 2 (x), and f 3 (x), respectively, which gives d dx where C 1 and C 2 are integration constants. The choices C 1 = C 2 = 0 are disallowed, as that will mimic the case Solving Eq.(27) one gets is an elliptic integral of first kind defined as and C 3 being an integration constant. One has to basically invert Eq. (28) to get the behavior of f 3 (x) as a function of x. This function is not infinite everywhere and concludes our analytical proof that there are no separable solutions in steady state. Still, this may be opaque, and to give a flavor for the solution, consider a special case with C 1 = C 2 , and the solution for Eq.(27) then becomes For any value of C 1 and C 3 , the solutions for f 1 (x), f 2 (x), and f 3 (x), are finite -with either oscillatory or constant behavior over all space. This implies that, for all the modes with v = 0 the solution for all the components of S v (x) will be exactly zero, giving rise to an unphysical solution. Another qualitative way to understand this is from Eq.(4), where one can see that due to the spatial dependence of the velocity dependent term, , in the right hand side of the equation, it can not be rotated away by a rotation of P v (x) in the position space. So, the equation of motion governing the spatial behavior for each P v (x) cannot be same for every velocity mode v, indicating a non-separable solution in x and v.

C. Approach to Steady State
In this section, we consider the approach to steady state, without making the steady state approximation. Rather, we investigate how the spatially averaged polarization vectors rearrange themselves at late time, conserving lepton number and giving rise to flavor depolarization over a range of velocities.

Dependence on Lepton Asymmetry
Polarization vectors cannot be completely depolarized if the lepton asymmetry is non-vanishing. To see this, we integrate both sides of Eq.(2) over all velocity modes to get Performing a spatial average on both sides of Eq.(31) over the entire box with the periodic boundary condition, This implies where M 0 (t) is understood as Using the approximate stationarity at late time, one can argue, from Eq. (31), that M 1 (x, t) is approximately constant all over space. The constant in Eq.(33) can be determined from the initial condition of the system, where all neutrinos are emitted as approximately flavor pure states from every point in space, i.e., s v (x, 0) = 1 for all x and v. One can write M 0 (t) as which remains true at all times. The above initial conditions, and defining gives rise to a set of three conditions that have to be satisfied even in the nonlinear regime: One way that Eqs.(36a) and (36b) can be satisfied at late time for systems with zero lepton asymmetry, i.e., = 0 for every velocity mode v. This indicates a depolarization of the system. However, for systems with A = 0 such a simple solution can exist if s (3) v = 1 for every v. This is possible only for an inert system where fast flavor conversion do not occur. Therefore, fast flavor conversion in the nonlinear regime can happen for A = 0, if different velocity modes s behave in such a velocity-dependent way that it preserves the value of the lepton asymmetry A. We will find, numerically, in case of A being +ve (or −ve), s will remain much closer to 1 for modes for which G v is +ve (or −ve).

Fast Depolarization
Now we will show that a set of polarization vectors can undergo a "fast depolarization" -i.e., the system exhibits a "fast spectral split" like evolution, similar to what appears in the foundational study by Raffelt and Smirnov [17,18], that leads to s (3) v (x, t) flipping (and approximately vanishing) over a range of velocities, but constrained by conserved lepton numbers. Consider Eq. (3). After averaging over space, we get The averaging is over the cross-product of the vectors but, because the polarization vectors vary very fast over space, one expects that the averaging factorizes over the cross product, which gives Integrating the above equation over all velocities, one recovers that M 0 is time-independent. However, multiplying the above equation with v and then integrating over v gives where we define a new vector, Eq. (38) ensures that each P v (t) precesses around its Hamiltonian, Initially all polarization vectors are either aligned or antialigned toê 3 , and over time M 1 (t) has dynamics, while the polarization vectors are dragged by H v (t) .
We assume that the polarization vectors stay close to their Hamiltonians, and thus remain in the plane formed by M 0 (t) and M 1 (t) , ignoring the precession around their Hamiltonian. The vector X(t) can be decomposed as which, upon insertion in Eq.(39), gives (43) Note that X(t) typically has dynamics at the same frequency as M 1 (t) . So M 1 (t) , that has the dynamics of a gyroscopic pendulum, can not only precess around M 0 (t) but can also have bipolar nutations.
In a frame that co-rotates with the plane formed by M 0 (t) and M 1 (t) , the Hamiltonian shifts by −(1 + α(t)) M 0 (t) due to pure precession of M 1 (t) . Focussing on theê 3 component of the Hamiltonian, we find that The components of M 0 (t) and M 1 (t) alongê 3 remain conserved if M 1 (t) has no nutation. In any case, in addition to Eq.(36b), we can define another useful lepton number,

M
The above stays constant while M 1 (t) precesses, but may flips to a different value if a bipolar instability is triggered instead. Theê 3 component of the co-rotating Hamiltonian is therefore, Depending on how α(t) and B(t) vary with time, the above can change its sign between some initial time and final time. The condition for such a sign-flip to occur, is where (. . .) ini, fin are at the beginning and end of the evolution. If the Hamiltonian for a velocity mode v can change sign by fulfilling Eq. (47), and obeying the constraints on A and B, then that velocity mode following its Hamiltonian may flip its sign as well. So far, our discussion closely follows the discussion in Refs. [17,18]. There are a few subtle differences. Unlike in the case of bipolar swaps where only the lepton number A needs to remain conserved, here one has two constraint equations, i.e., conservation of A and conservation or flip of B, depending on whether nutations occur. Further, this derivation was for spatially averaged polarization vectors, P v (t) , that do not necessarily maintain the same length as their unaveraged counterparts P v (t). This point is quite crucial, as the astute reader will notice, otherwise there is no obvious source of irreversibility: unlike in the bipolar spectral swap where that is provided by a decreasing µ, the irreversibility is provided by the relative dephasing, i.e., kinematic decoherence, of the polarization vectors, which is responsible for the irreversibility of α(t). Further, with simple choices of G v , the swapping function of the s (3) v cannot be a mere sign-flip over a block of velocities, in general. Such a sign-change across a crossing can preserve A, but unless this block of G v is antisymmetric in v such a flip doesn't obey the constraint on B.
More detailed exploration of fast depolarization will be published separately, but here we note a main qualitative feature: over time the s for a range of velocities become close to zero, even flipping their sign, but in the remaining range stay close to their initial state s Consider a simple case, where A > 0, B > 0 and G v has a single crossing at zero. It is possible that the modes with v ≈ +1 do not flip sign or become small, whereas those with smaller v either become very small or flip their sign (with change in magnitude, as well). The converse however is not possible because such a configuration cannot preserve A and B. On the other hand, when A < 0 and B < 0, the modes with v ≈ −1 cannot flip sign. We will see this pattern in our numerical calculations.

III. NUMERICAL STUDY
In the following subsections, we now discuss our own numerical strategy to solve the equations of motion, the late-time behavior of the solution, and their relation to the analytical claims made in Sec. II.
We developed a code for solving Eq.(2). In our code, we discretize each of the spatial directions as well as the velocities, considering N x spatial modes and N v velocity modes, to get a total of 3 N x N v coupled nonlinear ODEs, where the factor of 3 comes due to the 3 elements of each polarization vector. It solves these coupled nonlinear ODEs as a function of time using python's zvode solver, which is a complex-valued variable-coefficient ordinary differential equation solver in python which implements Backward Differentiation Formula for doing numerical integration. The spatial derivatives at each spatial point are computed using python's scipy.fftpack.diff package, which uses the Fast Fourier Transform method for calculating derivatives.
For the computations shown in this paper, we choose a periodic boundary condition in space such that, for any time t and any velocity mode v, our solution satisfies We choose initial conditions such that all the neutrinos, with any velocity, are emitted as purely electron flavored states at every point in space. To trigger the flavor evolution, a perturbation of 10 −6 is used as an initial seed for the transverse components of the polarization vectors, for every velocity mode at the midpoint of the 1D box. This choice is arbitrary, but motivated by trying to seed all wavelength modes of the instability equally, as opposed to using, say, homogeneous initial conditions that favor a specific mode. We ensure that the size of the box L and the maximum time t final up to which we solve the equations, are such that 2 t final < L, as a result of which the mode with the largest velocity does not see the boundary of the box when emitted from the centre of the box, where we seed the instability. Specifically, we choose a 1D box of size L = 115 and discretize it into 2 12 spatial bins to solve Eq.(2) numerically upto a time of t final = 50. We discretize the ELN distributions into 2 7 velocity modes, giving rise to a total of 3 × 2 12 × 2 7 = 1572864 coupled ODEs. These choices are optimized to obtain sufficient accuracy and precision, as we show in Appendix A. To show the dependence of our results on the ELN, we consider three qualitatively different continuous ELN distributions shown in Fig. 1. We choose these ELNs in such a way that they all have a zero crossing in their velocity distribution for fast flavor conversion to occur, but they have lepton asymmetries that are either zero, or negative, or positive. Fig. 2 clearly indicates that the system becomes approximately steady already at t NL ≈ 10, and to ensure that we were deeply in the nonlinear regime, we ran our code until t = 50.
We define some notation for convenience: where φ v (x, t) encodes the correlation between the two transverse components of the polarization vector for each velocity mode, with S ⊥ v (x, t) describing the length of the vector in theê 1 −ê 2 plane. We will denote all the spatially averaged quantities over the whole box of size L as · · · , the velocity averaged quantities by · · ·, the probability distribution for the spatial distribution of some quantity as PDF(· · · ) and the standard deviation for that distribution as · · · .

A. Late Time Behavior of Parallel and Perpendicular Components
In Fig. 3, we show the time evolution of the spatial average and standard deviation of s (3) v (x, t), for all three types of lepton asymmetry. In Fig. 4, we show the histograms of s (3) v (x, t) and |S ⊥ v (x, t)|, at a time t = 50 and for various velocity modes color-coded by their velocity. We have checked that these histograms are themselves quasi-stationary at late time. These results indicate approximate stationarity in time in the nonlinear regime for all velocities. Note how s (3) (t) separates into two cohorts, dictated essentially by the sign of v, at late time in the cases A > 0 and A < 0. This behavior is very similar to spectral swaps. As explained in Sec. II, which modes remain close to their initial state and which move away depends on A, B, and the higher moments.  (1) v (x, t) s (2) v (x, t) s (3) v (  (1) v (x, t) s (2) v (x, t) s (3) v (  (1) v (x, t) s (2) v (x, t) s (3) v (x, t) ⟨s (1) v ⟩ ⟨s (2) v ⟩ ⟨s (3) v ⟩ ⟨s (1) v ⟩ ⟨s (2) v ⟩ ⟨s (3) v ⟩ ⟨s (1) v ⟩ ⟨s (2) v ⟩ ⟨s ( The spatial distribution of s (3) v , i.e., PDF(s v ), and the distribution of |S ⊥ v |, i.e., PDF(|S ⊥ v |), shown in Fig. 4, are obviously related because |S v | = 1. The latter is obtained from the former by multiplying with the Jacobian of the transformation from s Consider the case with A = 0. In this case, PDF(s v ) is approximately uniform (see top left panel of Fig. 4) and as a result, PDF(|S z ⊥ v |) should be a sharply peaked distribution for every v with a peak very close to |S ⊥ v | = 1 and a tail at |S ⊥ v | = 0. In bottom left panel of Fig. 4, we recover this unsurprising feature. For cases with A < 0 (A > 0), the negative (positive) velocity modes show a peaked distribution in s (3) v almost close to 1, implying a more broadened distribution for |S ⊥ v | between 0 to 1, as seen in the middle and right panels.
Top panels of Fig. 5 show the oscillatory nature of the spatial variation of s (3) v (x, t), at time t = 50 when the system becomes fully nonlinear. The bottom panels of Fig. 5, show that the spatially averaged polarization vectors are depolarized for A = 0, with s v ≈ 0, but not so in case of nonzero lepton asymmetry. For A = 0, the polarization vectors have not flipped close to v ≈ −1 (v ≈ +1) for A < 0 (A > 0), as predicted by our analysis of fast spectral swaps. One interesting point to note from the bottom-middle and bottom-right panels of Fig. 5, with spectra A = 0, is that the functional behavior in the velocity space shows a mirror symmetry when we switch the sign of A. This is not surprising, because our chosen ELNs have a symmetry G A<0 (v) = −G A>0 (−v).

B. Behavior of the Phase φv(x, t)
We obtain the phase φ v (x, t), for a given space-time point (x, t) with a velocity v, in the range [−π, π] using the python package math.atan2 which takes s (1) v (x, t) and s (2) v (x, t) as inputs, obtained from our numerical simulation, and returns the principal value of arg S ⊥ v (x, t) as the output. Fig. 6 shows that the phase is almost uniformly random, and its mean and variance are consistent with that expected of a uniformly random distribution over space, for every velocity mode and for any type of lepton asymmetry. Initially, for each velocity mode the transverse components of a polarization vector S v (x, t) are corre- lated, but in the extreme nonlinear regime such a kinematic phase coherence is lost.

C. Conserved Quantities
In Sec. II we showed that M 1 is almost constant in space in the steady state approximation, which allowed us to go in a rotating frame where M 0 has a motion equivalent to a gyroscopic pendulum. Further, the pendulum's length M 0 , spin σ, and energy E can be seen to be spatially conserved. In this subsection, we undertake a numerical survey to verify if these quantities indeed satisfy the above analytical claim.
To calculate the different multipole moments numerically we used the late-time solution for P v (x, t) from our code and integrated over all the velocity modes from -1 to 1 at each spatial point with an appropriate weightage factor, corresponding to the Legendre polynomial of a degree equivalent to the multipole number of that particular moment. One technical issue is that the quantities σ and E involve infinite series in the multipole moments which are numerically impossible to calculate. Thankfully, the spatially averaged components of the polarization vectors, shown in the bottom panel of Fig. 5, behave very smoothly as a function of velocity and can be adequately approximated by the first few multipole moments. Fig. 7 shows the numerical results for the spatial variation of the conserved quantities at late time. Most of the conserved quantities are approximately constant over space, but the energy E of the gyroscopic pendulum does not appear to be sufficiently constant (relative to, say, the variation of individual polarization vectors). Even at late time the stationarity of the components of S v (x, t) is approximate, which reflects in the degree of variation in the spatial behavior of the conserved quantities. In con- trast, if one were to look at the velocity averaged P(x, t), that is much noisier even at late times, not unlike what is shown in the upper panels of Fig. 5.
To verify if Eq.(2) has a separable solution for the components of S v (x, t) in position and velocity space we define a ratio as R for the i th component of the polarization vector at a time t in the full nonlinear regime for two different velocity modes v 1 , v 2 . This will be spatially constant if the solution is separable in x and v. To calculate this quantity we considered the solution at t = 50 with the modes v 1 = 0.5 and v 2 = 1 for all the three components (i = 1, 2, 3) of the polarization vector, and plotted R (i) v1,v2 (x, t) as a function of x. Fig. 8 shows a large deviation from a constant, confirming a non-separable solution in the nonlinear regime.

IV. CONCLUSIONS
Our goal in this paper was to study the late-time behavior of fast oscillations. To this end, we considered a model in 1+1D and gave an analytical as well as a numerical understanding of its flavor dynamics. Here we summarize our main results: • The system reaches an approximately steady state in time in the extreme nonlinear regime.
• In a rotating frame the spatial evolution of all the polarization vectors is a velocity-dependent precession about a common axis.
• This common axis acts as a gyroscopic pendulum with a fixed length, spin, and energy leading to an oscillatory behavior in position space.
• The steady state solution for the components of the polarization vector is not separable in position and momentum space.
• The phase of the transverse component of the polarization vector at different spatial locations becomes randomly distributed over the interval [−π, π] at late time, for all velocity modes and with any value of lepton asymmetry.
Abs max(|P|-1) FIG. 9. Left: Precision of our calculation, estimated by convergence of the computed | P ⊥ | for different discretizations Nx ×Nv (to be read as the exponents of 2, i.e., 12 × 7 refers to Nx = 2 12 and Nv = 2 12 ), as a function of time t, shown for the A < 0 case. Right: Accuracy of our calculation, estimated by maximum departure (over all space) of the magnitude of Pv=0.5 from unity, as a function of time t, shown for the A < 0 case and different discretizations. Absolute and relative tolerances were taken to be same and mentioned alongside.
• The velocity dependence of the components of the polarization vector is controlled by the lepton asymmetry of the system. Systems with zero lepton asymmetry have a "decoherent" behavior for all velocity modes, with s = 0 for every v. In case of nonzero lepton asymmetry, the decoherence is not complete, but s show some dependence on velocity to conserve the lepton asymmetry.
• The behavior of s (3) v in the velocity space shows a mirror symmetry if the sign of the lepton asymmetry is flipped.
• s (3) v change with time in a way reminiscent of spectral swaps. At late time, the configuration becomes steady and is given by a velocity dependent "swap" function. The constraints on A and B must be obeyed for such a "swap".
We hope that these results provide some insight into the late-time flavor dynamics associated with fast flavor conversions of self-interacting neutrinos. Although we have chosen to study a simple system in 1+1D, some of these physics results may be useful to understand the more realistic scenario associated with neutrinos in a core-collapse supernova.

Appendix A: Error Estimate
In Fig. 9 (left panel), we illustrate the precision to be expected of our calculation. We solved the problem, outlined in the main text, for a sequence of increasingly finer discretizations for space (N x ) and velocity (N v ), with the log 2 of the numbers of divisions noted in the legend of the figure. We checked for convergence by comparing the results for these discretizations with the discretization N x = 2 12 and N v = 2 8 . Our results indicate that a discretization of N x = 2 12 and N v = 2 7 is at most O(10 −8 ) off from yet finer discretizations.
In Fig. 9 (right panel), we show the accuracy expected of our calculation. We check if the length of the polarization vectors remain fixed at 1. The error we incur on this is a lower bound on the error in our calculations. We find that our chosen discretization, N x = 2 12 and N v = 2 7 , does as well as finer discretizations, making an error of O(10 −4 ) at late time. Considering that the quantities we are interested in are close to O(1), this error is tolerable. Also, we find that the tolerance dictates how fast the errors increase initially, but the errors plateau out at roughly ∼ 100 times the tolerance. Whereas, if one is interested in a solution that is accurate to O(10 −4 ), it is less expensive to raise the tolerance to 10 −6 , without significant penalty on accuracy and with significant gain in speed.