A few moments to diagnose fast flavor conversions of supernova neutrinos

Neutrinos emitted from a supernova may undergo flavor conversions almost immediately above the core, with possible consequences for supernova dynamics and nucleosynthesis. However, the precise conditions for such fast conversions can be difficult to compute and require knowledge of the full angular distribution of the flavor-dependent neutrino fluxes, that is not available in typical supernova simulations. In this paper, we show that the overall flavor evolution is qualitatively similar to the growth of a so-called `zero mode', determined by the background matter and neutrino densities, which can be reliably predicted using only the second angular moments of the electron lepton number distribution, i.e., the difference in the angular distributions of $\nu_e$ and $\bar{\nu}_e$ fluxes. We propose that this zero mode, which neither requires computing the full Green's function nor a detailed knowledge of the angular distributions, may be useful for a preliminary diagnosis of possible fast flavor conversions in supernova simulations with modestly resolved angular distributions


I. INTRODUCTION
The interior of a supernova (SN) hosts a unique laboratory to probe quantum correlations between neutrinos. For instance, at distances r O(10 2 ) km from the centre of the SN, the neutrino density n ν is so high that they themselves produce a collective potential µ = √ 2G F n νe , defined in terms of the electron neutrino density n νe . This potential, being much larger than the neutrino oscillation frequency in vacuum, ω vac = ∆m 2 /(2E) for a typical neutrino energy E, leads to correlated neutrino flavor evolution [1][2][3]. The past decade of research in this area has unearthed many fascinating features in the collective oscillations of neutrinos [4][5][6][7][8][9][10]. See refs. [11][12][13][14] for recent reviews.
In a series of papers [4,8,10], Ray Sawyer has pointed out a new mechanism for self-induced flavor conversions called "fast" instabilities. These are expected to develop at very short distances, r O(1) m, from the neutrinosphere and grow with a rate µ, i.e., not only faster than the usual neutrino oscillations but also than the relatively slower collective oscillations, that lead to spectral splits and swaps [7,9,15,16], growing at a rate √ ω vac µ [6].
Recently, several groups have confirmed these results and further developed the original insights [17][18][19][20][21][22][23]. In particular, it has been understood that a necessary condition for fast conversions is that there is a "crossing" in the electron lepton number (ELN) angular distribution, i.e., the difference of the ν e andν e densities must change its sign as a function of emission angle. This crossing condition is similar to how collective spectral swaps require a crossing in the energy spectrum [9]. In the neutrino decoupling region inside a SN, where the different flavors have significantly different angular distributions, a crossing in the angular spectrum could be present. As a result, fast conversions may occur and lead to potentially radical changes in SN dynamics and neutrino signals.
The possibility of fast conversions need to be explored systematically in SN simulations. First steps in this direction were taken recently [24], where a dedicated analysis of the angular distributions of the neutrino radiation field for several spherically symmetric (1D) supernova simulations has not found any crossing in the ELN near the neutrinosphere. Ref. [23], on the other hand, found an instability in a 8.8 M electron capture SN simulation by the Garching group. More generally, 2D or 3D models can exhibit Lepton-Emission Self-sustained Asymmetry (LESA) [25], i.e., a large-scale dipole in the ELN emission, which also makes a crossing more likely to occur. Unfortunately, a study of fast oscillations in 2D or 3D simulations has been lacking for two reasons. Firstly, the study of fast neutrino instabilities requires characterizing the singularities of the Green's function of the system [20]. This is a computationally demanding task even for the simplest toy models, and perhaps prohibitively difficult for the multidimensional continuous angular distributions found in SN simulations. Secondly, most of the state-of-the-art simulations [25][26][27][28][29][30][31][32][33], maintain only the moments of the angular distributions of fluxes, and not the full distributions. This lack of information seems to preclude even a linear stability analysis that requires knowing these distributions. One may in fact worry, whether these coarse-grained distributions can correctly capture the physics of fast oscillations. Therefore, it is necessary to consider an alternative approach that uses available simulations in an optimal fashion.
In this work, we propose a simple analytical tool to diagnose fast instabilities. Our proposal is based on identifying a specific Fourier mode of the flavor instability field, that we call the "zero mode". The growth rate of this mode, calculated from the stability analysis, crudely approximates the growth of flavor conversions in detailed numerical calculations, essentially all the way until the instability saturates. The zero mode has an easily calcu-lable growth rate, that depends only on the second moments of the ELN. Thus, the proposed method has the dual advantage of not requiring complete knowledge of the neutrino distributions and being computationally far less expensive compared to a full-fledged numerical solution or a full characterization of the Green's function. In fact, in the absence of more detailed knowledge of the ELN distributions, as appears to be the case for available 2D and 3D SN simulations, this may be the only practical recourse to look for possible instabilities. Therefore, we expect that this method will be useful to scan the different regions of a SN in multidimensional simulations and study the possibility of fast flavor conversions therein.
We discuss these issues in the following sections. In Sec. II, we write down the equations of motion (EoMs) and review the framework for studying fast neutrino oscillations. In Sec. III we present our method for diagnosing instabilities and in Sec. IV perform numerical tests of the same, for simple box-like angular distributions for the ν e andν e as well as realistic angular distributions inspired by 1D SN models. Finally, in Sec. V, we conclude with a brief summary.

II. EQUATIONS OF MOTION
Neglecting collisions, the dynamics of p , the matrices of neutrino phase space occupation number densities for the momentum p, is described by the following EoMs [34][35][36][37][38][39] where, in the Liouville operator on the left-hand side, the first term accounts for explicit dependence on time t, while the second term, proportional to the neutrino velocity v p , encodes the dependence on position x due to particle free streaming. The right-hand-side contains the oscillation Hamiltonian where in a two-flavor approximation, one has in the mass basis, and in the weak interaction basis, contains the refractive effect of charged leptons in the medium, while contains the similar effect due to background neutrinos. Antineutrinos are described similarly using¯ p , with H vac replaced by −H vac .
The matrix can be written in the weak-interaction basis as [19,20,40,41] = f νe + f νx 2 where f νe and f νx are the occupation number densities at momentum p, and ν x is the relevant linear combination of ν µ and ν τ . Here and onwards, we drop the subscript p, which indicated that the were indexed by their momenta, to lighten the notation. One focuses on length and time scales over which f νe and f νx can be taken to be homogeneous and static, and thus the spatial and temporal dependence of is contained in S and s. The complex scalar field S(t, x) encodes the ν e ν x mean-field flavor coherence, and measures the extent of flavor conversion. To begin with, the neutrinos are initially in their unoscillated states, hence the initial condition is S(0, x) = 0. The real field s(t, x) encodes flavor occupation number, and satisfies s 2 + |S| 2 = 1, for each momentum p.
In the context of fast conversions, the effect of background neutrinos via H νν far exceeds that of the vacuum Hamiltonian H vac , which mainly plays the role of generating an initial disturbance to seed the oscillations. Hence, H vac can be neglected and the explicit dependence on energy E, via ω vac , disappears from the EoMs. The neutrino and antineutrino modes then enter the Hamiltonian in Eq. (5) only via the difference of occupation number densities integrated over energy, i.e., the ELN angular distribution [17], where we assume ν x andν x have identical distributions. We will often use the "4-vector" notation, e.g., a µ = (a 0 , a), advocated in Ref. [19]. For the familiar quantities, i.e., position x µ = (t, x), momentum p µ = (E, p), and wavevectors k µ = (ω, k) and K µ = (Ω, K), the zeroth component is denoted by its more recognizable symbol instead. The neutrinos are taken to be ultra-relativistic, with E = |p|, so v µ = (1, p/E), i.e., the zeroth component of their velocity is 1 and the spatial components are given by a unit vector v = p/E. In this notation, one can define a matter current which contain the effect of H mat and H νν , respectively.
The key feature which determines if the initial flavor composition is unstable and can undergo fast conversions, is if the ELN distribution G v crosses zero as a function of any of its arguments. This essentially requires that the flux of neutrinos is larger than that of antineutrinos in some direction, while being smaller in another direction. If the lepton asymmetry, ε = (n νe − nν e )/n νe = Φ 0 /µ is small, then the ELN distribution G v can have a crossing. This is because the density of forward-goingν e can exceed that of ν e , between theν e and ν e neutrinospheres.
The onset of the conversions can be examined by linearizing Eq. (1), using that initially |S| 1 and s 1. Starting from the linearized equations of motion [40], one seeks plane wave solutions obeying [42] S A specific eigenmode of flavor conversion can be labeled by its frequency and wavevector, Ω and K, respectively. If there are modes that have a complex Ω, such modes may lead to exponentially growing instabilities. The currents Φ µ and Λ µ lead to a common rotation for all modes, which does not lead to any instabilities as such. So it is more convenient to work in a rotating coordinate system where this common rotation is not present. In such a co-rotating frame, the different modes of flavor conversions are labeled by the shifted frequency and wavevectors, respectively. Note that the frequency and wavevector of the modes in the co-rotating frame, ω and k, have the same imaginary parts as in the non-rotating frame. Hence this shift, while simplifying the EoMs, does not give rise to any extra spurious instabilities. Of course, one must be careful of the shift when identifying a specific mode of the flavor conversion field S, e.g., the homogeneous mode, previously labeled by K = 0, now corresponds to the mode k = −(Λ + Φ).
In the remainder of the paper we will refer to the k = 0 mode as the "zero mode". This mode will be the focus of our work, and we will come back to it in the next section. We end this section with a brief remark about the role of ordinary matter density. From the definition ω = Ω−(Λ 0 +Φ 0 ), one sees that Im(ω) has no dependence on the ordinary matter density encoded in Λ 0 , which merely leads to a shift in Re(ω), as noted in Refs. [42,43]. Presence of a finite matter density only delays the onset by suppressing the mixing angle, keeping the growth rate same. A nonzero ordinary matter current Λ on the other hand can change Im(ω), but is expected to be negligible in a SN-like environment where the ordinary matter has small velocity anisotropy. In the remainder of this paper, we will ignore effects of ordinary matter.

III. ZERO MODE AND MOMENTS
Our proposal in this paper is to focus on the zero mode, labeled by k = 0. This is motivated by the fact that the calculation of ω for this mode is significantly simpler than a full characterization of the roots of the dispersion relations, D(ω, k) [20]. In fact, for this mode the ω in Eq. (12) can be pulled out of the integrals, and Eq. (11) becomes i.e., D(ω, 0) is a polynomial in ω. The specific model of SN neutrino fluxes and their angular distributions, encoded in the ELN, only enters the equation through the tensor V µν that contains the second moments of velocity, namely, This, in turn, depends on the second moments of neutrinos velocities evaluated using the flavor-dependent phase space distributions, i.e., where the notation . . . να refers to For spherically symmetric SN simulations, which are effectively 1D, Eq. (13) is explicitly quadratic in ω , with the solution The condition for the zero mode to become unstable is simply that the discriminant become negative, i.e., If this condition is satisfied, the mode grows at a rate For multidimensional models, i.e., for 2D or 3D simulations, Eq. (13) is a cubic or quartic equation for ω, respectively. In either case, the instability growth rate is similarly calculable as the imaginary part of ω, if the flavor-dependent second moments v µ v ν να are known. We propose that one should check for fast instabilities in a SN simulation, by testing the condition in Eq. (19) locally in each simulation cell. It is our understanding that SN simulations often do not track the complete distribution f (E, v). This precludes an exhaustive search for fast instabilities by studying the solution of the dispersion relation. However, one can learn about the stability of the zero mode without such detailed information. If the ELN distribution G v is spherically symmetric, the tensor V µν has no cross terms and only depends on 1 να , v r να , and v 2 r να , i.e., the zeroth, first, and second moments of the radial velocity. Such information is readily available even in the spherically symmetric SN simulations and a search for instabilities using Eq. (19) is straightforward. In general, the terms like v µ v ν να are important. Some multidimensional SN simulations can provide these cross moments and may allow one to search for fast instabilities. In these cases, if Eq. (13) has complex solutions for ω in some region in a SN simulation, it indicates that fast conversions should occur.
Finally, we note that, besides the zero mode we have identified, there are two other important modes. The true homogeneous mode of the system is given by K = 0, which is now labeled by k = −(Λ + Φ). This mode, that has conventionally been studied in calculations that enforce an evolution in time and ignore spatial variations, need not have an instability as will be clear from some of the examples we study in the next section. Nonetheless, our method cannot be used to predict the behavior of this mode. On the other hand, the mode with the maximum growth rate can be determined using the condition that the group velocity of that flavor wave is zero, i.e., ∂ω/∂k| kmax = 0 [20]. Predicting this mode and its growth rate also requires knowing the full dispersion relation. Although our method is not useful to study these modes directly, we will find that the exponential growth of the zero mode, accurately predicted by Eq. (20), is a good proxy for the overall flavor evolution.

IV. NUMERICAL TESTS IN 1D
In this section, we demonstrate our proposed method using the flavor evolution of models with continuous ELN, in one spatial dimension z and time t, neglecting ordinary matter (i.e., Λ = 0). In one spatial dimension, the k and Φ vectors can be labeled by their magnitudes k and Φ, respectively. We numerically solve the nonlinear partial differential EoMs, i.e., Eq. (1), for several models and compare the flavor evolution, thus obtained, with the growth rate predicted by the corresponding moments. First, we will consider a few toy examples with box-like angular distributions for the ν e andν e , with a crossing at v = v c , and then show a calculation with more realistic distributions inspired by SN simulations.
We work in units such that the neutrino potential µ = 1, and times and lengths are expressed in units of µ −1 . In vacuum, the oscillation frequency is taken to be ω vac = 9 × 10 −5 while the effective mixing angle is ϑ = 10 −3 . We assume an inverted neutrino mass ordering, but the results are insensitive to this choice. The solution is found over the z-t plane, allowing z to take values in the interval (0 : z max ) and t in (0 : t max ). The boundary conditions are chosen such that flavor-pure modes are emitted at all z when t = 0, with S(z, 0) being a Gaussian wavepacket centered around z = 100 with small width σ = 1 and initial amplitude 10 −6 for both the real and imaginary part of S. Initially, the k = 0 mode peaks, and the seeds for all other k modes are smaller. As the first case, we take a box-like distribution given by (G νe , Gν e ) = (0.3, 0.5) for −1 < v < 0 and (G νe , Gν e ) = (1.2, 0.5) for 0 < v < 1. The ELN, G νe −Gν e , changes sign and presents a crossing at v c = 0. In this case, there are counter-going neutrinos, and one expects the instability to be absolute, spreading around its origin without drifting [20].
In the left panel of Fig. 1 we show the numerical evolution of |S(z, t)|, in the z-t plane. This is obtained by numerically solving the nonlinear partial differential EoMs for the model. An instability corresponds to a growth of |S(z, t)|. Here, we see that an absolute instability is generated at t 30 and gradually spreads over space without drifting. The nonlinear regime is reached at t 60, when |S(z, t)| ∼ O(10 −1 ).
In order to estimate the temporal growth of the instability, we look at the Fourier transform of S(z, t) as a function of t, The S(z, t) we obtain from the numerical solution of the EoMs is not in the co-rotating frame, and its Fourier modes correspond to the unshifted wavenumbers labeled by K, as denoted above. However, we will continue to work in the co-rotating frame and relabel the modes using k = K − Φ to obtain S k (t) = 1 z max zmax 0 dz e i (k+Φ)z S(z, t) .
We remind, the zero mode is labeled by k = 0, the homogeneous mode by k = −Φ, and the mode with maximum growth by k = k max .
In the right panel of Fig. 1 we plot log 10Ŝk (t) versus t as obtained from the numerical data for different k modes. The zero mode, labeled by k = 0 (continuous red line) shows an initial slow growth until t 30, followed by an exponential rise until t 60, after which the conversions saturate. The initial slow phase has been identified as an "onset" phase and depends logarithmically on the initial seed [22]. In the regime with exponential growth, the agreement between the numerical solution and the analytical prediction of the growth rate (dashed red line) is excellent, both showing a growth rate Im(ω) = 0.30. For comparison, we also show the true homogeneous mode, given by k = −Φ (unbroken green line) and the mode with the largest growth, given by k = k max (unbroken blue line). Note that in this example, the homogeneous mode does not have a linear instability. It is marked by a larger seed and a longer onset period, before non-linearity sets in. The k = k max mode, on the other hand, clearly has a larger growth rate than the zero mode and becomes non-linear earlier. However, there is no simple analytical expression for the growth rate of this mode. One needs to find the complete dispersion relation, setting k = k max in Eq. (13). We have identified this mode by numerically searching for the highest growth rate across all k.
As a second case we consider a box-like distribution with only forward-going neutrinos (v > 0). In particular, we take (G νe , Gν e ) = (0.6, 1.0) for v < 0.5 and (G νe , Gν e ) = (2.4, 1.0) for v > 0.5. In this case there is a crossing in ELN at v c = 0.5 but no counter-going neutrinos and we expect a convective instability, where the instability convects away from its point of origin [20].
The numerical solution of the EoMs for this case is shown in left panel of Fig. 2. One finds that the instability is indeed convective, drifting away from its origin at z ∼ 100 as it grows. In the right panel of Fig. 2, we compare the growth rate for the zero mode (continuous red line), obtained from the numerical solution of the EoMs, with the analytical growth rate Im(ω) predicted by Eq. (20) (dashed red line). Clearly, in the exponential growth regime, starting at t 40, the agreement is again excellent, with a growth rate Im(ω) = 0.17, all the way until saturation. The true homogeneous mode (unbroken green line), and the mode with the maximum growth (unbroken blue line) are also shown for comparison. As a final example, we consider realistic angular distributions, inspired by 1D SN models simulated by the Garching group, shown in the left panel of Fig. 3. As discussed in the introduction, most 1D SN simulations do not present a crossing in the ELN, unlike what may be expected in multidimensional SN models. However, the angular distributions are expected to be similar, and we only change the relative weights of ν e andν e fluxes within the range predicted by models exhibiting LESA to get a crossing in ELN at v c = 0.4, as shown in the right panel of Fig. 3. This ensures that the model shows fast conversions.
In the left panel of Fig. 4 we show the numerical solution of the EoMs for these realistic angular distributions. The dense neutrino cloud has counter-going neutrinos and the instability is absolute, spreading across space at t 50, without drifting away completely. For this model we found, using Eq. (20), the growth rate to be Im(ω) = 0.17. As shown in the right panel, the numerically computed growth rate in the true model for the zero mode (continuous red line) agrees well with the analytical prediction using the moments (dashed red line). The development of the true homogeneous mode (unbroken green line) and the mode with the maximum growth (unbroken blue line) are similar. In fact, in this particular example, the maximum growth rate is roughly the same as that of the zero mode.

V. DISCUSSION AND CONCLUSIONS
Fast flavor conversions can occur close to the neutrino decoupling region in a SN, where ELN angular distributions might have crossings. If these conversions take place, they would bring into question the current paradigm of SN simulations that do not include neutrino flavor conversions in the spectra formation and in SN dynamics [44]. Therefore, it is imperative to scan over a large sample of multidimensional SN simulations and search for fast neutrino flavor instabilities. Unfortunately, the relevant length scales for fast conversions are much smaller than the resolution of SN simulations and most multidimensional supernova simulations do not provide detailed neutrino angular distributions, but rather only their integrated moments. Also, numerically identifying the singularities of the Green's function, in order to find all possible instabilities, is difficult and time consuming. Therefore, it is perhaps necessary to adopt a schematic implementation of these effects. We argue that it is possible to analytically calculate the growth of the zero mode, determined by the background matter and neutrino densities, using only the first three moments of the angular distributions. While the zero mode does not necessarily have the largest growth rate, it allows an easy estimate of possible fast flavor conversions in a supernova. Using simple toy examples of box spectra, as well as realistic angular distributions inspired by SN simulations, we have demonstrated that the numerically computed growth rate for the zero mode exactly matches the predictions from the moments of the angular distributions. For completeness, we have also shown the homogeneous and fastest growing modes in all these examples. In the cases we have checked, the zero mode gives a good indication of the timescale over which fast instabilities lead to large flavor conversions.
We believe that this simplified approach to fast flavor conversions may allow rapid progress in this line of research. Indeed, with our simple recipe given in Sec. III, it should be possible to perform a preliminary scan for possible fast flavor instabilities in 2D and 3D supernova (and neutron star merger) models. SN simulators are likely to find the computational cost of this method to be significantly lower and might want to use it as a consistency check on their simulations. If unstable cases are found, this would have a profound impact on SN simulations and one would be forced to include the effect of fast conversions into state-of-the-art SN simulations in order to obtain a correct description of the SN dynamics and of the observable SN neutrino fluxes.