Stationary solutions for fast flavor oscillations of a homogeneous dense neutrino gas

We present a method to find the stationary solutions for fast flavor oscillations of a homogeneous dense neutrino gas. These solutions correspond to collective rotation of all neutrino polarization vectors around a fixed axis in the flavor space on average, and are conveniently studied in the co-rotating frame. We show that these solutions can account for the numerical results of explicit evolution calculations, and that even with the simplest assumption of adiabatic evolution, they can provide the average survival probabilities to good approximation. We also discuss improvement of these solutions and their use as estimates of the effects of fast oscillations in astrophysical environments.


Introduction
Experiments with solar, atmospheric, reactor, and accelerator neutrinos have established that neutrinos produced in a specific flavor state oscillate among all three flavor states [1]. Neutrino oscillations studied by the above experiments depend on the vacuum neutrino mixing properties and forward neutrino-electron scattering in the relevant matter. Because flavor evolution of individual neutrinos can be treated separately in these cases, the theory is well understood and calculations are straightforward. In astrophysical environments such as supernovae and neutron star mergers, however, neutrino emission is so intense that forward scattering among neutrinos becomes important, which causes flavor evolution to be coupled for neutrinos emitted with different energies and directions. Consequently, collective flavor oscillations (see e.g., [2] for a review of the early studies) may occur for the dense neutrino gas in these environments. The treatment and understanding of such oscillations are still evolving and under intensive study (see e.g., [3][4][5][6][7][8] for some recent developments).
Many insights can be obtained by studying oscillations of astrophysical neutrinos in terms of mixing between ν e and ν x , which represents an appropriate linear combination of ν µ and ν τ . In this case, the flavor field can be conveniently described by the polarization vector P(ω, v) for a neutrino (antineutrino) emitted with energy E and velocity v. Here ω = ±δm 2 /(2E) is the vacuum oscillation frequency, δm 2 > 0 is the mass-squared difference between the two vacuum mass eigenstates, and ω > 0 (ω < 0) denotes neutrinos (antineutrinos). For clarity, the sans serif font and numeral indices are used for quantities in the flavor space, while the bold face and (x, y, z) indices are used for vectors in the Euclidean space. The probability of being Email address: z.xiong@gsi.de (Zewei Xiong) the electron flavor is (P 3 + 1)/2, where P 3 = P ·ê 3 andê 3 is the unit vector in the third direction of the flavor space. For coherent flavor evolution (i.e., in the absence of collisions) that starts with all neutrinos in pure flavor states, the P(ω, v) for the initial ν e and ν x (ν e andν x ) will differ only by an overall sign. Hereafter, P(ω, v) specifically refers to the polarization vector of the initial ν e orν e .
The general equation governing the spatial and temporal evolution of P(ω, v) is where H(ω, v) = H vac (ω)+H mat +H νν (v) is the total Hamiltonian. The term H vac (ω) = ωB accounts for the vacuum mixing, where B = (sin 2θ V , 0, − cos 2θ V ) for the normal mass hierarchy and B = (− sin 2θ V , 0, cos 2θ V ) for the inverted mass hierarchy with θ V being the vacuum mixing angle. The contribution H mat = λê 3 originates from forward neutrino-electron scattering, where λ = √ 2G F n e , G F is the Fermi coupling constant, and n e is the net electron number density. The contribution H νν (v) from forward scattering among neutrinos is our main concern and is discussed below. We is the fourvector corresponding to v with ρ running over (t, x, y, z), is the neutrino polarization current, µ = √ 2G F n ν e , n ν e = dv ∞ 0 dω F ν e (ω, v) is the ν e number density, F ν e (ω, v) with ω > 0 (ω < 0) is the ν e (ν e ) spectral and angular distribution function, and gives the factor (1 − v · v ) that is explicitly included in the usual expression of H νν (v). Here we emphasize the physical importance of J ρ and will discuss neutrino flavor evolution using its 12 components.
We focus on the so-called fast flavor oscillations [9,10]. When the angular distribution +∞ −∞ dω G(ω, v) of the electron lepton number (νELN) carried by a dense neutrino gas has a zero-crossing, i.e., switches from being positive to negative, an instability may be triggered, which could result in fast flavor conversion on length scales as short as ∼ O(1m) for conditions in neutron star mergers [11,12] and supernovae [13][14][15]. Whereas the above flavor instability can be identified by a linear stability analysis [5,6,16], the eventual outcome of fast oscillations is much harder to ascertain for realistic astrophysical environments. Apart from the uncertainties in modeling the neutrino emission in supernovae and neutron star mergers, numerical treatment of neutrino flavor evolution is further hampered by the complicated geometry and intrinsically dynamic nature of such environments. Consequently, studies of fast oscillations beyond the linear regime have been restricted to greatly simplified models so far. Specifically, Eq. (1) was solved for artificial νELN angular distributions by keeping the spatial derivative only in one [17? , 18] or two directions [19]. Other studies dropped either the time [20] or spatial derivative [21? , 22].
In this letter, we assume a homogeneous neutrino gas for which the spatial derivative can be ignored. We find the stationary solutions and compare them to the results from evolution calculations. We also discuss improvement of these solutions and their use as estimates of the effects of fast oscillations in astrophysical environments.

Stationary solutions
Under our assumption of homogeneity, Eq. (1) becomes The evolution of P(ω, v) for the initial ν e (ω > 0) differs from that for the initialν e (ω < 0) only through the vacuum term H vac (ω) in H(ω, v). For the dense astrophysical environments of interest to us, H vac (ω) can be ignored because the magnitude of ω is far less than that of λ or µ associated with H mat or H νν (v), respectively. The effect of H vac (ω) is to initiate neutrino flavor mixing, which can be approximated by allowing P(ω, v) to have small initial deviations from the pure flavor states. With this prescription, the evolution of P(ω, v) no longer depends on ω. We write P(ω, v) = P(v) and solve the evolution equation where is the angular νELN distribution. From Eq. (4), we obtain ∂ t J t = λê 3 × J t . For the small initial deviations of P(v) from the pure flavor states, J t is essentially parallel toê 3 initially (at t = t 0 ). Therefore, ∂ t J t ≈ 0 and J t ≈ J t t 0 ≈ µ effê3 , where µ eff = µ dv g(v) is specified by the νELN.
Based on the above discussion, we take The components of the polarization current J x , J y , and J z are vectors in the flavor space. At any specific time, the range of H(v) is determined by these vectors and the constraint (v x ) 2 + (v y ) 2 +(v z ) 2 = 1, and each P(v) precesses with the corresponding angular velocity H(v). We seek stationary solutions for which all polarization vectors collectively precess with the same angular velocity Ω on average. In the frame that co-rotates with these vectors, the evolution equation becomes We assume adiabatic evolution, for which the angle between P(v) and H (v) = H(v) − Ω stays fixed (see similar approach in [23,24]). Because of the constraint J t ≈ µ effê3 , Ω is approximately parallel toê 3 for the stationary solutions. For the small initial deviations of P(v) from the pure flavor states, J x , J y , and J z are also approximately parallel toê 3 at time t 0 . So P(v) is parallel to H (v) initially, and the adiabatic condition can be written as whereĤ (v) is the unit vector in the direction of H (v), and The angular velocity Ω and the polarization current components J x , J y , and J z for the stationary solutions can be solved iteratively by using Eqs. (7) and (8) in the definition of J ρ and applying the constraint J t ≈ µ effê3 . Clearly, it is most convenient to carry out the above procedure in the co-rotating frame.

Example solutions
We now illustrate the stationary solutions with specific examples. Because the angular velocity Ω is approximately parallel to the matter term λê 3 in H, λ effectively shifts the magnitude of Ω. We take λ = 0 and µ = 10 4 /(4π) km −1 . Motivated by the conditions in supernovae [25], we assume azimuthally sym- The range of γ = 1-4 allows the above g(v z ) to have any zerocrossing between v z = 1 and −1 (see Fig. 1). We focus on the cases of γ = 2 and 3 with a crossing at v z ≈ 0.38 and 0.07, respectively. We first solve Eq. (4) numerically using 600 bins for v z and 320 bins for the azimuthal angle of v. We have checked that convergence is achieved for this angular resolution. All polarization vectors P(v) are assigned random deviations between −10 −3 and 10 −3 from P 3 (v) = 1 at t = t 0 . Their evolution is followed up to t − t 0 = 2 km, when an approximately stationary state has been reached (see Appendix A). The results on the survival probability (P 3 + 1)/2 for γ = 2 and 3 are shown in Figs. 2(b) and 3(b), respectively. Because the azimuthal symmetry is broken by the initial conditions adopted to approximate the effects of the vacuum term in H(ω, v), these results depend  on the azimuthal angle of v for each bin of v z . The corresponding values of Ω, J x , J y , and J z are given in Table 1.
We next calculate the stationary solutions as presented in Section 2. For all practical purposes, we can assume pure fla-   vor states to obtain J t = µ effê3 , J x = J y = 0, and J z ê 3 at t = t 0 , which gives H t 0 (v) ê 3 . Specifically, the initial conditions are J t 3 = 2812 km −1 , J z 3 = −3030 km −1 for γ = 2 and J t 3 = −1028 km −1 , J z 3 = −3286 km −1 for γ = 3. Conservation of J t requires Ω = Ωê 3 . Due to the azimuthal symmetry of g(v z ) around the z axis in the Euclidean space, we can specify the x and y axes by setting J y 3 = 0 for the stationary solutions. In addition, the rotational symmetry aroundê 3 in the flavor space allows us to specify the directionsê 1 andê 2 of the corotating frame by setting J x 2 = 0. The remaining components of the neutrino polarization current J x 1 , J y 1 , J z 1 , J y 2 , J z 2 , J x 3 , and J z 3 for the stationary solutions can be solved from the iterative procedure in Section 2.
We find that the stationary solutions can be classified into the following types, with the subscript ⊥ denoting eitherê 1 or e 2 : (I) J x,y,z ⊥ ∼ 0, which represents the initial configuration, and therefore, is trivial, (II) J . While the above classification reflects specific relations among the neutrino polarization current components, such relations are not used to find the stationary solutions but are observed after the solutions are obtained. Note that large magnitudes of J x,y,z ⊥ correspond to substantial overall flavor transformation with large magnitudes of P ⊥ for wide ranges of v.
Because the procedure outlined in Section 2 involves solving integral equations, the results are not unique in general and are only candidate stationary solutions. We perturb the polarization vectors P(v) of a candidate solution with random deviations of magnitude up to 10 −3 and evolve them with Eq. (4) until the stability of the solution can be established or otherwise. The stable ones are selected as true stationary solutions. Note that the results may not be unique even after the above stability test.
For γ = 2, we find seven candidate stationary solutions, two of which are stable (types IIIa and IVa). In contrast, only one of the five candidate solutions are stable (type IIIa) for γ = 3 (see Appendix A). The survival probabilities (P 3 +1)/2 corresponding to the true stationary solutions are shown in Figs. 2 and 3 for comparison with the numerical results. The corresponding values of Ω, J x , J y , and J z are given in Table 1. For both γ = 2 and 3, the stationary solution of type IIIa has J x 3 = J y 3 = J z ⊥ = 0 and J x ⊥ ⊥ J y ⊥ with J x ⊥ = J y ⊥ . Consequently,Ĥ 3 and hence P 3 are independent of the azimuthal angle [see Eqs. (5) and (7)] with the corresponding survival probabilities exhibiting azimuthal symmetry. In contrast, for γ = 2, the stationary solution of type IVa has J x 3 0 and J x,y,z ⊥ 0. Therefore, the corresponding survival probabilities depend on the azimuthal angle.
The νELN distribution g(v z ) crosses zero at v z ≈ 0.38 and 0.07 for γ = 2 and 3, respectively. It can be seen from Figs. 2 and 3 that except for v z near the zero-crossings, the survival probabilities given by the stationary solutions show the same trends as the numerical results. For more detailed comparisons, we calculate the average survival probability for each bin of v z by averaging the numerical results over the azimuthal angle. Figure 3(a) shows that the stationary solution for γ = 3 describes the average survival probabilities very well away from the zero-crossing. Because there are two types of stationary solutions for γ = 2, we calculate the corresponding average survival probability for each bin of v z by weighing each type of solution equally and then averaging the results over the azimuthal angle. Figure 2(a) shows that these average survival probabilities are also in good agreement with the numerical results away from the zero-crossing.
The large deviations of the survival probabilities for the stationary solutions from the numerical results for v z near the zerocrossings are caused by the breakdown of adiabatic evolution assumed in deriving these solutions. The corresponding neutrinos have small values of |H t 0 (v)| = |H 3,t 0 (v)| µ that result in nonadiabatic flavor evolution. The effect of this nonadiabaticity can be approximated by including a component of P(v) that is perpendicular to and rapidly rotating around H (v) in addition to the component aligned with the latter and given by In the above equation,Ĥ (v) is taken from the stationary solution and eff (v) is taken as the error function is the mean value due to the component aligned with H (v). The corresponding survival probabilities are shown as the blue, green (limits), and red (mean) symbols in Figs. 2(d), 2(f) for γ = 2 and in Fig. 3(c) for γ = 3. It can be seen that these results describe both the trends and the range of the numerical results rather well. Because the rapidly rotating component of P(v) perpendicular to H (v) is essentially averaged out, only P align (v) is effectively used to find the stationary solutions. We can choose eff (v) and use Eq. (10) to replace Eq. (8) in the procedure to find these solutions. As an example, we choose eff (v) = erf 1.5H 3,t 0 (v)/µ (14) for γ = 3 and obtain a new stationary solution, which is also of type IIIa (denoted as IIIa ). The corresponding survival probabilities (mean and limits) are calculated in the same way as for Fig. 3(c) and shown in Fig. 3(d). The corresponding values of Ω, J x , J y , and J z are given in Table 1. It can be seen that the new stationary solution is much closer to the numerical results. While the effect of nonadiabatic evolution can be well approximated by choosing an appropriate eff (v) as in the above example, we note, however, that our original procedure to find the stationary solutions assuming adiabatic evolution is more straightforward and can already produce average survival probabilities to good approximation. As discussed above, nonadiabatic evolution appears to be associated with the zero-crossing of the νELN distribution g(v z ).
On the other hand, if g(v z ) has no zero-crossing as for γ = 1 or 4 (see Fig. 1), the initial values of P 3 (v) = 1 make J t 3 reach the most positive or negative value for γ = 1 or 4, respectively. Conservation of J t 3 then ensures P 3 (v) = 1 subsequently, and therefore, there is no flavor evolution at all [5].

Discussion and conclusions
We have presented a method to find the stationary solutions for fast flavor oscillations of a homogeneous dense neutrino gas whose angular νELN distribution has a zero-crossing. These solutions correspond to collective precession of all neutrino polarization vectors around a fixed axis in the flavor space on average, and are conveniently studied in the co-rotating frame. We have shown that these solutions can account for the numerical results of explicit evolution calculations, and that even with the simplest assumption of adiabatic evolution, they can provide the average survival probabilities to good approximation. These solutions can be further improved by including the effect of nonadiabatic evolution, which in turn, can be approximated by choosing an appropriate ansatz for the alignment of the polarization vectors with the effective Hamiltonian in the co-rotating frame. While we have focused on specific νELN distributions for our examples, we show in Appendix B that our method applies to other νELN distributions as well.
Besides shedding light on the nonlinear regime of fast oscillations, the stationary solutions discussed here provide physically motivated estimates of the average survival probabilities beyond the simple limit of complete flavor equilibration. Because these solutions can be efficiently calculated, they may be incorporated in simulations of dynamical astrophysical environments such as supernovae and neutron star mergers, for which the computational resources must be almost exclusively devoted to hydrodynamics and regular neutrino transport by necessity. We note that our method assumes a homogeneous neutrino gas but the environments where fast oscillations may occur are most likely inhomogeneous. Consequently, the stationary solutions discussed here may only provide a crude yet efficient way to explore the effects of fast oscillations on neutrino-related processes in supernovae and neutron star mergers. The general problem of flavor evolution of an inhomogeneous dense neutrino gas is beyond our scope here and requires further dedicated studies.
Appendix A. Numerical tests Figure A.4 shows the snapshots of our explicit evolution calculations at t = 1.6, 1.8, and 2 km for the assumed νELN distributions with γ = 2 and 3. It can be seen that an approximately stationary state has been reached by t = 2 km in both cases. Figure A.5 shows the evolution of J z 3 for some candidate stationary solutions after they are perturbed. Note that the results for type I are calculated from the initial configurations and become approximately constant at late times as expected of the approximately stationary states shown in Fig. A.4. For γ = 2, only types IIIa and IVa are the true stationary solutions. Type IVb is quasi-stable for t < 0.3 km but deviates from the initial state subsequently. Therefore, it is not a true stationary solution. Note that the average value of J z 3 for types IIIa and IVa is close to the late-time value of J z 3 for type I. For γ = 3, only type IIIa is the true stationary solution.

Appendix B. Other νELN distributions
We have repeated the calculations for two more νELN distributions taken from [? ] (case A) and [6] (case B): The zero-crossing occurs at v z ≈ 0.38 and 0.6 for g A (v z ) and g B (v z ), respectively. Figure B.6 shows that an approximately stationary state has been reached by t = 2 km for the explicit evolution calculations for both cases A and B, and compares the survival probabilities from the stationary solutions with the numerical results. For case A, we find three nontrivial candidate stationary solutions of types II, IIIb, and IVb. Only type IVb is stable (see  The flavor evolution for the νELN distribution of case B differs from all the other cases considered in that it is associated with a sign flip of J z 3 . As J z 3 changes from the initial value of −803 km −1 to ≈ 450 km −1 for the approximately stationary state (see Fig. B.7), neutrinos with a wide range of v z undergo nonadiabatic evolution. This effect cannot be captured by our approximate treatment of nonadiabatic evolution, which is appropriate only for a limited range of v z around the zero-crossing.
We note, however, that for the likely inhomogeneous astrophysical environments, the dominant instability associated with the νELN distribution of case B grows differently on different spatial scales [6,17]. Therefore, a realistic treatment of this case must address both the temporal and spatial flavor evolution of neutrinos, which is beyond our scope here. Symbols at the same v z are for different azimuthal angles. In panels (a), (b), and (c), snapshots of the numerical results are shown for t = 1.6, 1.8, and 2 km, respectively, for case A. The magenta curves are the survival probabilities averaged over the azimuthal angle. In panel (d), the symbols show the stationary solution for case A, the orange curve is obtained by averaging these symbols over the azimuthal angle, and the magenta curve is the same as that in panel (c). In panel (e), the symbols show the survival probabilities (red: mean, blue and green: limits) including the approximate effect of nonadiabatic evolution for the stationary solution for case A. Panels (f)-(j) are similar to panels (a)-(e), but for case B.