Fast Flavor Oscillations of Astrophysical Neutrinos with 1,\,2,\,\ldots,\,\infty Crossings
PPrepared for submission to JCAP
TIFR/TH/20-47
Fast Flavor Oscillations ofAstrophysical Neutrinos with1, 2, . . . , ∞ Crossings
Soumya Bhattacharyya ID and Basudeb Dasgupta ID Tata Institute of Fundamental Research,Homi Bhabha Road, Mumbai 400005, IndiaE-mail: [email protected], [email protected]
Abstract.
In the early Universe, as well as in supernovae and merging neutron stars,neutrinos have such high densities that they affect each other and exhibit collectiveflavor oscillations. A crucial ingredient for fast collective flavor oscillations is that theelectron lepton number (ELN) distribution changes its sign as a function of direction,i.e., has a zero crossing. We present a study in two dimensions and show how fastflavor oscillations depend on the ELN and its crossings. We show that a large numberof crossings can inhibit flavor oscillations. This may be a natural self-limiting mecha-nism that stabilizes the flavor content of the dense neutrino gas in a vast majority ofscenarios, especially the early Universe, where the angular distributions for all flavorsare very similar and crossings occur mainly due to fluctuations. a r X i v : . [ h e p - ph ] J a n ontents ∞ Crossings 52.2.1 True vs. Apparent Dimensionality 62.2.2 ELN Models 72.3 Linear Stability Analysis 82.4 EoM Solver 102.5 Dispersion Relation Solver 11
Astrophysical neutrinos are a valuable probe of fundamental physics and astrophysics.The focus of this paper is the so-called fast collective flavor evolution of neutrinos,which can have a dramatic effect in dense astrophysical environments, e.g., on theexplosion of core-collapse supernovae, nucleosynthesis in supernovae and neutron starmergers, and perhaps for neutrinos in cosmology.Neutrino oscillation in ordinary matter is governed through two frequency scales –one is the vacuum oscillation frequency, ω E = | ∆ m | / (2 E ), related to the mass-squareddifference between two different mass eigenstates, ∆ m = m − m , and the other is thematter potential, λ = √ G F n e , coming from the coherent forward scatterings of neutri-nos with electrons with possibly spatially varying density n e in the medium. Broadly,high density, λ (cid:29) ω E , suppresses flavor mixing. As a result, large flavor conversionin ordinary matter can occur either after λ drops below ω E , whence ordinary neutrinooscillations ensue, or when these two frequency scales match, i.e., λ ≈ ω E , giving riseto the well known phenomenon of matter-enhanced oscillations [1, 2]. Thus one mightthink that neutrinos cannot oscillate deep inside stars or in the early Universe. The quantum mechanical amplitude of forward scattering interferes with free propagation andgives a potential for flavor oscillations. Momentum and number changing collisions do not produceinterference effects. – 1 –eutrino-neutrino scattering in a neutrino gas gives rise to a new scale, µ = √ G F n ν , proportional to the neutrino density n ν [3]. This new potential makes theflavor evolution nonlinear, allowing novel collective flavor oscillations in dense environ-ments. A salient feature of these collective oscillations is that neutrinos of differentenergies oscillate approximately at the same rate, i.e., the average ω E for synchronizedoscillations [4], √ ω E µ for bipolar or slow collective oscillations [5, 6], and as rapidlyas µ for fast collective oscillations [7]. The slow variant of these collective oscillationswas studied deeply starting from the mid-2000s, exposing a sequence of new effectssuch as swapping of the flavor-dependent spectra as well as momentum, space, andtime-dependent flavor transformations [8–19].Fast flavor oscillations became a topic of wide interest only around 2015, startingwith the influential paper by Sawyer [20], which pointed out that the angular distri-butions for the different neutrino flavors are different in the decoupling region in asupernova and that can cause oscillations with rate ∝ µ . This was followed by twocrucial studies which developed theoretical understanding of fast oscillations in thelinearized regime [21, 22]. The first elucidated an underlying instability, showing howit has little dependence on ω E [21], and the second clarified the role of homogeneityand stationarity, and showed that fast oscillations may occur for the the presumablyrealistic neutrino distributions inspired by state-of-the-art simulations [22]. Subse-quently, this subject has seen significant progress [23–32], and it appears that fastoscillations can occur in many environments [33–40] and drastically alter the flavorcomposition therein [41, 42]. In particular, they may lead to almost complete flavordepolarization [41, 42], limited by conservation of lepton asymmetry.Properly understanding these new effects requires solving a set of coupled nonlin-ear partial integro-differential equations in a 3 (spatial) + 3 (momentum) + 1 (time)= 7 dimensional space, which has remained unachievable till date. As a result, moststudies assume a high degree of symmetry and limit the flavor evolution to be alonga single coordinate — be it temporal evolution of a homogeneous gas or the spatialevolution of a stationary dense neutrino gas. Previous experience with slow collectiveoscillations [16–19] suggests that such restrictions artificially exclude allowed modes offlavor evolution. Calculations in 1+1+1 dimensions are superior [41, 42], but as we willdiscuss in Sec. 2.2.1 they too are not fully satisfactory. While a calculation in 3+3+1dimensions remains the holy grail, going to at least 2+2+1 dimensions is necessary tosatisfy the constraints on neutrino velocities, include fully self-consistent space-timeevolution, and not artificially exclude allowed modes of flavor evolution.In this paper, we present a study of fast flavor oscillations of a dense neutrinogas in 2+2+1 dimensions. Our study explores how the flavor evolution is related tothe symmetries and the number of zero crossings of the electron lepton number (ELN)distribution, i.e., the difference of ν e and ¯ ν e angular distributions; a more precisedefinition appears later. We perform detailed comparisons of the numerical resultswith those from a linearized analysis, finding excellent agreement. The main insightsobtained from our study are that the flavor evolution inherits symmetries of the ELN,and why ELNs with a large number of zero crossings can in fact lead to a much weakerflavor instability. – 2 –he paper is organized as follows: Sec. 2.1 sets up the problem. Sec. 2.2 introducesthe types of neutrino angular distributions that we consider in this study. In Secs. 2.3,2.4, and 2.5, we describe the analytical and numerical methods we use for solving theequations. In Sec. 3, we discuss our numerical and analytical results for all examplesin a systematic way. In Sec. 4, we conclude with a brief summary. Consider an effective two-flavor framework with each neutrino being a superpositionof the e and µ flavors. Neglecting momentum-changing collisions, the flavor evolutionof the neutrinos at position (cid:126)x and time t , with momentum (cid:126)p ≈ E(cid:126)v , is given by [43] (cid:0) ∂ t + (cid:126)v · (cid:126)∂ (cid:1) P [ E, (cid:126)v ] = H [ E, (cid:126)v ] × P [ E, (cid:126)v ] . (2.1)We use the notation of Ref. [42]. The polarization vector P [ E, (cid:126)v ] for each momentumlabeled by (
E, (cid:126)v ) encodes the flavor state. Antineutrinos are represented by ¯ P [ E, (cid:126)v ].However, polarization vectors for antineutrinos behave as if they were polarizationvectors for neutrinos with negative E , so it is convenient to define P [ − E, (cid:126)v ] := − ¯ P [ E, (cid:126)v ],where now the argument E in P [ E, (cid:126)v ] takes values between −∞ to + ∞ . The overallminus sign is a notational foresight that makes Eq.(2.5) simpler.Neutrino oscillations do not change the total occupation numbers of neutrinos(or antineutrinos), but only the difference between flavors. One defines P [ E, (cid:126)v ] = g [ E, (cid:126)v ] S [ E, (cid:126)v ], with S [ E, (cid:126)v ] having a unit length and g [ E, (cid:126)v ] = (cid:40) + f ν e [+ E, (cid:126)v ] − f ν µ [+ E, (cid:126)v ] for
E > − f ¯ ν e [ − E, (cid:126)v ] + f ¯ ν µ [ − E, (cid:126)v ] for
E < E , but together they are packaged into a single function g [ E, (cid:126)v ] which spans over E ∈ ( −∞ , + ∞ ). Typically one has an excess of ν e over ν µ (resp. ¯ ν e over ¯ ν µ ) at anymomentum, because the electron flavors can be preferentially produced via chargedcurrent processes. This means that g [ E, (cid:126)v ] is positive (resp. negative) for neutrinos(resp. antineutrinos). If instead there is an excess of ν µ over ν e (resp. ¯ ν µ over ¯ ν e ) itsimply means that g [ E, (cid:126)v ] becomes negative (resp. positive). A few words about notation: Functional dependence will be denoted by square brackets [ . . . ],though dependence on space-time ( (cid:126)x, t ) is implicit. We reserve the parenthesis ( . . . ) for groupingterms together or to denote vectors written as a tuple of their components. Sans-serif letters suchas S represent three-component vectors in flavor space. Symbols in the usual italic with an arrow ontop, e.g., (cid:126)x , represent three-dimensional vectors in real space. Later in the paper we will introducematrices that are space-time tensors; these are written as bold symbols, e.g., Π , and their componentsas Π ij with subscripts ij being indices that run over space and time. Angle brackets, i.e., (cid:104)· · · (cid:105) , willstand for spatial averaging over all spatial coordinates. – 3 –n the flavor basis, the orientation (0 , , +1) corresponds to a purely electronflavor and (0 , , −
1) to a muon flavor. All neutrinos and antineutrinos start out asflavor eigenstates and the initial state of all Bloch vectors is chosen to be S ini [ E, (cid:126)v ] =(0 , , +1). In this convention, at any space-time point, the third component of theBloch vector S [ E, (cid:126)v ] is equal to twice the survival probability minus one.The Bloch vector for the hamiltonian has the form H = H vac + H mat + H self .Neutrino mass-mixing gives H vac [ E ] = ± ω E (sin 2 θ, , cos 2 θ ) , (2.3)with the plus (resp. minus) sign chosen for the normal (resp. inverted) mass ordering.Note that antineutrinos having been defined as neutrinos with negative E accountsfor the sign-flip needed in the mass-mixing hamiltonian. Further, the effect of forwardscattering on electrons in the background matter is encoded in H mat = λ (0 , , , (2.4)and H self [ (cid:126)v ] = √ G F (cid:90) d (cid:126)v (cid:48) (2 π ) (cid:0) − (cid:126)v · (cid:126)v (cid:48) (cid:1) (cid:90) + ∞−∞ E (cid:48) dE (cid:48) g [ E (cid:48) , (cid:126)v (cid:48) ] S [ E (cid:48) , (cid:126)v (cid:48) ] (2.5)is the neutrino-neutrino interaction term that depends on the flavor states of otherneutrinos and antineutrinos, and causes collective flavor oscillations. Note that theoverall minus sign chosen in the definition the polarization vectors for antineutrinoshas converted (cid:82) ∞ E dE ( P − ¯ P ) to (cid:82) + ∞−∞ E dE g [ E, (cid:126)v ] S [ E, (cid:126)v ].We will consider astrophysical scenarios where the neutrino density is large, i.e., µ (cid:29) ω E , λ , and neglect the vacuum and the matter terms from the hamiltonian. Inthis limit, inspection of the above equations shows that the dependence of S [ E, (cid:126)v ] onenergy drops out. So, it makes sense to rewrite Eq.(2.1) as (cid:16) ∂ t + (cid:126)v · (cid:126)∂ (cid:17) S [ (cid:126)v ] = µ (cid:90) d (cid:126)v (cid:48) (1 − (cid:126)v · (cid:126)v (cid:48) ) G [ (cid:126)v (cid:48) ] S [ (cid:126)v (cid:48) ] × S [ (cid:126)v ] , (2.6)where we note that S no longer depends on E but only on (cid:126)v , and the collective potentialis some constant µ = √ G F n ν over length and time scales of interest. We use µ − = 1as our unit of distance.The function G [ (cid:126)v ] is called the electron lepton number (ELN) distribution as afunction the direction (cid:126)v [21], G [ (cid:126)v ] = 1(2 π ) n ν (cid:90) ∞ E dE (cid:16) f ν e [ E, (cid:126)v ] − f ¯ ν e [ E, (cid:126)v ] − f ν µ [ E, (cid:126)v ] + f ¯ ν µ [ E, (cid:126)v ] (cid:17) (2.7a)= G ν e [ (cid:126)v ] − G ¯ ν e [ (cid:126)v ] − G ν µ [ (cid:126)v ] + G ¯ ν µ [ (cid:126)v ] . (2.7b)Here each term on the right hand side is the ratio of the number density of neutrinosor antineutrinos with that flavor and velocity (cid:126)v to the neutrino density n ν . Strictly speaking, one must keep ω E (cid:54) = 0 in H vac ; otherwise one can set the mass and flavor basisto be identical, once and for all, and there are no oscillations. Thus ω E affects the kickstarting of theoscillations, but the subsequent evolution very weakly [21, 23]. In practice, one imagines setting ω E (or an external perturbation used as its proxy) to zero immediately after the evolution begins. – 4 – igure 1 : Schematic representations of contributions to the ELN, with the four kindsangular distributions of the ν e and ¯ ν e and thus different numbers of crossings. ∞ Crossings
In a typical astrophysical scenario, where the ambient temperature is lower than themuon mass, there is no significant difference in the occupations for ν µ and ¯ ν µ . Thus G [ (cid:126)v ] is approximately equal to G ν e [ (cid:126)v ] − G ¯ ν e [ (cid:126)v ], and can be called the electron leptonnumber distribution. Integrating it over velocity, one finds A = (cid:90) d (cid:126)v G [ (cid:126)v ] = n ν e − n ¯ ν e n ν , (2.8)which is the net lepton asymmetry in neutrinos.We now focus on a broad feature of this function G [ (cid:126)v ] — the possibility thatthis function changes sign as a function of direction — a feature that is referred toas an ELN crossing. The reason we focus on this feature, is that it appears to benecessary (and perhaps sufficient) for causing fast flavor oscillations. A crossing canoccur if G ν e [ (cid:126)v ] and G ¯ ν e [ (cid:126)v ] have different velocity dependence and are of a comparablemagnitude. Broadly, one can think of four possible scenarios:– 5 – No crossing : If the density of ν e far exceeds that of ¯ ν e , or vice versa, then oneof the terms in G [ (cid:126)v ] dominates and there is no crossing in the ELN. The function G [ (cid:126)v ] remains positive or negative everywhere. This situation is shown as thetop-left schematic in Fig. 1. • One crossing : In the neutrino decoupling region in a SN, ν e is expected to havea higher number density compared to ¯ ν e . However, ¯ ν e kinematically decouples ata smaller radius, and in the forward direction one may expect the ¯ ν e contributionto the ELN to exceed the ν e [20], if there are directions along which the leptonasymmetry is not too large [22]. In this case, the function G [ (cid:126)v ] is negativein the forward direction, but positive elsewhere, and there is a single crossing.Hydrodynamic fluctuations of the lepton asymmetry can allow asymmetry tobe small along some direction(s) and ELN crossing(s) could occur [34–36]. Asanother example, in the disk of a neutron star merger [33, 38] (or if a quark-hadron phase transition occurs in a SN [44]) the number density of ¯ ν e can exceedthat of ν e , and all else remaining the same, a ¯ ν e excess in the radially outwarddirection is even more likely. This situation is shown as the top-right schematicin Fig. 1. • Two crossings : If the ¯ ν e contributions to the ELN exceed the ν e in the forwardand backward directions, the function G [ (cid:126)v ] is negative in the forward and back-ward directions, but positive in the directions tangential to the radial direction.A possible cause of larger ¯ ν e numbers in the backward direction, in addition toa forward excess as discussed above, could be that ¯ ν e have a larger cross sec-tion to backscatter off nuclei, thus creating a dominantly ¯ ν e back-flux [34]. Thissituation with two crossings is shown as the bottom-left schematic in Fig. 1. • Many crossings : If the ν e and ¯ ν e contributions to the ELN are almost equalbut fluctuate independently, perhaps due hydrodynamic waves and instabilities,one would expect G [ (cid:126)v ] to be changing sign frequently as a function of direction.In this scenario, the ELN has many crossings. This is likely in the convectivelayer of a SN [35–37], and perhaps in the early Universe [40]. A schematic of thissituation is shown on the bottom-right in Fig. 1. The figures in the schematic in Fig. 1 are best visualized as sections of the correspondingsurfaces in the three-dimensional velocity space. However, studying fast flavor oscil-lation in three spatial dimensions is numerically challenging. Typically, one assumesthat G [ (cid:126)v ] is azimuthally symmetric, i.e., invariant to its rotations in velocity-spaceabout the radially outward direction (here coinciding with ˆ x ). This is the picture thathas been adopted in almost all studies. Mathematically, this requires dropping v z and v y everywhere and setting (cid:126)v · (cid:126)∂ → v x ∂ x and (cid:126)v · (cid:126)v (cid:48) → v x v (cid:48) x in Eq. (2.6). Thus, in moststudies, even when one solves the counterpart of Eq.(2.1) in some lower number of spa-tial dimensions, say d = 1 assuming azimuthal symmetry about the radially outwarddirection, what one has in mind is the higher dimensional version, say with D = 3,– 6 –nd the solution is assumed to be strictly symmetric with respect to the remaining D − d = 2 components of the velocity, v y and v z . In this approach v x does not have aunit magnitude, and one only requires | v x | <
1. Although it is prevalent and appealingas a modeling simplification, it is an uncontrolled approximation. In previous studieson slow collective oscillations [16–19], it was shown that additional instabilities canbe generated by spontaneous symmetry breaking along the assumed-to-be-symmetricdimensions.In this paper we will consider a neutrino gas in strictly two spatial dimensions, x and y , with velocities restricted to a two-dimensional plane, i.e., (cid:126)v = v x ˆ x + v y ˆ y with v x + v y = 1. This should not be thought of as a projection of a three dimensionalproblem on to two dimensions, as described in the paragraph above. In this way, thevelocity vector has a unit length and no instabilities are ignored by fiat. The downsideis that we are solving a two-dimensional toy problem, and it may not be obvious howit applies to the three-dimensional real-world. Although our calculation in 2+2+1dimensions is the state-of-the-art, it should really be taken as a step towards morerealistic studies in full 3+3+1 dimensions. Still, this set-up has its merits as discussedabove, and we will find important insights by undertaking this calculation. In this strictly two-dimensional approach, the ELNs in the figure are functions of asingle independent variable, say θ , and one writes G [ (cid:126)v ] = G [ θ ], where v x = cos θ and v y = sin θ . In the remainder of our study, we consider three out of the four zero crossingscenarios shown in Fig.1. We ignore the case with no crossings because one finds noinstabilities in that case. We further simplify our ELNs to be piecewise constant orsinusoidal, for concreteness, and consider • Type I (single crossing): G [ θ ] = (cid:40) A − π , if v x = cos θ > π , if v x = cos θ < . • Type II (two crossings): G [ θ ] = A − π , if v x = cos θ > v y = sin θ > π , if v x = cos θ < v y = sin θ > A − π , if v x = cos θ < v y = sin θ < π , if v x = cos θ > v y = sin θ < . • Type III (many crossings): G [ θ ] = A π + c cos mθ + c sin mθ .The above choices are made in a way such that G [ θ ] remains unchanged under v y → − v y for Type I, and a simultaneous exchange v y → − v y and v x → − v x in caseof Type II. These choices are made to explore the dependence of the solution on the One point of semantics before we proceed further: For a three-dimensional but azimuth-symmetricELNs these correspond to one closed curve , two closed curves , and many closed curves worth ofcrossings. For our strictly two-dimensional ELNs we will continue to call the above ELNs as havingone, two, or many crossings, though strictly speaking they have two, four, and twice as many crossings,respectively. – 7 –ature of symmetry of G [ θ ]. For Type III, m will be taken to be large and G [ θ ] canhave O ( m ) number of zero crossings as a function of θ . In all these cases, A = (cid:90) π dθ G [ θ ] (2.9)denotes the lepton asymmetry. Initially the transverse component of the Bloch vector is small, i.e., S ⊥ [ v x , v y ] (cid:28) ∂ t + v x ∂ x + v y ∂ y ) S ⊥ [ v x , v y ] = i (cid:90) +1 − (cid:90) +1 − dv (cid:48) x dv (cid:48) y δ [ v (cid:48) − (cid:0) − v x v x (cid:48) − v y v (cid:48) y (cid:1) G [ v (cid:48) x , v (cid:48) y ] × (cid:0) S ⊥ [ v (cid:48) x , v (cid:48) y ] − S ⊥ [ v x , v y ] (cid:1) . (2.10)where v = (cid:112) v x + v y . This is a linear equation in S ⊥ . It is thus natural to decomposeit in the Fourier basis, S ⊥ [ v x , v y ] = (cid:88) (cid:126)K, Ω Q ⊥ (cid:126)K, Ω [ v x , v y ] e i ( K x x + K y y − Ω t ) , (2.11)which gives a linear algebraic equation connecting Ω with (cid:126)K . Once this relationshipΩ( (cid:126)K ) is known, it gives a set of basis functions Q labelled by ( (cid:126)K, Ω) which can belinearly superposed to describe any solution of S ⊥ [ v x , v y ] in the linear regime.Concretely, one finds the dispersion relation [21, 22, 46–52] D = det Π [ k x , k y , ω ] = 0 , (2.12)where Π [ k x , k y , ω ] = η + (cid:90) − (cid:90) − dv x dv y G [ v x , v y ] ω − v x k x − v y k y δ [ v − W [ v x , v y ] . (2.13)If there is a solution that grows exponentially, e.g., Im ω >
0, that solution is saidto be unstable [47]. A more detailed classification can be found in Refs. [49, 50]. InEq.(2.13) the following definitions have been used: η = − − , (2.14)and W [ v x , v y ] = v x v y v x v x v x v y v y v x v y v y , (2.15)– 8 – able 1 : Elements of φ for various types of ELNsELN φ tt φ tx φ ty φ xx φ yy φ xy γ γ γ Type I Aπ (2 − A )4 A π A π A π + ( A −
16 2 A π - A ( A − π Type II Aπ A π A π A − π − A +4 A − π A ( A +4 A − π Type III A 0 0 A A − A A and ω = Ω − φ tt , (2.16a) k x = K x − φ tx , (2.16b) k y = K y − φ ty , (2.16c)wherein φ = (cid:90) − (cid:90) − dv x dv y G [ v x , v y ] δ [ v − W [ v x , v y ] . (2.17)Thus, in the linear regime, allowed solutions are given by Eq.(2.12), which uponafter expanding gives − (Π ty ) Π xx + 2Π tx Π ty Π xy − Π tt (Π xy ) − (Π tx ) Π yy + Π tt Π xx Π yy = 0 (2.18)For each pair ( k x , k y ) one needs to solve Eq.(2.18) for ω , whose imaginary part describesthe growth of the flavor instabilities. Note that because the only dimensionful quantityin the problem is µ , these are all fast instabilities, i.e., Im ω ∝ µ . In general, theabove equation is transcendental and analytical solution is impossible. However, for( k x = 0 , k y = 0), Eq.(2.18) becomes a simple cubic equation in ω = ω [ k x = 0 , k y = 0], ω + γ ω + γ ω + γ = 0 , (2.19)where γ = φ tt − φ xx − φ yy , (2.20a) γ = ( φ tx ) + ( φ ty ) − φ tt φ xx − ( φ xy ) − φ tt φ yy + φ xx φ yy , (2.20b) γ = − ( φ ty ) φ xx + 2 φ tx φ ty φ xy − φ tt ( φ xy ) − ( φ tx ) φ yy + φ tt φ xx φ yy . (2.20c)We remind the reader that the ω here is not | ∆ m | / (2 E ), but merely the zerothcomponent of the Fourier mode in Eq.(2.16a). The values of ( γ , γ , γ ) and φ ij for ourchosen Type I, II, and III neutrino angular distributions are listed in Table 1.– 9 – .4 EoM Solver We developed our own numerical routines for solving Eq.(2.6). Our approach involvesdiscretizing the spatial directions into N x and N y uniformly spaced bins, resulting in N x N y number of coupled nonlinear ODEs in time for each momentum mode labeled byits ( v x , v y ) pair. A periodic boundary condition in each spatial direction is assumed. Wealso discretize the velocity modes in one direction (either v x or v y ) into N vel uniformlyspaced bins, and due to the restriction v = 1 the binning of the velocity modes inthe other direction gets fixed resulting in total of (3 × N x N y N vel coupled nonlinearODEs. The factor of 3 comes from the three components of the polarization vectorand the factor of 2 comes from the fact that for each choice of v x one has two allowedchoices of v y .We solve the system equations in a 2D square box of area L × L , with L = 18 inunits of µ − . We choose µ = 3 π cm − for Type I ELNs, and for Type II ELNs we takeeither µ = 3 π cm − or µ = 17 π cm − which correspond to a neutrino number densityof O (10 ) cm − . Periodic boundary conditions are assumed on both spatial directions,i.e., on x, y ∈ (cid:0) − L , L (cid:1) . This periodic boundary condition physically represents thatwe are treating this box as a part of a larger system. The finiteness of the box affectsthe smallest ∆ (cid:126)k we can distinguish in our calculation. We discretize x and y into N x = N y = 480 bins, enough to trigger as many Fourier modes as possible, and wellabove what is needed to trigger all unstable (cid:126)k modes, limited only by CPU hours. Thevelocity of outgoing neutrinos are in the range v x , v y ∈ ( − ,
1) with N vel = 32. Intotal, we solve a system of 6 N x × N y × N vel = 44236800 coupled nonlinear ODEs intime up to t fin = 3 . µ − . The choices for N x , N y , N vel are optimized toobtain sufficient precision and accuracy as shown in Appendix A.The initial conditions are that all Bloch vectors S [ (cid:126)v ] are equal to (0 , , G [ (cid:126)v ], the polarizationvectors P [ (cid:126)v ] = G [ (cid:126)v ] S [ (cid:126)v ] points along (or opposite to) the vertical in flavor space. Forour chosen set of ELNs, as discussed in Sec. 2.2, this means that the P [ (cid:126)v ] start withone, two, or many crossings as a function of θ . Normally H vac ω would start the flavorevolution by tilting the Bloch vectors away from their initial positions. However, forour calculations, we have set H vac ω and H mat to zero for numerical convenience, andinstead provide an external perturbation of O (10 − ) to both transverse components ofthe polarization vectors at ( x = 0 , y = 0) to kickstart the evolution.The code is written in Python and uses the zvode solver, a variable-coefficientdifferential equation solver in
Python , to solve the system of ODE as a functionof time. This solver implements the backward differentiation formula for numericalintegration. Our technique of converting a set of coupled nonlinear PDEs into ODEsallows easy use of existing ODE libraries and makes the numerical integration muchfaster. The spatial derivatives are computed using a Fast Fourier Transform employing
Python ’s scipy.fftpack.diff package. See the footnote in Sec. 2.2, clarifying what we mean by one, two, and many crossings in thecontext of our ELNs. – 10 – .5 Dispersion Relation Solver
To understand the nature of the flavor evolution in the linear regime one needs toknow the behavior of ω as a function of (cid:126)k , which we represent by its magnitude k andargument as β , i.e., k = (cid:112) k x + k y , (2.21) β = arctan (cid:16) k y k x (cid:17) . (2.22)This ω [ k, β ] can only be obtained after solving the transcendental equation describedby Eq.(2.18). This is not easy, even numerically, as one needs to scan over a largespace spanned by ( ω, k, β ), and brute force root-finding is inefficient.To speed up our root-finding, we use the method of iterative solving where webegin at a known solution (or initial guess) and then use it to propagate the solutionfurther in ( ω, k, β ) space. As our initial guess we use the k = 0 solution, ω , which fora given ELN can be determined easily from Eq.(2.19). This is simply the zero-modesolution that was advocated in Ref.[25], and even in the most general case with 3+3+1dimensions Eq.(2.19) is analytically tractable. Then we define a circular boundary ofradius r , chosen to be sufficiently small and close to k = 0 point in the k − β plane.Using w as our initial guess we numerically solve Eq.(2.18) to calculate ω for different β directions within the region 0 < k < r in the k − β plane. Then we proceed to a newpoint on the boundary defined by r , where we already have a solution, to define anothercircle of radius r . At each new point within this new circular region we can start withprevious solution as a starting guess, and find the updated solutions. Note that wechoose r in a way such that the previous guess works reasonably well. Repeating this,we can find the solution on the entire k − β plane.The roots of Eq.(2.18) are obtained by Python’s fsolve package which usesPowell’s conjugate direction method to find the local minima of a nonlinear equa-tion. The method requires an initial guess, but does not require differentiability ofthe underlying complex function because no derivatives are computed in order to findthe solution. The integrals in Eq.(2.13) are evaluated using the numerical routine foradaptive quadrature implemented in
Python’s quad solver.
For our Type I ELNs, G [ θ ] has a reflection symmetry for v y = sin θ , i.e., G [ θ ] remainsinvariant under v y → − v y . This symmetry leads to a similar symmetry in the k x − k y plane, i.e., D [ k x , − k y , ω ] = D [ k x , k y , ω ]. This can be understood considering theinterchange v y → − v y and k y → − k y in Eq.(2.13):Π ij [ k x , k y , ω ] k y → − k y −−−−−→ η ij + (cid:90) − (cid:90) − dv x dv y G [ v x , v y ] ω − k x v x + k y v y δ [ v − W ij [ v x , v y ] v y → − v y −−−−−→ ± Π ij [ k x , k y , ω ] , (3.1)– 11 – ype I, A = 0 β = 0 β = π β = π β = π β = πβ = π β = π β = π Im ω = 8
Im ω = 6
Im ω = 4
Im ω = 2 k = 6 k = 9 k = 13 I m ω [ k , β ] β = π2 , β = 0β = π Type I, A (cid:54) = β = 0 β = π β = π β = π β = πβ = π β = π β = π Im ω = 8
Im ω = 6
Im ω = 4
Im ω = 2 A = 0 . A = 0 I m ω [ k , β ] A = 0.4A = 0
Figure 2 : Left: Angular variation of Im ω with respect to β . Right : Radial variationof Im ω with respect to k . The top plots are done with A = 0 for different values of k (left) and β (right) whereas the bottom ones with fixed k = 9 (left) and β = 0 (right)for A = 0 . A = 0 case are shown in blue dashedlines in the bottom panel plots.where we have used the fact that G [ v x , v y ] = G [ v x , − v y ] in the last step. Eq.(3.1)basically says that Π ij remains invariant under the above two operations up to a ± sign. The minus sign occurs only for Π ty and Π xy while all others come with a plussign. However, Π ty and Π xy always come in pairs in Eq.(2.18), i.e., as (Π ty ) , (Π xy ) or Π ty Π xy , which immediately says that D [ k x , k y , ω ], as well as the solution for Im ω ,remains invariant under k y → − k y . – 12 –q.(3.1) implies for k = 0 or ( k x = k y = 0) mode:Π ty [0 , , ω ] = 0 (3.2)and Π xy [0 , , ω ] = 0 . (3.3)Eq.(3.2) and Eq.(3.3) help us to write the full dispersion relation in Eq.(2.18) as twoseparate equations: Π yy [0 , , ω ] = 0 (3.4)and Π tx [0 , , ω ]Π tx [0 , , ω ] − Π tt [0 , , ω ]Π xx [0 , , ω ] = 0 . (3.5)Eq.(3.4) implies ω = φ yy which is a real solution. Eq.(3.5) can be simplified to obtaina quadratic equation in ω as, ω + ω (cid:16) φ tt − φ xx (cid:17) − (cid:16) φ tt φ xx − ( φ tx ) (cid:17) = 0 (3.6)The solutions determined by Eq.(3.6) can be complex only if (cid:0) φ tt + φ xx (cid:1) − φ tx ) < . (3.7)Eq.(3.6) can have complex solutions, in general, leading to the k = 0 mode becomingunstable for Type I cases. For instance, in our numerical examples with A = 0 (resp. A = 0 .
4) the LHS of Eq.(3.7) becomes − − . k = 0 solution, ω , we can compute ω for othervalue of (cid:126)k using the iterative method described in Sec. 2.5. On the other hand, we canalso numerically simulate Eq. (2.6) to obtain S [ (cid:126)v ] as a function of ( x, y, t ). We thentake spatial Fourier transforms of this solution, and ask how the amplitude of each (cid:126)k mode changes with time. For some modes, we find the mode-amplitude increasesexponentially, and we extract the imaginary part of ω ( (cid:126)k ) from numerical data. Thesetwo methods give results in excellent agreement, as shown in Fig. 2.For Type I ELNs, G [ θ ] is symmetric between the regions sin θ > θ < v y → − v y symmetry. This gives rise to a similar symmetryin the angular variation of Im ω with respect to β , i.e., Im ω | k y > = Im ω | k y < as canbe seen in the top left panel plot of Fig.2. The radial variation of Im ω with respectto k in the top right plot of Fig.2 clearly shows a much larger growth for the modesvery close to k = 0, and then the growth rate decreases as a function of k for all β directions. The k = 0 mode being unstable for this case is understood from ouranalytical arguments. We find the decrease is much slower along the k y axis, aboutwhich there is a k y → − k y symmetry. Interestingly the position of the Fourier modewith the maximum linear growth rate is aligned along k x axis, as shown in the blackstarred point in top left plot of Fig. 2.Even for A (cid:54) = 0 a similar kind of k y → − k y symmetry in the angular variation ofIm ω is shown in the bottom left plot of Fig. 2. The k = 0 mode is unstable for this– 13 –ase as well. All the results for the radial and angular variation of Im ω are similarto the A = 0 case with only an exception that the overall growth rate as well as themaximum growth rate decreases for larger (positive) A . This is seen in the bottomleft panel of Fig. 2, where the position of the maximum growth rate is indicated bydifferent starred points that correspond to specific choices of A , as indicated by thecolor code. This effect of non-zero lepton asymmetry also results in a faster decreaseof Im ω with larger k , as shown in bottom right plot of Fig. 2. Type II ELNs have a symmetry under the joint operations v x → − v x and v y → − v y ,for any value of A . This results in a D [ − k x , − k y , ω ] = D [ k x , k y , ω ] type of symmetryin k x − k y plane. This can be understood throughΠ ij [ k x , k y , ω ] k y → − k y −−−−−→ k x → − k x η ij + (cid:90) − (cid:90) − dv x dv y G [ v x , v y ] ω + k x v x + k y v y δ [ v − W ij [ v x , v y ] v y → − v y −−−−−→ v x → − v x ± Π ij [ k x , k y , ω ] . (3.8)In the last step of Eq.(3.8), G [ − v x , − v y ] = G [ v x , v y ] has been used. Eq.(3.8) saysthat Π ij remains invariant under the joint operations of k x → − k x and k y → − k y except a minus sign that only occurs for Π tx and Π ty . But interestingly again theycome in pairs in Eq.(2.18), such as (Π tx ) , (Π ty ) or Π tx Π ty , leaving the dispersionrelation as well as as the solution for Im ω invariant under the above operations. Aninteresting special case occurs for Type II ELNs with A = 0 where now the dispersionrelation can have two more symmetries. For example, D [ − k x , k y , ω ] = D [ k x , k y , ω ] and D [ k x , − k y , ω ] = D [ k x , k y , ω ] along with the previous one. These extra symmetries for A = 0 case can be easily understood using similar arguments as above.Eq.(3.8) for k x = k y = 0 gives Π tx [0 , , ω ] = 0 (3.9)and Π ty [0 , , ω ] = 0 . (3.10)This further simplifies Eq.(2.18) into two separate equations:Π tt [0 , , ω ] = 0 (3.11)and Π xy [0 , , ω ]Π xy [0 , , ω ] − Π xx [0 , , ω ]Π yy [0 , , ω ] = 0 . (3.12)Eq.(3.11) implies ω = − φ tt which is a real solution. Eq.(3.12) simplifies to giverise to a quadratic equation, ω − ω (cid:16) φ xx + φ yy (cid:17) + (cid:16) φ xx φ yy − ( φ xy ) (cid:17) = 0 . (3.13)– 14 – ype II, A = 0 β = 0 β = π β = π β = π β = πβ = π β = π β = π Im ω = 8
Im ω = 6
Im ω = 4
Im ω = 2 k = 18 k = 22 k = 30 I m ω [ k , β ] β = π2 β = π6 β = 0 Type II, A (cid:54) = β = 0 β = π β = π β = π β = πβ = π β = π β = π Im ω = 8
Im ω = 6
Im ω = 4
Im ω = 2 A = 0 . A = 0 I m ω [ k , β ] A=0.1A=0
Figure 3 : Left: Angular variation of Im ω with respect to β . Right : Radial variationof Im ω with respect to k . The top plots are done with A = 0 for different values of k (left) and β (right) whereas the bottom ones with fixed k = 30 (left) and β = π (right)for A = 0 .
1. In these plots the continuous lines represent the linear stability solutionwhile the dotted points the solution of the full equation of motion. For comparisonpurpose the solution for A = 0 is also shown in green dashed lines in the bottom panelplots.The solutions of Eq.(3.13) are complex only if (cid:0) φ xx − φ yy (cid:1) + 4 (cid:0) φ xy (cid:1) < . (3.14)The condition in Eq.(3.14) can never be fulfilled, as the LHS is a sum of two perfectsquares, thus implying a stable zero-mode for Type II ELNs.Fig. 3 shows the angular variation of Im ω predicted by our previous analyticalarguments for this case. The stability of k = 0 mode as shown in the radial variation– 15 –n top-right plot of Fig.3, confirms our analytical claim. Already one finds that havingtwo crossings leads to less instability in some sense. Interestingly, the radial variationof Im ω for this case shows an approximately Lorentzian shape as a function of k ,i.e., large wavelength or small k modes are inert then the growth rate increases as weincrease k with a maximum around k = 30 −
40 and then starts to decrease as a resultvery small wavelength or very large k modes again become inert. The Lorentzian ismuch wider closer to the β = π/ β = 0. The Fourier modes closeto k y ( β = π/
2) axis have much larger growth rate compared to modes close to k x ( β = 0) for this case. In contrast with the Type I case, the Fourier mode with thelargest growth rate in this case lies along the k y axis ( β = π ) about which there is a k y → − k y symmetry.For A (cid:54) = 0, the angular variation of Im ω shown in the bottom left panel ofFig.3 shows a skewed symmetry. The non-zero value of lepton asymmetry breaks thesymmetry along k x and k y axes keeping the symmetry along the diagonals intact andalso tilts the overall angular distribution towards one of the diagonals. The non-zerovalue of A also shifts the position of the Fourier mode with the largest growth ratefrom the k y axis to along one of the diagonals. This plot also indicates that the Fouriermodes along β = π/ β = 0direction but with an exception that in this case the overall growth of the system issuppressed compared to zero lepton asymmetry case. The k = 0 mode is stable in thiscase also, as shown in bottom right plot of Fig. 3. The same plot also indicates thatIm ω as a function of k for specific β direction shows a similar Lorentzian nature butit is slightly shifted towards k = 0 with lower width compared to A = 0 case. Thediagonal symmetry, the shift in the position of the maximum of Im ω and the decreasein overall growth of the system become more pronounced as we increase the value oflepton asymmetry. Type III ELN corresponds to a scenario where, e.g., due to fluctuations of the neutrinodistributions, the angular distributions of ν e , ν e are rapidly changing as a function of θ and their difference can go through large number of zero crossings. To mimic thisscenario we considered the Type III G [ θ ] with m being quite large. We will show howthe maximum growth of such systems depends on the number of zero crossings, closelyrelated to m . To understand this let us consider Eq.(2.13) in terms of the angularvariable θ and the k − β coordinates, Π [ k, β, ω ] = η + ψ [ ω, k, β ] = η + (cid:90) π dθ f [ θ, ω, k, β ] G [ θ ] W [ θ ] , (3.15)where we have defined ψ [ ω, k, β ] as the matrix of integrals on the RHS. As before, butnow explicitly in terms of θ , one has f [ θ, ω, k, β ] = 1 ω − k cos ( θ − β ) (3.16)– 16 – m . . . . I m e ω m a x Numerical m fit (Analytical) Figure 4 : Im (cid:101) ω max as a function of the crossing number m is shown for G [ θ ] = 6 sin mθ .Analytical arguments show that Im (cid:101) ω max decreases as 1 /m .and W [ θ ] = θ sin θ cos θ cos θ sin θ cos θ sin θ sin θ cos θ sin θ . (3.17)Eq.(3.16) indicates that the integrands in each component of ψ [ ω, k, β ] has a saddle-point at say θ [ i, j ]. In general, depending on W ij [ θ ], the stationary points θ [ i, j ]will be shifted away from cos − (cid:0) Re ωk (cid:1) + β , and will not be the same for the different( i, j ) components. But from now on for convenience we will stop explicitly writingthe functional dependence on i, j in θ . Our objective, will be to compute the integralin Eq.(3.15) in the saddle-point approximation, and show how the Fourier mode withmaximum growth, labeled by ( ω max , k max , β max ), has smaller instability for larger m .First we rewrite the ij th component of ψ [ ω max , k max , β max ] in terms of the loga-rithm of its integrand as, ψ ij [ ω max , k max , β max ] = (cid:90) π dθ f [ θ, ω max , k max , β max ] G [ θ ] W ij [ θ ] = (cid:90) π dθ exp[ F ij [ θ ]] . (3.18)We then expand F ij [ θ ] about its stationary point θ up to second order to obtain, F ij [ θ ] = F ij [ θ ] + ( θ − θ ) d F ij [ θ ] dθ (cid:12)(cid:12)(cid:12)(cid:12) θ = θ . (3.19)Note F ij [ θ ] here is a complex function. Eq.(3.18) and Eq.(3.19) allow us to performthe saddle-point integral around θ to get ψ ij [ ω max , k max , β max ] = f [ θ , ω max , k max , β max ] G [ θ ] W ij [ θ ] (cid:115) π − d F ij [ θ ] dθ (cid:12)(cid:12) θ = θ . (3.20)– 17 –ow, one can insert the expression of Type III ELN and use the large- m limit to write − d F ij [ θ ] dθ (cid:12)(cid:12)(cid:12)(cid:12) θ = θ = m (cid:16) c + c + A π G [ θ ] − A π (cid:17) G [ θ ] + O ( m ) , (3.21)where c , are the coefficients of the sin and cos terms in the Type III ELN. The O ( m )can be neglected for Type III ELNs using the fact that m is large. Eq.(3.21) allows usto approximately simplify Eq.(3.18) to ψ ij [ ω max , k max , β max ] = (cid:101) G [ θ ] W ij [ θ ] m (cid:101) ω max , (3.22)where (cid:101) G [ θ ] = √ π (cid:0) c + c + A π G [ θ ] − A π (cid:1) / G [ θ ] , (3.23)and (cid:101) ω max = ω max − k max cos ( θ − β max ) (3.24)is the complex growth rate for the mode ( k max , β max ). Note (cid:101) ω max in principle canhave dependence on i, j indices via θ but as it appears only in the real part of (cid:101) ω max and in the limit θ [ i, j ] for different i, j indices are close to cos − (cid:16) Re ω max k max (cid:17) + β max orRe (cid:101) ω max ≈
0, it can be ignored.Equation (3.22) shows that (cid:101) ω max , and thus ω max , always appears multiplied by m . It is then obvious that Im ω max must scale as 1 /m . We checked this behaviorby numerically solving the dispersion relation and then locating the maximum of thissolution in the k − β plane. For Type III ELN with A = 0 , c = 1 , c = 6, and m inthe range 20 −
40 for which the ELN has many zero crossings, this numerical result isshown as the blue continuous line in Fig. 4. The black dashed line is the best fit of thatnumerical data with a m + a where the fitted parameter values are ( a , a ) = (13 . , . /m . This is of course obtained witha very particular form of the Type III ELN, which is purely sinusoidal with a single m .Based on numerical experiments, we conjecture that growth rates with many crossingsshould be small in general, all else being equal. Generalizing the above analysis for anarbitrary ELN, say written using a sine and cosine series in θ , does not seem straight-forward. In this paper we explored, analytically and numerically, how the initial growth of flavorinstabilities of a dense neutrino gas depends on the number of zero crossings in theELN. Improving upon previous lower-dimensional studies, this is the first study in 2(space) + 2 (momentum) + 1 (time) dimensions. We developed our own code thatsolves for the flavor evolution of dense neutrinos. We also presented a new strategy to– 18 –olve the transcendental equations appearing in the dispersion relation for the linearevolution of such systems. With these new tools we explored the linear behavior bylooking at the different radial and angular distributions of the dispersion relation inthe k − β plane. Our main results are • The symmetries of the Im ω in the k − β plane, its radial (i.e., vs. k ) andangular (i.e., vs β ) variation, the stability of the k = 0 mode, overall lineargrowth and the position of the Fourier mode with the highest growth rate, etc.,all have an intimate connection with the various symmetries of the neutrinoangular distributions and one can analytically understand them in great detail.Figs. 2, 3, 4 show the exquisite match between growth rates predicted by linearstability analysis and fully numerical evaluation of the solutions. This matchingand understanding over a variety of lepton asymmetries, A , shapes of ELNs,different wavelength k and directions β , shows the extraordinary power of lineartheory and a testament to the fidelity of our numerics. • ELNs with large number of zero crossings lead to a relatively smaller growthrate, essentially decreasing as 1 /m where m is the number of crossings. Wespeculate that this may be important for many realistic environments where ν e and ¯ ν e distributions are close to each other and crossings occur in the ELNs dueto noise or fluctuations. It seems that the growth rates for such instabilities willbe relatively suppressed. Acknowledgements
The work of B.D. is supported by the Dept. of Atomic Energy (Govt. of India) re-search project under Project Identification No. RTI 4002, the Dept. of Science andTechnology (Govt. of India) through a Swarnajayanti Fellowship, and by the Max-Planck-Gesellschaft through a Max Planck Partner Group.
A Error Estimate for Numerical Solutions
In Fig. 5 (left panel), we show the accuracy expected of our calculation. We check ifthe length of the polarization vectors remain fixed at unity or not. The error we incuron this is a lower bound on the error in our calculations. We find that our chosendiscretization, N x = N y = 500 and N v = 64, does as well as finer discretizations,incurring an error of O (10 − ) at t ≈ − . For illustration,here we have chosen a Type II ELNs with nonzero lepton asymmetry.In Fig. 5 (right panel), we illustrate the precision to be expected of our numericalsolutions of the equations of motion. We check for convergence by comparing thelength of the spatially averaged version of the polarization vector perpendicular tothe z -axis between two different discretizations: one with N x = N y = 500 and theother N x = N y = 1000. The computations are shown with N vel = 64. Our results– 19 – . . . . . t − − − − − m a x ( x , y ) { | S v x , v y | − } N vel = 64 N x = N y = 500 N x = N y = 1000 . . . . . t − − − − − D | S ⊥ N s p = − S ⊥ N s p = | E N vel = 64 Figure 5 : Left: Accuracy of our calculation, estimated by maximum departure of | S | from unity in the x − y plane, as a function of time t . Right: Precision of ourcalculation, estimated by convergence of (cid:104)| S ⊥ |(cid:105) with respect to our best discretization,i.e., ( N x = N y = 1000), shown as a function of time t . In all these calculations, wechoose N vel = 64 and for these plots we show the mode { v x = 0 . , v y = 0 . } . t − L − L L L y x = − L − − − − log | S ⊥ | t − L − L L L y x = 0 − − − − log | S ⊥ | t − L − L L L y x = L − − − − log | S ⊥ | t − L − L L L x y = − L − − − − log | S ⊥ | t − L − L L L x y = 0 − − − − log | S ⊥ | t − L − L L L x y = L − − − − log | S ⊥ | Figure 6 : Top : Flavor evolution in the y − t plane obtained from the numericalsolution of Eq.(2.6) for Type II ELNs with A = 0 .
1, shown at three different x positions.Bottom : Same in the x − t plane at three y positions. The color-coding in the colorbardepicts the log of the magnitude of S ⊥ .indicate that a discretization of N x = N y = 500 is at most O (10 − ) off from yet finerdiscretizations in the linear regime which ends almost at t = 2. This result is alsoshown for the same velocity mode with the same choice of ELN.For completeness, we show the flavor evolution on the x − t plane (resp. y − t – 20 –lane) at three y (resp. x ) positions for the above-considered case in Fig. 6. Note theoverall growth in flavor along y direction is much larger compared to x direction and itdecreases from the center towards the edge of the box in either directions, as dictatedby the fastest growing (cid:126)k mode. The box has been chosen to be much larger than theregion where the solution is nontrivial; this is to avoid artifacts of the finiteness of thebox. References [1] L. Wolfenstein,
Neutrino oscillations in matter , Phys. Rev. D (1978) 2369.[2] S. P. Mikheev and A. Yu. Smirnov, Resonance Oscillations of Neutrinos in Matter , Sov. Phys. Usp. (1987) 759.[3] J. T. Pantaleone, Neutrino oscillations at high densities , Phys. Lett. B (1992) 128.[4] V. A. Kostelecky and S. Samuel,
Selfmaintained coherent oscillations in denseneutrino gases , Phys. Rev.
D52 (1995) 621 [ hep-ph/9506262 ].[5] S. Pastor, G. G. Raffelt and D. V. Semikoz,
Physics of synchronized neutrinooscillations caused by selfinteractions , Phys. Rev.
D65 (2002) 053011[ hep-ph/0109035 ].[6] H. Duan, G. M. Fuller, J. Carlson and Y.-Z. Qian,
Simulation of Coherent Non-LinearNeutrino Flavor Transformation in the Supernova Environment. 1. CorrelatedNeutrino Trajectories , Phys. Rev. D (2006) 105014 [ astro-ph/0606616 ].[7] R. Sawyer, Speed-up of neutrino transformations in a supernova environment , Phys.Rev. D (2005) 045003 [ hep-ph/0503013 ].[8] S. Hannestad, G. G. Raffelt, G. Sigl and Y. Y. Wong, Self-induced conversion in denseneutrino gases: Pendulum in flavour space , Phys. Rev. D (2006) 105010[ astro-ph/0608695 ].[9] G. G. Raffelt and A. Y. Smirnov, Self-induced spectral splits in supernova neutrinofluxes , Phys. Rev. D (2007) 081301 [ ].[10] G. L. Fogli, E. Lisi, A. Marrone and A. Mirizzi, Collective neutrino flavor transitionsin supernovae and the role of trajectory averaging , JCAP (2007) 010 [ ].[11] A. Esteban-Pretel, S. Pastor, R. Tomas, G. G. Raffelt and G. Sigl, Decoherence insupernova neutrino transformations suppressed by deleptonization , Phys. Rev. D (2007) 125018 [ ].[12] B. Dasgupta and A. Dighe, Collective three-flavor oscillations of supernova neutrinos , Phys. Rev. D (2008) 113002 [ ].[13] B. Dasgupta, A. Dighe, G. G. Raffelt and A. Y. Smirnov, Multiple Spectral Splits ofSupernova Neutrinos , Phys. Rev. Lett. (2009) 051105 [ ].[14] H. Duan, G. M. Fuller and Y.-Z. Qian,
Stepwise spectral swapping with three neutrinoflavors , Phys. Rev. D (2008) 085016 [ ]. – 21 –
15] B. Dasgupta, A. Dighe, A. Mirizzi and G. G. Raffelt,
Spectral split in promptsupernova neutrino burst: Analytic three-flavor treatment , Phys. Rev.
D77 (2008)113007 [ ].[16] G. Raffelt, S. Sarikas and D. de Sousa Seixas,
Axial Symmetry Breaking inSelf-Induced Flavor Conversion of Supernova Neutrino Fluxes , Phys. Rev. Lett. (2013) 091101 [ ].[17] G. Mangano, A. Mirizzi and N. Saviano,
Damping the neutrino flavor pendulum bybreaking homogeneity , Phys. Rev.
D89 (2014) 073017 [ ].[18] B. Dasgupta and A. Mirizzi,
Temporal Instability Enables Neutrino FlavorConversions Deep Inside Supernovae , Phys. Rev. D (2015) 125030 [ ].[19] S. Abbar and H. Duan, Neutrino flavor instabilities in a time-dependent supernovamodel , Phys. Lett. B (2015) 43 [ ].[20] R. F. Sawyer,
Neutrino cloud instabilities just above the neutrino sphere of asupernova , Phys. Rev. Lett. (2016) 081101 [ ].[21] S. Chakraborty, R. S. Hansen, I. Izaguirre and G. Raffelt,
Self-induced neutrino flavorconversion without flavor mixing , JCAP (2016) 042 [ ].[22] B. Dasgupta, A. Mirizzi and M. Sen, Fast neutrino flavor conversions near thesupernova core with realistic flavor-dependent angular distributions , JCAP (2017)019 [ ].[23] B. Dasgupta and M. Sen, Fast Neutrino Flavor Conversion as Oscillations in aQuartic Potential , Phys. Rev. D (2018) 023017 [ ].[24] S. Abbar and H. Duan, Fast neutrino flavor conversion: roles of dense matter andspectrum crossing , Phys. Rev. D (2018) 043014 [ ].[25] B. Dasgupta, A. Mirizzi and M. Sen, Simple method of diagnosing fast flavorconversions of supernova neutrinos , Phys. Rev. D (2018) 103001 [ ].[26] F. Capozzi, B. Dasgupta, A. Mirizzi, M. Sen and G. Sigl, Collisional triggering of fastflavor conversions of supernova neutrinos , Phys. Rev. Lett. (2019) 091101[ ].[27] S. Abbar and M. C. Volpe,
On Fast Neutrino Flavor Conversion Modes in theNonlinear Regime , Phys. Lett. B (2019) 545 [ ].[28] F. Capozzi, G. Raffelt and T. Stirner,
Fast Neutrino Flavor Conversion: CollectiveMotion vs. Decoherence , JCAP (2019) 002 [ ].[29] J. D. Martin, C. Yi and H. Duan, Dynamic fast flavor oscillation waves in denseneutrino gases , Phys. Lett.
B800 (2020) 135088 [ ].[30] L. Johns, H. Nagakura, G. M. Fuller and A. Burrows,
Neutrino oscillations insupernovae: angular moments and fast instabilities , Phys. Rev. D (2020) 043009[ ].[31] F. Capozzi, M. Chakraborty, S. Chakraborty and M. Sen,
Fast flavor conversions insupernovae: the rise of mu-tau neutrinos , . – 22 –
32] L. Johns, H. Nagakura, G. M. Fuller and A. Burrows,
Fast oscillations, collisionlessrelaxation, and spurious evolution of supernova neutrino flavor , Phys. Rev. D (2020) 103017 [ ].[33] M.-R. Wu and I. Tamborra,
Fast neutrino conversions: Ubiquitous in compact binarymerger remnants , Phys. Rev. D (2017) 103007 [ ].[34] T. Morinaga, H. Nagakura, C. Kato and S. Yamada, Fast neutrino-flavor conversionin the preshock region of core-collapse supernovae , Phys. Rev. Res. (2020) 012046[ ].[35] M. Delfan Azari, S. Yamada, T. Morinaga, H. Nagakura, S. Furusawa, A. Haradaet al., Fast collective neutrino oscillations inside the neutrino sphere in core-collapsesupernovae , Phys. Rev. D (2020) 023018 [ ].[36] S. Abbar, H. Duan, K. Sumiyoshi, T. Takiwaki and M. C. Volpe,
Fast Neutrino FlavorConversion Modes in Multidimensional Core-collapse Supernova Models: the Role ofthe Asymmetric Neutrino Distributions , Phys. Rev. D (2020) 043016[ ].[37] R. Glas, H. T. Janka, F. Capozzi, M. Sen, B. Dasgupta, A. Mirizzi et al.,
FastNeutrino Flavor Instability in the Neutron-star Convection Layer of Three-dimensionalSupernova Models , Phys. Rev. D (2020) 063001 [ ].[38] M. George, M.-R. Wu, I. Tamborra, R. Ardevol-Pulpillo and H.-T. Janka,
Fastneutrino flavor conversion, ejecta properties, and nucleosynthesis in newly-formedhypermassive remnants of neutron-star mergers , .[39] S. Abbar, F. Capozzi, R. Glas, H. T. Janka and I. Tamborra, On the characteristics offast neutrino flavor instabilities in three-dimensional core-collapse supernova models , .[40] R. S. Hansen, S. Shalgar and I. Tamborra, Neutrino flavor mixing breaks isotropy inthe early universe , .[41] S. Bhattacharyya and B. Dasgupta, Late-time behavior of fast neutrino oscillations , Phys. Rev. D (2020) 063018 [ ].[42] S. Bhattacharyya and B. Dasgupta,
Fast Flavor Depolarization of SupernovaNeutrinos , .[43] G. Sigl and G. Raffelt, General kinetic description of relativistic mixed neutrinos , Nucl. Phys. B (1993) 423.[44] B. Dasgupta, T. Fischer, S. Horiuchi, M. Liebendorfer, A. Mirizzi, I. Sagert et al.,
Detecting the QCD phase transition in the next Galactic supernova neutrino burst , Phys. Rev. D (2010) 103005 [ ].[45] A. Banerjee, A. Dighe and G. Raffelt, Linearized flavor-stability analysis of denseneutrino streams , Phys. Rev. D (2011) 053013 [ ].[46] I. Izaguirre, G. Raffelt and I. Tamborra, Fast Pairwise Conversion of SupernovaNeutrinos: A Dispersion-Relation Approach , Phys. Rev. Lett. (2017) 021101[ ].[47] F. Capozzi, B. Dasgupta, E. Lisi, A. Marrone and A. Mirizzi,
Fast flavor conversions – 23 – f supernova neutrinos: Classifying instabilities via dispersion relations , Phys. Rev. D (2017) 043016 [ ].[48] T. Morinaga and S. Yamada, Linear stability analysis of collective neutrinooscillations without spurious modes , Phys. Rev. D (2018) 023024 [ ].[49] S. Airen, F. Capozzi, S. Chakraborty, B. Dasgupta, G. Raffelt and T. Stirner, Normal-mode Analysis for Collective Neutrino Oscillations , JCAP (2018) 019[ ].[50] C. Yi, L. Ma, J. D. Martin and H. Duan, Dispersion relation of the fast neutrinooscillation wave , Phys. Rev. D (2019) 063005 [ ].[51] M. Chakraborty and S. Chakraborty, Three flavor neutrino conversions in supernovae:slow \ & fast instabilities , JCAP (2020) 005 [ ].[52] C. Doering, R. S. Hansen and M. Lindner, Stability of three neutrino flavor conversionin supernovae , JCAP (2019) 003 [ ].].