Elaborating the Ultimate Fate of Fast Collective Neutrino Flavor Oscillations

Dense clouds of neutrinos and antineutrinos can exhibit fast collective flavor oscillations. Previously, in Phys. Rev. Lett. 126 (2021) 061302, we proposed that such flavor oscillations lead to depolarization, i.e., an irreversible mixing of the flavors, whose extent depends on the initial momentum distributions of the different flavors. In this paper, we elaborate and extend this proposal, and compare it with related results in the literature. We present an accurate analytical estimate for the lower resting point of the fast flavor pendulum and underline the relaxation mechanisms, i.e., transverse relaxation, multipole cascade, and mixing of flavor-waves, that cause it to settle down. We estimate the extent of depolarization, its dependence on momentum and net lepton asymmetry, and its generalization to three flavors. Finally, we prescribe approximate analytical recipes for the depolarized distributions and fluxes that can be used in supernova/nucleosynthesis simulations and supernova neutrino phenomenology.

Neutrinos change their flavor with a time-periodic probability, e.g., sin 2 2θ sin 2 [ωt/2] in vacuum, due to quantum interference of two eigenstates evolving with a frequency difference ω = ∆m 2 /(2E) [1][2][3]. In ordinary matter, forward-scattering off the background particles modifies the mixing angle θ and the oscillation * soumyaquanta@gmail.com † bdasgupta@theory.tifr.res.in rate ω [4,5]. However, often one is interested in the flavor composition after a sufficiently long time, when the flavor conversion probability is found to become timeindependent, e.g., 1 2 sin 2 2θ for averaged oscillations in vaccum [6][7][8]. Generally this is because of decoherence, which can occur in two ways: Either the flavor evolution of an individual neutrino becomes stochastic, e.g., due to collisions [9,10]; Or due to observational limitations -such as spatial, temporal, or energy resolution -which result in a pooling together of many neutrinos, each with a slightly different relative phase between its two components [11]. See Ref. [12] for a clear exposition. It should be noted that the nature and extent of the late-time neutrino mixing, even after the oscillations have ceased, can have nontrivial dependence on energy/momentum and can encode information about the system [6][7][8].
Flavor oscillations of neutrinos from dense astrophysical sources, e.g., deep in a supernova, merging neutron stars, or the early Universe, exhibit an additional novelty. These neutrinos can frequently forward-scatter off other oscillating neutrinos, leading to novel collective flavor oscillations [13]. The effect depends on the neutrinoneutrino forward-scattering rate µ = √ 2G F n ν [14,15], which typically exceeds the average oscillation rate ω in these environments. Under its influence, neutrinos can collectively oscillate at a synchronized rate ω [16], or the bipolar/slow rate µ ω [17], or the fast rate µ [18]. Remarkably, collective oscillations are predicted to occur with large amplitudes even for the matter-suppressed mixing angles expected in dense regions [19][20][21][22].
Collective oscillations display a rich phenomenology, but most remarkably they can lead to novel signatures of flavor mixing at late times. For slow collective effects a prominent signature is a set of energy-dependent swaps between the flavor spectra [23][24][25], with subleading decoherence effects [26,27]. For fast oscillations, the signature is less clear, but it is plausible that the decoherence effects are more important. There are two noteworthy issues: Firstly, one cannot therefore derive the late-time decoherent limit by straight-forwardly averaging over a known coherent oscillation probability. This is because all neu-arXiv:2205.05129v3 [hep-ph] 27 Dec 2022 trinos evolve interdependently, with unusually weak dependence on both θ and ω, and an analytical solution is not available in general; see however Ref. [28]; Secondly, fast oscillations can occur very quickly. While slow instabilities develop over a few 100 km, or more, fast oscillations and their associated decoherence effects can occur over much smaller distances ∼ 10 −4 m. Thus, their impact can be important already inside the star. E.g., stellar heating and nucleosynthesis could be affected [29,30]. See Refs. [31][32][33][34][35][36][37] for studies in this direction. As a result, it is both challenging and important to understand the late-time behavior of fast oscillations.
Starting with the first explorations of fast flavor oscillations in the nonlinear regime [38], it was seen that the survival probability eventually stops oscillating, and instead approaches a quasi-steady state [39]. The phase space distributions of the different flavors get irreversibly mixed [40,41]. We note that this is because of dephasing, and not because of collisions which help kick-start but do not overwhelm fast oscillations [42]. To emphasize this distinction, we denote this as depolarization [40,41]. The moniker is borrowed from optics, where it is used to refer to the shrinking of the polarization sphere (i.e., the Stokes parameters Q, U, V get smaller) without dissipation (i.e., loss in intensity I). We use it also to draw attention to the novel associated flavor conversion -full/partial equilibration of the flavor spectra, depending on velocity, while conserving the lepton asymmetry.
In two previous papers, hereafter B20a [40], and B20b [41], we have explored this phenomenon in detail. The purpose of this work is to elaborate these results, and to compare them with several closely related works. First, we compare our predictions for the behavior of the so-called fast flavor pendulum, with those by Johns et al. [43], hereafter J20, and the recent work by Padilla-Gay et al. [44], hereafter PG21. Then we compare our depolarization proposal with work by Wu et al., hereafter Wu21 [45], and by Richers et al., hereafter R21a [46] and R21b [47], respectively, which contain closely related results. We also compare and contrast our results with those by Martin et al., hereafter M20 [48] and M21 [49], where they do not find a depolarized steady state. Although our study is not intended to supplant a systematic code-comparison, the comparisons provided here should clarify a number of conceptual issues. See Refs. [50][51][52][53][54][55][56] for related studies of fast oscillations in the nonlinear regime. A separate body of work has focussed on the initial growth of fast instabilities; see Refs. [57][58][59][60][61][62][63][64][65][66][67][68].
This paper is structured as follows: We outline our setup in Sec. II. Sec. III presents an accurate estimate for the lower resting point of the fast flavor pendulum, and analyses of transverse relaxation, cascading of multipoles, and mixing of flavor waves. Sec. IV gives an estimate of the extent of depolarization and its generalization to three flavors. Sec. V contains recipes for the depolarized distributions and fluxes in a form that is usable for supernova simulations or neutrino phenomenology. Finally, in Sec. VI, we conclude with a summary and outlook.

II. FRAMEWORK AND METHODS
We use natural units throughout, with = c = 1. In each phase space cell d 3 p d 3 x, one has [69] i where ρ p is the matrix of densities and H p is the flavor Hamiltonian matrix. The phase space cells are taken to be sufficiently large, so that position and momentum can be simultaneously specified [70]. We ignore momentum changing collisions, external forces, and neutrino mass-mixing [71], which are typically negligible on timescales of the fastest neutrino oscillations. The velocities v = p/|p| and energies E = |p| do not change and serve as immutable labels. The range of E is from −∞ to +∞, to include antineutrinos of energy E by writing them as if they were neutrinos of energy −E. Axisymmetry restricts that the flavor evolution depends on a single spatial coordinate z, a single momentum coordinate v, and of course on time. This is a simple model for neutrino flavor evolution in a supernova, after it starts free streaming. Under these assumptions, the flavor content encoded in each ρ v evolves as [40,41] Here S v is the Bloch vector encoding the flavor state for neutrino modes with velocity v, with |v| < 1. We denote flavor space vectors by sans-serif letters, e.g., S, and the components parallel to theê 3 direction by (. . .) . The transverse vector confined to theê 1 −ê 2 plane for any flavor space vector, for e.g., S, is defined through the following vector formula Magnitudes are shown in the usual font, e.g., S = |S|.
The ELN distribution function G v is the excess of the phase space distribution of ν e over ν µ (andν µ overν e ), integrated over E 2 dE and divided by a typical density, say n ν . Only the product of µ 0 and G v appears; though, one defines a rate µ 0 ∝ G F n ν as the collective potential. Hereafter, we set µ 0 = 1, and express z and t in units of µ −1 0 . The ELN becomes dimensionless in these units. For this paper, we will mainly consider three families of ELN distributions, shown in Fig.1: G B v is a 'Box' spectrum, piecewise constant in v on either side of the crossing at v = 0. G L v is piecewise 'Linear', with an adjustable crossing at v c . Similarly, G C v is 'Cubic'. In addition, we will also study the ELNs G 3a v and G 4b v , as defined in Refs. [48,49], as well as the ELN in Ref. [45], to compare our results. For all the ELNs, the lepton asymmetry is denoted by A = dv G v . All these ELNs are inspired by SN ELNs with a single crossing, where ν e dominates overν e in the forward direction v > 0 and vice versa, and we restrict our study to A > 0. We study the dependence of the flavor state on lepton asymmetry A and on the crossing velocity v c .
In all our numerical computations, we solve Eq. (2) with initial conditions that all neutrinos -with any velocity v and at all locations in the periodic one dimensional box of length L -are emitted in the electron flavor, i.e., S v [t = 0, z] = +ê 3 . The numerical set up, i.e., discretizations, dimensionality, tolerances, etc., have been kept exactly the same as in B20b [41].
To start the flavor evolution, we supply tiny initial perturbations to the transverse components of the Bloch vectors. These are referred to as seeds, and are a numerically efficient means of initiating the flavor evolution. In reality, neutrino mass terms would provide the initial misalignment from a pure flavor state, but as we have set them to zero in the fast oscillation limit we resort to this numerical alternative. See the Supplemental Material of B20b [41] for details, including a discussion of the dependence on seeds. Unless stated otherwise, we will assume a spatially extended seed with transverse perturbations of amplitude 10 −6 and random relative phases.

III. FLAVOR PENDULUM AND RELAXATION
as the n th moment of the Bloch vector S v , with L n [v] being the n th Legendre polynomial in v, we can rewrite Eq. (2) in multipole space as [26,72], where Using periodic boundary conditions, and the approximation that spatial averaging factorizes over the dot and cross products of vectors, one can write the spatially averaged or coarse-grained version of Eq.(5) as For brevity, hereon we will mostly omit writing ... for the spatially averaged quantities. Instead, when we occasionally need to refer to quantities which are not spatially averaged, we will explicitly show the z-dependence, e.g., , as opposed to the averaged version S v . Hopefully, the distinction will also be clear from the context. Eq. (7) represents an infinite tower of equations. We will truncate this tower beyond n = 3, effectively assuming that the n ≥ 4 multipoles are negligible. This gives a set of four coupled ODEs: where Note that Eq.(8a)-(8d) are written in a frame rotating around M 0 with a frequency √ M 0 · M 0 , so that the common rotation of all M n around the axis M 0 , encapsulated in the second term on the left side of Eq.(7), is undone.
Eq.(8a) is the usual lepton number conservation which gives M 0 = constant, whose 3rd component is which is the lepton asymmetry. Eq.(8b) and Eq.(8c) are similar to Eq.(7) in Ref. [22], and can be combined to get which is the familiar pendulum equation for M 1 with length √ M 1 · M 1 . However, the vector B that acts akin to gravity is not a constant, and instead obeys Eq.(8d). We remind that M 1 is spatially averaged. Fig.2 shows the numerical solution of Eq.(2) for the parallel component of M 1 . The first thing to note is that it does not continue to oscillate forever. Rather, it comes to rest after a few cycles of oscillations. The late-time resting point depends on the lepton asymmetry (A), zero crossing position (v c ), and the nature of ELN. Note also that relaxation leads to a lower final resting point than the lower turning point of the first few oscillations, especially for the smaller values of A.

A. Resting point for M1
We now compute the resting point of the M 1 pendulum starting from the equations of motion. We will not assume that the moment vectors have constant lengths, and instead assume that certain phases randomize. In doing so, our approach departs from J20 [43], or the more recent PG21 [44], where spatial dependence and relaxation are absent. Thus, rather than deriving the lower turning point of the periodic M 1 pendulum, we focus on deriving the lower resting point of the relaxed M 1 pendulum.
According to Eqs.(8a)-(8d) the energy and spin of M 1 pendulum are which are conserved quantities in time t. The motion of B in Eq.(8d) allows us to write two more conserved quantities in time, i.e., We give a name to the parallel component of the pendulum vector M 1 : We also use some temporary shorthand notation to eliminate excess clutter in the derivation to follow, until Eq.(21d): The quantities in Eqs.(13)-(14c) will be denoted with subscripts, i at the initial time t = 0, and f at the late time when the system becomes steady. We do not use subscripts for quantities that are constant in time, e.g., for b and k.
Our aim is to derive m f , i.e., theê 3 component of the steady-state relaxed M 1 pendulum. The key idea is to use the steady state condition and to eliminate any unknown late-time perpendicular components in terms of conserved quantities.
At the resting point of M 1 one must have Note that this resting point, defined above, allows for relaxed solutions arising due to assumption of dephasing in deriving our approximate Eq.(7) from the exact Eq.(5).
In contrast, the turning point of D 1 (as given in Eq.14 of J20 or Eq.13 of PG21) explicitly excludes the spatial dependence of D 1 . Equations (15) and (8b) imply D f = (0, 0, u f ). This along with the conservation of σ between the initial and final positions of M 1 pendulum implies Equation (16) and the conservation of E and B 2 2 + K · D predict that At the resting position of the M 1 pendulum, ∂ t m| f ≈ 0, and from Eq.(8b) one has D ⊥ f ≈ 0. Then, using the conservation of σ per Eq.(11b) implies Using Eqs.(17a) and (17b), and ignoring the trivial solution m f = m i , help to simplify Eq. (19) in terms of the desired variable m f : The solutions of Eq. (20) are Only m +− f in Eq.(21c) has the correct qualitative behavior with A and v c to qualify as a solution. Note that m i , u i , b, and k, and thus m +− f , are known from the ELN. In Fig. 3 we plot m +− f , as obtained from Eq.(21c), in In all these plots, the blue disks  The qualitative dependence of the resting point on A or v c can be understood as follows: The kinetic energy of the M 1 pendulum is Clearly σ is a constant and can be determined from the initial conditions. For example if we consider For the case with A = 0 and v c = 0, one has σ ∼ A/2, whereas for A = 0.8, v c = 0, one has σ ∼ 0.4 (1 + v c ). As a result, in the limit A → 0 or v c → −1, we have σ → 0 so that the impact of internal spin σ is small. Thus, the pendulum swings like an ordinary pendulum resting at a smaller m f . In the other limit when σ is large one can approximately neglect the M 1 × ∂ 2 t M 1 term in Eq.(10) so that the M 1 pendulum equation becomes a simple spinprecession equation indicating m f ≈ m i . This is roughly the case with large A or v c → +1.
The resting point need not coincide with the lower turning point given in Eq. (14) of J20 [43] (the factor of 9 therein should be 5/4 -which has been corrected here). For the ELNs G B v and G L v , which do not have a cubic term, one does not expect a sensible estimate from Eq.(14) of J20 (since it formally diverges for D 3 = 0). However, even for the ELNs where a cubic term is present, and one ought to get a sensible estimate, we see that the turning point is not the resting point (gray line marked as J20 in the top right plot of Fig. 3). PG21 [44] solves the fast flavor pendulum, assuming homogeneity. Eq.(13) therein is an accurate description of the strictly homogeneous evolution, but it cannot be applied in more general inhomogeneous settings, e.g., to our Eq.(2). In fact, the homogeneous mode is typically stable for our ELNs. They also note that the truncated multipole approach is not accurate for homogeneous evolution. However, note that in our case inhomogeneity and dephasing are present. For all ELNs we have checked, the m f in Eq.(21c), derived by assuming dephasing and a truncated multipole tower, agrees well with the spatial , at 4096 different spatial locations at three different instants in time t = 0, t = 8 and t = 12. For this calculation, we used a non-random seed (see text). All panels show the data for v = −0.5 for G L v [vc = 0] and A = 0.8. One sees that at late times, here after t ≈ 12, the transverse Bloch vectors become large and random.

B. Transverse relaxation
The spatially averaged version of the flavor evolution given by Eq.(2) can be derived using approximations similar to those used in deriving Eq.(7) from Eq.(5), to get in a special corotating frame where we ensure that Hamiltonian for M 1 is purely transverse.  [40] for a derivation. In the remainder of this paper, corotating frame will refer to this special frame.
The length of each coarse-grained Bloch spin S v is predicted to remain constant according to Eq. (24). This isn't borne out by numerical calculations. The reason is simply that the spatially averaged equations are approximate in the first place. To understand this analytically, one needs to study the pre-coarse-grained partial differential equation in Eq.(5). Here we draw an analogy to the nuclear magnetic resonance of macroscopic samples, to obtain a semi-quantitative understanding.
In Eq.(24), the Bloch vector S v can be interpreted as the net spin of a macroscopic sample volume being acted on by a magnetic field equivalent to the corotating Hamiltonian h v . In reality, the macroscopic spin is composed of several microscopic spins bunched together, similar to how we have defined the coarse-grained S v from the precoarse-grained S v [z]. Initially, h v is along theê (3) direction and the Bloch spins remain aligned with h v . However, as M 1 tips over, h v decreases and concomitantly h ⊥ v increases. As a result, for some velocity modes, the transverse component of h v can become of similar size as its parallel component. We remind, the spatially averaged lengths of the parallel and transverse components  of h v , i.e., h v and h ⊥ v , respectively, are defined as follows: In Eq.(25b), we have used the fact that the length of the transverse vector remains invariant under the rotation aboutê 3 . The Bloch spins for those velocity modes develop a large precession angle to reach the transverse plane. At this juncture, the dispersion of the magnetic field h v [z], within the coarse-graining volume, can lead to the constituent microscopic spins to precess at different rates at different locations within the coarse-graining volume. This causes the transverse component of the macroscopic spins to become smaller over a timescale T2 -a process known as T2 relaxation. The transverse components of M 1 relax in the same way. As M 1 oscillates, for some velocity modes the relaxation turns on and off repeatedly.
The above analogy predicts that relaxation is strongest when the co-rotating Hamiltonian develops a large transverse component, i.e., when h ⊥ v ∼ h v . Roughly, this must coincide with M 1 developing a large transverse component. Further, one expects that transverse relaxation is prominent for those velocity modes for which h ⊥ v becomes comparable to h v . Conversely, for velocity modes whose transverse corotating Hamiltonians never grow too large, relaxation should be less efficient.
We will demonstrate the development of transverse relaxation using our numerical results for G L v [v c = 0], for two values of lepton asymmetry A = 0.8 (Figs. 4 and 5) and A = 0.2 (Figs. 6 and 7). The former case shows a slower rate and lesser degree of relaxation, while the latter shows faster and more extensive relaxation. For these four plots, we use a localized non-random seed: , i.e., the initial transverse components are taken to be ≈ 10 −6 , with fixed relative phase, localized around the centre of the box. This choice of seed (i.e., not random, unlike elsewhere in this paper) is to emphasize that even if we do not put random relative phases by hand, the system generates effective random phases on its own. The long-term results will be at best mildly sensitive to the seeds.
We begin with G L v [v c = 0] with A = 0.8. In Fig. 4, at t = 0 the transverse components of the Bloch vectors at all locations start out in phase (as set by the initial seeds in this case). By t ≈ 8 they begin to get dephased relative to each other, though the transverse vector is still very small at most locations. By t ≈ 12, this dephasing is essentially complete. In other words, the transverse components of S ⊥ v [z] become large and randomized across different locations z, as was shown in the bottom panel results of Fig. 6 in Ref. [40]. One thing to note is that S ⊥ v [z] = 0 without coarse-graining in z, but vanishes upon coarse-graining in z. Obviously, the transverse components of the coarse-grained multipole moments of S ⊥ v also decay due to this relaxation. In Fig. 5 the left panel shows the growth of the transverse components for v = −0.5. One sees that up to t ≈ 8 the evolution in linear. Yet, relative dephasing causes | S ⊥ v | to become smaller than |S ⊥ v | . Around t ≈ 10, close to the time of the first dip of the M 1 pendulum, nonlinearity sets in. The transverse components quickly grow to O(10 −2 ) and saturate. In the middle panel one sees that

intermittently or permanently, leading to the relaxation of corresponding modes (as seen in the right panel). For v > 0 one has h
For G L v [v c = 0] with A = 0.2 the relaxation is quicker, stronger, and more ubiquitous. In Fig. 6, one can see, the transverse components of the Bloch vectors at all locations start out in phase (as set by the initial seeds in this case). By t ≈ 4 they start to get dephased relative to each other, though the transverse vector is still very small at most locations. By t ≈ 6, this dephasing is essentially complete. In Fig. 7 the left panel shows the growth of the transverse components for v = −0.5. One sees that up to t ≈ 4 the evolution in linear. Yet, relative dephasing causes | S ⊥ v | to become smaller than |S ⊥ v | . Around t ≈ 4, close to the time of the first dip of the M 1 pendulum, nonlinearity sets in. The transverse components quickly grow to O(10 −2 ) and saturate. In the middle panel one sees that h v − h ⊥ v starts decreasing around t ≈ 4, owing to the growth of h ⊥ v . After t ≈ 8 one has h v − h ⊥ v < 0 for all the velocity modes, leading to the relaxation of corresponding modes (as seen in the right panel). In other words, one finds that the depolarization of S v [z] occurs if and when one has h ⊥ v h v . We have found this expected correlation for all the ELNs we have considered in this paper.
In Wu21 [45] doubts were raised whether the transverse relaxation mechanism holds in general, as they failed to find the correlation noted above. We investigated by repeating the computations of Wu21. Our results for the survival probability P ee , as shown in Fig. 8, are in quite  good agreement for the survival probabilities, showing partial depolarization. However, unlike Wu21, we clearly see the correlation expected from T2 relaxation. As shown in Fig. 9, for v < 0 modes one has h ⊥ v ∼ h v at around t ∼ 250, and around the same time but the difference is less compared to v = +1 mode and thus one finds almost no depolarization for v = +1 mode whereas partial depolarization for v = +0.5 mode. In Wu21 this comparsion was not made in the corotating frame and h ⊥ v was computed as the magnitude of the spatial average of the vectors h ⊥ v [z] (which is always close to zero due to dephasing), as opposed to average of the magnitudes, leading to their conflicting observation. Correcting for these misunderstandings, we find that the computations in B20a, B20b as well as in Wu21 are consistent with each other. The final depolarized state is almost entirely identical in both computations, and more importantly, the mechanism of T2 relaxation seems to work as we predicted when their example is analyzed as we recommended. There are minor differences because our code uses Fast Fourier Transform for differentiation in a way that creates ring-down effects around features that are sharp on the scale of the discretization scale 1 .
Smoother initial conditions do not get affected by this. Despite this difference, our predictions for the final survival probability agree to better than 5% r.m.s. error for the tested example.

C. Multipole cascade
The discussion in the previous subsection was limited to the first four multipole moments. In this subsection we review the nonlinear behavior of the higher multipole moments, as given in B20b [41]. Spatially averaging over Eq.(5) assuming periodic boundary conditions on z, and taking n 1, gives

Note in our approximation
where c 1 , c 2 are integration constants with Ei[x] = x −∞ dy e y /y. Eq.(26) and Eq. (27) together indeed indicate that there is a diffusion of the quantity M n [t] from low to high n multipoles as time passes causing irreversibility in the system. Due to such leakage of power from smaller moments, M n [t] for large n starting from some initial value grows exponentially to peak roughly around t peak n ≈ n 2 /(2 M 1 ) and then asymptotes to some steady final value at late times. Note t peak n increases with n. In B20b [41], we had shown this for the box-type ELNs denoted here by G B v . We have now verified it hold for all the ELNs considered in this paper. An example is shown in the top panel plot of Fig.10 with G L v for v c = 0 and A = 0.2. The takeaway is that the flavor difference increasingly gets moved to high multipoles. If a physical process does not distinguish closely spaced momentum modes, it no longer sees the flavor difference stored in high-multipoles.

D. Mixing of flavor waves
Flavor waves also cascade to smaller distance scales, similar to the cascading to smaller momentum scales we just discussed. This was shown very clearly in R21a [46] and R21b [47]. To understand this we take the Fourier transform of Eq. [x, t]dx to rewrite the following equation:  Fig. 10. We see that for t ≤ 2 the system is in the linear regime and the power, defined to be S ⊥ v [k, t], for a specific k mode does not cascade to other k modes. Until about t = 2, each curve grows with time exponentially for each k but with its characteristic k-dependent linear growth rate Im Ω[k]. In the linear regime one can clearly see that the footprint of instability is limited to the k-modes between ∼ (5 − 15) for our chosen example. By t ≈ 2.5 the modes close to k ≈ 8 have become large and they start affecting the growth of modes close to k ∼ 0, enhancing them considerably. This sudden distortion is a signature of mode-coupling in Eq. (28). Further, mode-coupling also allows the large-|k| modes with smaller amplitude to grow in a cascade, at the expense of the modes that start with higher amplitudes ,and thus spread the flavor instability to almost all k modes. This moves the flavor differences to smaller and smaller distance scales. If a physical process does not distinguish closely spaced locations, it does not see the flavor difference that is now stored in very high-|k| modes. Multipole diffusion and mode-coupling, together create extremely fine structures in the phase space, which upon coarse graining present themselves as effective depolarization.

E. Flavor waves vs. depolarization
In M21 [49], the authors speculated that the simulation tools used in our previous work in B20a [40] and B20b [41] may have failed to maintain causality. This speculation stemmed from the persistence of wavelike numerical solutions found in M20 [48] and M21 [49], as opposed to a depolarized state. In the mean time, other groups have found results that are broadly consistent with depolarization seen in our previous works (see e.g., [45][46][47]). Here, we reproduce the key results of M20 [48], to show our code produces results identical to theirs, if restricted to the regime they have explored. If extended to longer times, one finds depolarization.
To benchmark our code against the calculation in M20 [48], we focus on their G 3a and G 4b ELNs. Our dv as a function of z at various time snapshots up to t = 900 are shown in Fig. 11. The results agree, with excellent fidelity, with their counterparts in Fig.3 of M20 [48]. One clearly notices flavor waves in space and the region over which they exist extends with time as they propagate. Note the flavor waves show convective and absolute nature for G 3a and G 4b , respectively.
However, we believe that two important issues were ignored in in M20 [48] and M21 [49]. Firstly, the numerical results were shown only up to t = 900 when the system does not reach sufficient nonlinearity. Secondly, the quantities were not coarse-grained over a spatial volume. Both of these were important to obtain the irreversible steady state depolarized solution in our previous work. To clarify these two points, we scale up the neutrino ELNs G 3a and G 4b by a factor G 0 = 100 (i.e., instead of G 3a and G 4b , we consider the ELNs to be 100 × G 3a and 100 × G 4b , respectively) and otherwise retain exactly the same specifications, i.e., same box size, spatial discretization, initial condition, boundary condition, and so on, as in M20 [48]. In the top most curves for the top and middle panels of Fig. 11 we show our numerical results for S v=1 [z] (dark gray) and S ⊥ v=1 [z] (dark or light gray lines) as a function of z choosing G 0 = 100 and t N L = 12. One can clearly see that flavor waves break down at after reaching nonlinearity. Note that this respects L/2 t N L , required to avoid boundary effects due to the periodic boundary condition at late times. We show S v=1 (after spatial averaging) vs. t in Fig. 12, which shows that the system indeed reaches a flavor depolarized steady state.

IV. FLAVOR DEPOLARIZATION
To quantify the amount of flavor depolarization we define a depolarization factor in the following way: Note that t f is chosen to be large enough, such that the system has reached steady state. Full flavor depolarization leads to f D v = 0.5, whereas no depolarization is given by f D v = 0, and partial depolarization by f D v between 0 to 0.5. Sometimes one may find f D v > 0.5. This happens because the system first changes flavor almost completely, corresponding to a flavor conversion probability of 1, and then depolarized partially.
We show our numerical solution for f D v as a function of different velocity modes v in Fig.13 considering G B v , G L v and G C v for various choices of A > 0. Our numerical analysis suggests that depolarization is velocity-dependent: the negative velocity modes are almost always fully flavor depolarized for A > 0, but the positive ones are partially flavor depolarized. The extent of partial flavor depolarization depends on lepton asymmetry A, zero crossing position of neutrino angular distributions v c .

A. Extent of depolarization
In this subsection we analytically explain the functional dependence of f D v on A, v, v c and give an explicit linearized formula for f D v in terms of quantities determined from initial conditions. To derive this we use the numerical observation that S v [t f ] ≈ 0, i.e., f D v ≈ 1 2 , for v < 0, in all four cases based on our numerical analysis in Fig.13. This assumption, for A > 0, is motivated by our qualitative understanding of which modes get more depolarized. Using this and enforcing lepton number con- for v > 0 modes where we define the "forward" moments of the ELN as To obtain the linear order correction to the above result, we expand S v [t f ] as a function of v as where s 0 , s 1 are space-time independent constants but can depend on A and the nature of G v . Note s 0 , s 1 can be determined from the following formula: For our chosen form of G v , with A > 0 and a forward excess, we use S v>0 [t f ] ≈ A γ0 and S v<0 [t f ] ≈ 0 to deduce s 0 , s 1 from Eq.(32) as and Using Eqs.(31), (33), and (34) we can write f D v as In case of G B v , G C v we find γ 0 = 1 but for G L v , γ 0 = 1 − 2v c . Clearly the functional dependence in Eq. (35) indicates f D v for v > 0 modes decreases with increase in A and decrease in γ 0 (or in other words v c → 1). Plugging in the values for γ 0 , A and v c we get a good agreement between the numerical and analytical solution of f D v as a function v with v > 0 modes for all the cases except G C v as shown in Fig. 13. For cubic ELN, our linear approximations used in the above derivation might be inapppropriate since G C v itself contains only terms higher than linear order. Also, even for v < 0 modes, S v [t f ] to 0 is not entirely correct as we see. However, for reasonable values of asymmetry A ≈ 0.2, our prescriptions seems to work quite reasonably because the naive equilibration hypothesis with f D v = 0.5 for all modes is already a good approximation, and one only needs to "fix" the lepton number conservation constraint that is violated by naive equilibration. A small linear correction, as provided by our approach, provides such an improved estimate.

B. Three flavor generalization
Now that we have an estimate of the depolarization for two flavors, we seek its generalization to the real-world situation with three flavors. In general, this requires a completely new analysis [63,65]. However, if µ and τ flavors are taken to behave identically, the treatment is very simple. In such a case, the three flavor oscillations are treated in a restricted manner -with the ν e oscillating to ν µ and ν τ , democratically, and the oscillations between ν µ and ν τ being very efficient. Here, one can guess the effective three-flavor depolarization factor from symmetry considerations alone.
In Fig. 14, we show a section of the three-flavor Bloch volume -the so-calledê (3) -ê (8) triangle [73] -on which lie the states corresponding to pure flavor states. This region is an equilateral triangle with sides of unit length, with the vertices corresponding to flavor states. The twoflavor depolarization factor f D v is the distance of the tip from the top vertex along the left (or right edge). For three flavors, assuming µ − τ symmetry, the tip of the Bloch vector lies along the vertical perpendicular bisector. Note that transverse components of the Bloch vector (i.e., components out of the plane; in analogy to components orthogonal to an edge of the triangle for a twoflavor scenario) get T2-relaxed. The three-flavor depolarization factor f D, 3 flav v , to be used in Eq. (37), is then easily recognized as in terms of the two-flavor depolarization factor. Note that our analytical estimate of the two-flavor f D v , as in Eq. (35), stays between 0 to 1/2, which corresponds to f D, 3 flav v being in the range 0 to 2/3, as one would expect. Numerically, one finds the two-flavor f D v can sometimes exceed 1/2. This corresponds to predominant flavor conversion from ν e to say ν µ , and then partial depolarization. Heres, one expects a similar transition to the third flavor ν τ as well. The combined action projects the Bloch vector as shown by the lighter dashed grey lines. It is easy to see why: if ν e almost fully convert to ν µ (while ν µ and ν τ are symmetric), in a three-flavor framework ν e has zero survival probability, with equal conversion probability of 1/2 to both ν µ and ν τ .

V. PRESCRIPTIONS
The takeaway is that we expect depolarization to be the end-state of neutrinos that have undergone fast oscillations. Below, we provide two easily usable set of expressions related to fast oscillated neutrinos. Our intended users are supernova simulators in the first instance, and supernova neutrino phenomenologists for the second.

A. Sub-grid recipe for SN simulations
In supernova simulations, one computes the neutrino distribution function -whether in detail or using moments. See, e.g., Refs. [74][75][76][77][78][79]. The finite elements for these simulations are about 0.1 km in size, and it is inconceivable for the foreseeable future how one could faithfully include fast oscillations occurring on sub-cm scales into these already hugely expensive supernova hydrodynamic calculations.
Our proposal is that one should first identify each 'pixel' in the star where fast instabilities can exist. This can be accomplished using a variety of ways, including computationally efficient and increasingly more reliable approximations involving the moments of the neutrino distributions [80][81][82] or simply applying the crossing criterion [67,68]. Therein, to obtain an estimate of the effect of flavor oscillations, one should replace the original phase space distributions F ini with the depolarized distributions F depol : where x = µ/τ . As fast oscillations are insensitive to neutrino energy E, the same f D, 3 flav v applies to neutrinos and antineutrinos. Note that this does not impose naive equalization of all flavors, but a much less extreme mixing consistent with conservation laws. Of course, if perfect depolarization is allowed then the ν e distribution becomes 1 3 . This is easily recognized as the usual 1 : 1 : 1 mixture of the three flavors.
The main advantage of this sub-grid prescription is that one can avoid performing the expensive fast oscillation calculation, using an analytically pre-computed look-up table instead. Further, it implements a meaningful estimate of the oscillated distributions -conserving the relevant lepton asymmetry and carrying nontrivial momentum dependence of the degree of depolarization.

B. Depolarized flavor-dependent flux
To compute the terrestrially observable neutrino fluxes, we need the fluxes at at radius of say about 100 km from the center of the star, where fast oscillations have ceased and one has to then include slower collective effects, MSW transitions, etc. The procedure to include these slower effects are by now well understood. But suppose we only have the undepolarized primary fluxes provided by existing supernova simulations. How can we include an estimate of the depolarization? In general, this is complicated. However, making some symmetry assumptions, a simple estimate is possible.
We assume that the neutrino emission is axially symmetric at each point in the star and that the star is axially symmetric about the axis joining the star and Earth. Thus, the net observable flux from all source regions is simply given by summing over the velocity modes that leave that region in the direction parallel to the axis. The appropriately velocity-weighed depolarization factor is then given by Since we are considering fast oscillations, we approximate the putative neutrino-sphere as an infinite wall. As a result, only the v > 0 modes can be observed, with f D v < 1/2 always. Note that γ 0 is the zeroth forward moment of the ELN, cf. Eq. (30). Similarly we can define the n th forward moment of f D 3 flav v in the following way: Putting n = 0 in Eq.(39) gives back Eq. (38). The total flux per unit energy detected at a distance r from the neutrino sphere of radius R is where α = e, µ, τ . If we consider no oscillation then F For our calculation we assume the following v dependence of F Thus, knowing the distributions one can compute the coefficients a n ,ā n and b n , as well as the moments of the depolarization factor f n . Together, these allow one to compute the depolarized fluxes from the unoscillated fluxes 2 . For multidimensional simulations, one may have 2 Note that the differential ELN distribution is more detailed information that allows summing over the momenta appropriately, and the recipe in the previous section is superior in that case. However, if a neutrino phenomenologist wants to approximately readjust the primary fluxes predicted by a supernova simulation to include for potential effects of fast depolarization, the above recipe gives a crude but meaningful estimate.

VI. SUMMARY AND OUTLOOK
In this paper, we have presented detailed analytical as well as numerical analysis of the late time nonlinear behavior of a dense neutrino gas undergoing fast collective oscillations in the collisionless quantum kinetics approximation. Our study includes time-dependence, but is restricted to one spatial dimension and one nontrivial momentum coordinate that we have taken to be the radial velocity. Unbroken azimuthal symmetry around the radial coordinate is assumed. Under these assumptions, we find the following results: 1. The evolution of the average flavor content is similar to the motion of a pendulum. However, this pendulum neither preserves its length nor retains its periodic motion, as seen in Fig. 2. It settles down to a resting point, which is analytically known in terms of the ELN and its moments, cf. Eq.(21c), and shown in Fig. 3. misunderstandings: the former applied our criterion of comparing the Hamiltonian components in a non-standard way (see Fig. 9), and the latter didn't show results after spatial averaging at sufficiently late time (see Fig. 12).
5. The flavor content eventually acquires an approximately time-independent character. This is called depolarization. The extent of depolarization is nonuniform over neutrino and antineutrino momentum, as shown in Fig. 13. In general, it depends on the ELN. This is essentially because the net lepton asymmetry needs to remain conserved.
6. The extent of depolarization, encoded in the depolarization factor, can be predicted -if the range of fully depolarized modes is assumed. The prediction is based on a series expansion of the final flavor composition, and enforces lepton number conservation. Equation (35) gives an estimate to linear order in v, in the two-flavor approximation.
7. The above result is in the two-flavor approximation. Equation (36) generalizes it to a restricted three-flavor scenario where the initial conditions and evolution of the µ and τ flavors are taken to be identical.
8. The depolarized flavor distributions (in Eq. (37)) and the depolarized fluxes (in Eqs.(44a)-(44d)) are given in terms of the original distributions (in Eqs.(43a)-(43c)) and forward moments f n of the depolarization factor (in Eq. (39)). These are approximate but readily usable ingredients for implementation in supernova/nucleosynthesis simulations and for computations of neutrino signals.
Dephasing leads to qualitatively different results than purely coherent evolution. This the fundamental result we hope to have conveyed. Our treatment of depolarization rests on the idea that there is dephasing of many modes. It is the dephasing assumption that allows going from Eq.(2) to Eq. (7), allows truncation of the multipole equations, introduces irreversibility, leads to the steadystate solution in Eq.(21c), and allows a description of depolarization. While we do not use the truncated or dephased equations for any numerical computations, rather preferring to solve Eq.(2) directly and then averaging the solutions appropriately, the analytical results of the relaxed and truncated multipole equations, e.g., Eq.(21c), provide remarkable agreement with the numerical solutions of the full equations at late times.
The reader may see parallels with the "derivation" of the Boltzmann equation [83,84]. Hamilton's equations for many interacting particles can be expressed as the BBGKY hierarchy, but there is no way to truncate that hierarchy without assuming something more, viz., molecular chaos, coarse graining, etc. These assumptions serve to introduce, by hand, the loss of correlation required to explain irreversibility. While the derivation continues to be a matter of discussion, there is no doubt that its end result, i.e., the Boltzmann equation, is extraordinarily useful and describes macroscopic reality much more appropriately than the technically better justified microscopic equations of motion.
We conclude this paper with our outlook for further work on this subject. We believe that an immediate task is to arrive at a better estimate of the range of depolarized modes. Perhaps the answer will lie in devising an improved criterion on the Hamiltonian, or finding the exact depolarization envelop. With that, the problem of computing the depolarized final state of fast oscillating neutrinos would be largely accomplished. It is our belief that this will be important and useful for any practical study accounting for the fast flavor oscillations of neutrinos in supernovae.