Self-induced flavor conversion of supernova neutrinos on small scales

Self-induced flavor conversion of supernova (SN) neutrinos is a generic feature of neutrino-neutrino dispersion. The corresponding run-away modes in flavor space can spontaneously break the original symmetries of the neutrino flux and in particular can spontaneously produce small-scale features as shown in recent schematic studies. However, the unavoidable"multi-angle matter effect"shifts these small-scale instabilities into regions of matter and neutrino density which are not encountered on the way out from a SN. The traditional modes which are uniform on the largest scales are most prone for instabilities and thus provide the most sensitive test for the appearance of self-induced flavor conversion. As a by-product we clarify the relation between the time evolution of an expanding neutrino gas and the radial evolution of a stationary SN neutrino flux. Our results depend on several simplifying assumptions, notably stationarity of the solution, the absence of a"backward"neutrino flux caused by residual scattering, and global spherical symmetry of emission.


Introduction
Core-collapse supernovae (SNe) or neutron-star mergers are powerful neutrino sources and probably the only astrophysical phenomena where these elusive particles are dynamically important and crucial for nucleosynthesis [1,2]. The low energies of some tens of MeV imply that β reactions of the type ν e + n ↔ p + e − andν e + p ↔ n + e + are the dominant chargedcurrent processes. Heavy-lepton neutrinos ν µ ,ν µ , ν τ , andν τ , in this context often collectively referred to as ν x , interact only by neutral-current processes. Therefore, neutrino energy transfer and the emitted fluxes depend on flavor and one may think that flavor oscillations are an important ingredient. However, the large matter effect in this dense environment implies that eigenstates of propagation and those of interaction very nearly coincide [3,4].
In spite of large mixing angles, flavor oscillations are irrelevant except for MSW conversion when neutrinos pass the resonant density as they stream away [5,6]. Therefore, the neutrino signal from the next galactic SN may carry a detectable imprint of the yet unknown neutrino mass hierarchy [7][8][9][10][11]. This picture can fundamentally change when the refractive effect of neutrinos on each other is included [12][13][14]. The mean field representing background neutrinos can have flavor off-diagonal elements ("off-diagonal refractive index") due to flavor coherence caused by oscillations and can lead to strong flavor conversion effects . (We ignore additional effects that would arise from non-standard neutrino interactions [55], spin-flip effects caused by neutrino magnetic dipole moments [56][57][58], by refraction in inhomogeneous or anisotropic media [59][60][61][62], or the role of neutrino-antineutrino pair correlations [63][64][65][66].) Self-induced flavor conversion preserves the global flavor content of the ensemble, but re-shuffles it among momentum modes or between neutrinos and antineutrinos. The simplest example would be a gas of ν e andν e converting to ν µ andν µ , leaving the overall flavor content unchanged. The interacting modes of the neutrino field can be seen as a collection of coupled oscillators in flavor space. The eigenmodes of this interacting system include collective harmonic oscillations, but can also include run-away solutions (instabilities) which lead to self-induced flavor conversion [32]. Under which physical conditions will instabilities occur, how can we visualize them, and how will they affect the flavor composition of neutrinos propagating in the early universe or stream away from a SN core?
To study these questions, many simplifications were used and especially symmetry assumptions were made to reduce the dimensionality of the problem. However, symmetry assumptions suppress those unstable solutions which break the assumed symmetry. Therefore, when instabilities are the defining feature of the dynamics, symmetry assumptions about the solutions can lead to misleading conclusions because, even if the system was set up in a symmetric state, the interacting ensemble can break this symmetry spontaneously. This behavior is analogous to the hydrodynamical aspects of SN physics which cannot show convective overturn if the simulation is spherically symmetric, yet such 3D effects are now understood to be crucial for SN physics [67][68][69][70][71][72][73][74][75][76][77][78][79][80][81][82].
Our present concern is the question of "spatial spontaneous symmetry breaking" in selfinduced flavor conversion. In previous studies, the flavor content of neutrinos streaming away from a SN core was taken to remain uniform in the transverse directions. However, recent -1 -

JCAP01(2016)028
studies of simplified systems suggest that this symmetry can be spontaneously broken [51][52][53][54]. To avoid the complication of global spherical coordinates, it is sufficient to model the neutrino stream at some distance with plane-parallel geometry, i.e., we can use wave vectors k in the transverse plane to describe small-scale spatial variations. In this terminology, traditional studies only considered k = 0 (global spherical symmetry). In analogy to this k = 0 case it was found that for any k there is some range of effective neutrino densities where unstable solutions exist [51][52][53][54]. We usually express the neutrino density in terms of an effective neutrino-neutrino interaction energy µ = √ 2G F n νe (r) (R/r) 2 , where n νe (r) is the ν e density at distance r and R is some reference radius playing the role of the neutrino sphere. The parameter µ varies with r −4 because the neutrino density decreases as r −2 with distance. In this terminology, for any k there is a range µ min < µ < µ max where the system is unstable. For larger k (smaller spatial scales), the instability range shifts to larger µ, i.e., to regions closer to the SN core. This finding suggests that the neutrino stream is never stable because at any neutrino density there is some range of unstable k modes.
However, such conclusions may be premature as one also needs to include the refractive effect of matter which also tends to shift the instability to larger µ, a phenomenon termed "multi-angle matter suppression" of the instability [26,37]. We usually express the multiangle matter effect in terms of the parameter λ = √ 2G F n e (r) (R/r) 2 , where n e (r) is the electron density at distance r. One should study the instability region in the two-parameter space of effective matter and neutrino densities, λ and µ, which we call the "footprint of the instability." In figure 1 we show as an example the footprint of the MAA instability for a schematic SN model (MAA stands for "multi azimuth angle," i.e., one type of instability). We always define the instability region by the requirement that the growth rate κ > 10 −2 ω 0 , where ω 0 is a typical vacuum oscillation frequency.
Our schematic SN model consists of neutrinos and antineutrinos with a single energy corresponding to a vacuum oscillation frequency ω 0 = ∆m 2 /2E = 1 km −1 . We assume they are emitted isotropically at the neutrino-sphere radius R = 30 km. We consider two-flavor oscillations and, after subtracting the ν x flux, we are left with a ν e andν e flux such that there emerge twice as many ν e thanν e . Finally, we assume that the effective neutrino-neutrino interaction energy, µ = √ 2G F n νe (R/r) 2 , at the neutrino sphere r = R is µ 0 = 10 5 km −1 . Under these circumstances, run-away solutions for k = 0 exist within the footprint area shown in figure 1 as a blue-shaded region. We also show, as a solid red line, a possible density profile, corresponding to electron density as function of radius. The sudden density drop at r ≈ 200 km is the shock wave. In this example, the density profile does not intersect with the instability footprint for radii below the shock wave. On the other hand, for larger radii, collective flavor conversion would begin. We further show the footprint for inhomogeneities with assumed wave-number k = 10 2 and k = 10 3 , measured in units of the vacuum oscillation frequency ω 0 . Notice that we define k in "co-moving" coordinates along the radius, i.e., a fixed k represents a fixed angular scale, not a fixed length scale. Notice also that for the MAA instability shown here, which is relevant for normal mass ordering, non-zero k values lead to two instability regions. This phenomenon does not arise for the bimodal instability.
The full range of all k-values inevitably fills the entire region below the k = 0 mode (blue shading) in this plot, i.e., the entire gray-shaded region is unstable, whereas the region above the blue-shaded part remains unscathed. In other words, essentially the largest-scale mode with k = 0 is the "most dangerous" mode. If this mode is stable on the locus of the SN density profile in figure 1, the higher-k modes are stable as well. Of course, if any instability is encountered by the physical SN density profile, these instabilities will span a range of scales and create complicated flavor conversion patterns.
where R is the neutrino-sphere radius, and matter density λ = √ 2G F n e (R/r) 2 for the schematic SN model described in the text. Because µ ∝ r −4 , the horizontal axis is equivalent to the distance from the SN as indicated on the lower horizontal axis. We also show a representative schematic SN density profile where the sharp density drop marks the shock wave. We also show the instability footprint explicitly for co-moving wave numbers k = 10 2 and k = 10 3 in units of the vacuum oscillation frequency. Notice that for the same value of k there are two separate instability strips. The collection of all small-scale instabilities fill the gray-shaded region below the traditional k = 0 (blue shaded) instability region, whereas they leave the space above untouched.
The rest of our paper is devoted to substantiating this main point and to explain our exact assumptions. We stress that our simplifications may be too restrictive to provide a reliable proxy for a realistic SN. In particular, we assume stationary neutrino emission and that the solution is stationary as well, i.e., we assume that the evolution can be expressed as a function of distance from the surface alone. We also ignore the "halo flux" caused by residual scattering which can be a strong effect. Our study would not be applicable at all in regions of strong scattering, i.e., below the neutrino sphere. We assume that the original neutrino flux is homogeneous and isotropic in the transverse directions, i.e., global spherical symmetry of emission at the neutrino sphere. It has not yet been studied if this particular assumption has any strong impact on the stability question, i.e., if violations of such an ideal initial state substantially change the instability footprint, or if such disturbances would simply provide seeds for instabilities to grow. It is impossible to understand and study all effects at once, so here we only attempt to get a grasp of the differential impact of including spatial inhomogeneities in the form of self-induced small-scale flavor instabilities. All the other questions must be left for future studies.

Equations of motion
Beginning from the full equation of motion for the neutrino density matrices in flavor space, we develop step-by-step the simplified equations used in our linearized stability analysis. In particular, we formulate the stationary, spherical SN problem, where the flavor evolution is a function of radius, as an equivalent time-dependent 2D problem in the tangential plane. A fixed neutrino speed in the tangential plane corresponds to the traditional "single angle" treatment, whereas neutrino speeds taking on values between 0 and a maximum, determined by the distance from the neutrino sphere, corresponds to the traditional "multi angle" case.

Setting up the system
We describe the neutrino field in the usual way by 3×3 flavor matrices (t, r, E, v), where the diagonal elements are occupation numbers for the different flavors, whereas the off-diagonal elements contain correlations among different flavor states of equal momentum. We follow the convention where antineutrinos are described by negative energy E and the corresponding matrix includes a minus sign, i.e., it is a matrix of negative occupation numbers.
We always work in the free-streaming limit, ignoring neutrino collisions. In this case, neutrino propagation is described by the commutator equation [14,45] i where and H are functions of t, r, E, and v. The Hamiltonian matrix is where M 2 , the matrix of neutrino mass-squares, is what causes vacuum oscillations. The matrix of charged-lepton densities, N , provides the usual Wolfenstein matter effect. The integration dΓ is over the neutrino and antineutrino phase space. Because antineutrinos are denoted with negative energies, we have explicitly dΓ = +∞ −∞ dE E 2 dv /(2π) 3 and the velocity integration dv is over the unit sphere. Because the neutrino speed |v| = 1 we were able, for later convenience, to write the current-current velocity factor in the unusual form Studying this 7-dimensional problem requires significant simplifications. For neutrino oscillations in the early universe, one will usually assume initial conditions at some time t = 0 and then solve these equations as a function of time. To include spatial variations, one may Fourier transform these equations in space, replacing the spatial dependence on r by a wave-number dependence k, whereas v · ∇ r → iv · k and the r.h.s. becomes a convolution of Fourier modes [51]. One can then perform a linearized stability analysis for every mode k and identify when modes of different wave number are unstable and lead to self-induced flavor conversion [52]. One can also use this representation for numerical studies [51,53,54].
The other relatively simple case is inspired by neutrinos streaming from a supernova (SN) core. One assumes that, on the relevant time scales, the source is stationary and that the solution is stationary as well, so that ∂ t → 0. In addition, one assumes that neutrinos stream only away from the SN so that it makes sense to ask about the variation of the neutrino flavor content as a function of distance, assuming we are provided with boundary conditions at some radius R which we may call the neutrino sphere. Actually, this description can be a poor proxy for a real SN because the small "backward" flux caused by residual neutrino -4 -

JCAP01(2016)028
scattering in the outer SN layers, the "halo flux," can be surprisingly important for neutrinoneutrino refraction because of its broad angular range [40][41][42]. Here we will ignore this issue and use the simple picture of neutrinos streaming only outward.
We stress that this simplification is the main limitation of our study and its interpretation in the physical SN context. If (self-induced) instabilities exist on small spatial scales, they could even exist below the neutrino sphere where the picture of neutrinos streaming only in one direction would be very poor. The approach taken here to reduce the 7-dimensional problem to a manageable scope may then hide the crucial physics. Therefore, our case study leaves open important questions about a real SN.

Large-distance approximation
The main point of our study is to drop the assumption of spatial uniformity, i.e., we include variations transverse to the radial direction. However, we are not interested in an exact description of large-scale modes. At some distance from the SN, outward-streaming neutrinos cannot communicate with others which travel in some completely different direction as long as we only include neutrino-neutrino refraction and not, for example, lateral communication by hydrodynamical effects. If we are only interested in relatively small transverse scales, we may approximate a given spherical shell locally as a plane, allowing us to use Cartesian coordinates in the transverse direction rather than a global expansion in spherical harmonics.
We now denote with z the "radial" direction, and use bold-faced characters to denote vectors in the transverse plane, notably a for the coordinate vector in the transverse plane and β the transverse velocity vector. In the stationary limit, the equation of motion (EoM) becomes where v z = 1 − β 2 and and H depend on z, a, E, and β. If the neutrino-sphere radius is R, then at distance r from the SN the maximum neutrino transverse velocity is β max ≈ R/r 1. This latter "large distance approximation" is the very justification for using Cartesian coordinates in the transverse direction.
Therefore, it is self-consistent to expand the equations to order β 2 . (We need to go to quadratic order lest the neutrino-neutrino interaction term vanishes.) The β expansion is not necessary for our stability analysis, but avoiding it does not provide additional precision and performing it provides significant conceptual clarity. Numerical precision for a specific SN model is not our goal, and in this case we would have to avoid modeling the transverse direction as a flat space anyway, especially when considering regions that are not very far away from the SN core.
In equation (2.3) we multiply with 1/v z ≈ (1 + 1 2 β 2 ) and notice that the gradient term remains unchanged if we expand only to order β 2 so that where the Hamiltonian matrix is the old one times (1 + 1 2 β 2 ) or explicitly The flux factor under the integral in the second expression is also an expansion to O(β 2 ) in the form Multiplying this expression with (1 + 1 2 β 2 ) makes no difference because it is already O(β 2 ).

JCAP01(2016)028
As a next step, we re-label our variables and denote the radial direction z as time t. Moreover, we rescale the transverse velocities as β = vβ max where v is now a 2D vector obeying 0 ≤ |v| ≤ 1. Coordinate vectors in the transverse plane are also rescaled as x = aβ max , i.e., the new transverse coordinate vector x is "co-moving" in that it denotes a fixed angular scale relative to the SN. After these substitutions, the EoMs are with Of course, the neutrino phase-space integration dΓ is understood in the new variables.
Our stationary 3D problem has now become a time-dependent 2D problem. In the SN interpretation, the aspect ratio of the neutrino sphere shrinks with distance and correspondingly, β max shrinks. In other words, the physics is analogous to neutrino oscillations in the expanding universe. In the SN case, linear transverse scales grow as r/R, where r is the distance to the SN, i.e., our "Hubble parameter" is R −1 where R is the neutrino-sphere radius and the scale factor grows linearly with "time." The physical neutrino density decreases with inverse-distance squared and in addition, the factor β 2 max accounts for the decreasing value of the current-current factor in the neutrino interaction. Therefore, the effective neutrino number density decreases with (scale factor) −4 in the familiar way.

Single angle vs. multi angle
In the SN context, one often distinguishes between the single-angle and multi-angle cases, referring to the zenith angle of neutrino emission at the SN core. If all neutrinos were emitted with a fixed zenith angle, their transverse speeds would be |β| = β max and in our new variables, |v| = 1. In this case (1 + 1 2 β 2 max v 2 ) → (1 + 1 2 β 2 max ) is simply a small and negligible numerical correction to the vacuum oscillation frequencies and the matter effect. Also, we can revert to the traditional form of the flux factor 1 2 . Therefore, the SN single-angle case is equivalent, without restrictions, to a 2D neutrino gas evolving in time. Therefore, neutrino oscillations in an expanding space ("early universe") is exactly equivalent to the single-angle approximation of neutrinos streaming from a SN core with properly scaled effective neutrino and matter densities.
It has been recognized a long time ago that in the single-angle case, the ordinary matter effect has no strong impact on self-induced flavor conversion [19]. As usual, one can go to a rotating coordinate system in flavor space. In this new frame, the matrix of vacuum oscillation frequencies, M 2 /2E, has fast-oscillating off-diagonal elements and, in a time-averaged sense, it is diagonal in the weak-interaction basis. These fast-oscillating terms are what kick-starts the instabilities at the beginning of self-induced flavor conversion but are otherwise irrelevant. For a larger matter effect, more e-foldings of exponential growth of the instability are needed to "go nonlinear." In this sense, matter has a similar effect concerning the onset of the instability that would be caused by reducing the mixing angle. These effects concern the perturbations which cause the onset of instabilities, not the existence and properties of the unstable modes themselves.
We may ignore the small correction to the vacuum oscillation frequency provided by the factor (1 + 1 2 β 2 max v 2 ). We need to keep terms of order β 2 in the context of the matter -6 -

JCAP01(2016)028
and neutrino-neutrino term which in the interesting case are large and, after multiplication with β 2 max , still larger than the vacuum oscillation term. Therefore, we find where the first term symbolizes the time-averaged vacuum term in the fast co-rotating frame.
In the single-angle case, where v 2 = 1 for all modes, the remaining matter term can be rotated away as well.
The multi-angle SN case, in this representation, corresponds to a 2D neutrino gas with variable propagation speed 0 ≤ |v| ≤ 1, i.e., the velocity phase space is not just the surface of the 2D unit sphere (a circle with |v| = 1), but fills the entire 2D unit sphere (a disk with |v| ≤ 1). There is no counterpart to this effect in a "normal" neutrino gas. The early-universe analogy does not produce multi-angle effects, but we can include them without much ado by allowing the neutrino velocities to fill the 2D unit sphere.
If neutrinos are emitted "black-body like" from a spherical surface, from a distance this neutrino sphere looks like a disk of uniform surface brightness, in analogy to the solar disk in the sky. Therefore, this assumption corresponds to the neutrino transverse velocities filling the 2D unit sphere uniformly. In earlier papers of our group, we have used the variable u = v 2 with 0 ≤ u ≤ 1 as a co-moving transverse velocity coordinate, representing the neutrino zenith angle of emission. In terms of this variable, the black-body like case corresponds to the familiar top-hat u distribution on the interval 0 ≤ u ≤ 1.
Our overall set-up was inspired by that of Duan and Shalgar [52], except that they consider only one transverse dimension with single-angle emission at the SN. In other words, their system is equivalent to two colliding beams evolving in time and allowing spatial variations. A numerical study of this case, in both single and multi angle, was very recently performed by Mirizzi, Mangano and Saviano [53], going beyond the linearized case. Not unexpectedly, the outcome of the nonlinear evolution is found to be flavor decoherence. In our approach, multi-angle effects in this colliding-beam system can be easily included by using a 1D velocity distribution that fills the entire interval −1 ≤ v ≤ +1 and not just the two values v = ±1. A black-body like zenith-angle distribution corresponds to a uniform velocity distribution on this interval.

Two-flavor system
As a further simplification we limit our discussion to a two-flavor system consisting of ν e and some combination of ν µ and ν τ that we call ν x , following the usual convention in SN physics. We are having in mind oscillations driven by the atmospheric neutrino mass difference and by the small mixing angle Θ 13 . The vacuum oscillation frequency is where we have used ∆m 2 = 2.5 × 10 −3 eV 2 . Henceforth we will describe the neutrino energy spectrum by an ω spectrum instead, with negative ω describing antineutrinos. The matrix of vacuum oscillation frequencies, in the fast-rotating flavor basis, takes on the diagonal form

JCAP01(2016)028
where we have removed the part proportional to the unit matrix which drops out of commutator expressions. We do not include the fast-oscillating off-diagonal elements which is irrelevant for the stability analysis. The matter effect appears in a similar form, where again we have removed the piece proportional to the unit matrix. The parameter λ describing the multi-angle matter effect at distance r from the SN with neutrino-sphere radius R and using β max = R/r is where n e (r) is the net density of electrons minus positrons, ρ(r) the mass density, and Y e (r) the electron fraction per baryon, each at radius r. The matter density drops steeply outside the neutrino sphere and jumps downward by an order of magnitude at the shock-wave radius. Therefore, we need to consider λ values perhaps as large as some 10 7 km −1 all the way to vanishingly small values.
Turning to the neutrino-neutrino term, notice that the matrices play the role of occupation numbers and that the dΓ integration includes the entire phase space of occupied neutrino and antineutrino modes. Therefore, N ν = dΓ is a flavor matrix of net neutrino minus antineutrino number densities, in analogy to the corresponding charged-lepton matrix N . It is less obvious, however, how to best define an effective neutrino-neutrino interaction strength µ which plays an analogous role to λ. If we were to study a system that initially consists of equal number densities of ν e andν e , the matrix N ν vanishes, but later develops off-diagonal elements. Therefore, we rather use the number density of ν e without subtracting the antineutrinos and define where L νe is the ν e luminosity and E νe their average energy. More precisely, n νe is the ν e density at radius r that we would obtain in the absence of flavor conversions after emission at radius R. Previously we have sometimes normalized µ to nν e instead, or to the difference between theν e andν x densities. However, in our schematic studies we assume that initially we have only a gas consisting of ν e andν e , again obviating the need for these fine distinctions.
The exact definition of µ has no physical impact because it always appears as a product with the density matrices.
In previous papers [32,45], a further factor 1/2 was included in the definition of the multi-angle λ and µ. We have kept this factor explicitly in equation (2.8) both in the matter term and in the flux factor 1 2 (v − v ) 2 to maintain its traditional form. In this way, the equations can be directly applied to a traditional "early universe" system. To make contact with previous SN discussions, one can always absorb this factor in the definition of λ and µ.
As a next step, we project out the trace-free part of the density matrices and normalize them to account for the above normalization of the effective neutrino-neutrino interaction strength µ, (2.14)

JCAP01(2016)028
With these definitions, the two-flavor EoMs finally become where we have included a possible spatial dependence of the electron density in the form of λ x depending on location in the 2D space. The neutrino velocity domain of integration is determined by the dimensionality of the chosen problem and if multi-angle effects are to be considered.

Mass ordering
In a two-flavor system, one important parameter for matter effects in general and for selfinduced flavor conversion in particular is the mass ordering. In our context the question is if the dominant mass component of ν e is the heavier one (inverted ordering) or the lighter one (normal ordering). Traditionally "mass ordering" is also termed "mass hierarchy" and we denote the two cases as IH (inverted hierarchy) and NH (normal hierarchy). We are concerned with 1-3-mixing, the corresponding mixing angle is not large, and so it is clear what we mean with the "dominant mass component." Our equations are formulated such that they apply to IH, the traditional case where self-induced flavor conversion is important in the form of the bimodal instability. Of course, it has become clear that NH is actually the more interesting case. For NH, ∆m 2 is negative, but we prefer to consider ∆m 2 a positive parameter. Therefore, NH is achieved by including explicitly a minus sign on the r.h.s. of equation (2.10). This change of sign translates into a minus sign for ω in the first bracket in equation (2.16).
For flavor conversion, it is irrelevant if neutrinos oscillated "clockwise" or "counter clockwise" in flavor space, i.e., in equation (2.15) we may change i → −i or H → −H without changing physical results. However, the relative sign between ω and λ and µ is crucial. Therefore, switching the hierarchy is achieved by In our stability analysis we will consider the parameter range −∞ < µ < +∞ and −∞ < λ < +∞ as these are simply formal mathematical parameters. Physically both parameters being positive corresponds to IH, whereas the quadrant of both parameters being negative corresponds to NH.

Linearization
As a next step, we linearize the EoMs in the sense that the complex off-diagonal element of every G is supposed to be very small compared to its diagonal part. We write these matrices explicitly as where g is a real and G a complex number and all quantities carry indices (t, x, ω, v). To linear order in G we then find the EoMs Up to normalization, the "spectrum" g t,x,ω,v is essentially the phase-space density of all neutrinos. It is not affected by flavor conversion, but evolves by free-streaming if it is not homogeneous.

Homogeneous neutrino and electron densities
We are primarily interested in self-induced instabilities. Disturbances in the neutrino density and/or the electron density will certainly exist in a real SN and can play the role of seeds for growing modes. However, if these disturbances are small, it is unlikely that they will be responsible for instabilities themselves. Henceforth we will assume that the neutrino and electron densities do not depend on the transverse coordinate x, although the flavor content may well depend on x. Free streaming does not change the density if it is uniform. As a consequence, g t,x,ω,v does not depend on x or t and likewise, λ x does not depend on x.
With this assumption, g ω,v describes the initially prepared neutrino distribution, i.e., their density in the phase space spanned by ω and v. Inspecting the first integral in equation (2.19b), we may write the three independent terms as Here, represents the "asymmetry" between neutrinos and antineutrinos. The second term, 1 , represents a neutrino current which exists if their distribution is not isotropic and not symmetric between neutrinos and antineutrinos. Overall, the first term in square brackets becomes The last term is simply a constant and can be removed by changing the overall frequency of the rotating frame. Definingλ the term in square brackets effectively becomes ω We may also return to the notation used in our previous papers and define The linearized EoM then takes on the more familiar form

JCAP01(2016)028
This equation corresponds to equation (6) of reference [45]. Besides the streaming term (the gradient term on the l.h.s.) that we have now included to deal with self-induced inhomogeneities, we have also found the additional term µ 1 · v which is unavoidable in a non-isotropic system, irrespective of the question of homogeneity. This neutrino flux term is missing in reference [45]. The presence of this term modifies the eigenvalue equation for a non-isotropic system.

Spatial Fourier transform
We can now perform the spatial Fourier transform of our linearized EoM of equation (2.24). It simply amounts to replacing the spatial dependence on x of S by it dependence on the wave vector k and v · ∇ x → iv · k, leading to wherek = k − µ 1 . Therefore, the wave vector k has the same effect as a neutrino current.
Including an electron current would appear in a similar way. However, in the following studies of explicit cases we will focus on isotropic distributions and worry primarily about self-induced anisotropies, not about the modifications caused by initially prepared anisotropies.
In equation (2.25) we have spelled out the meaning of the phase-space integral dΓ = dω dv. Notice that the meaning of dΓ has changed in the course of changing variables that describe the neutrino modes. All phase-space factors and Jacobians have been absorbed in the definition of the effective neutrino-neutrino interaction strength µ as well as the normalization of the "spectrum" g ω,v . In particular, if we begin with an ensemble consisting of only ν e andν e and no ν x orν x , then our normalizations mean that ∞ 0 dω dv g ω,v = 1. In this latter integral, we have only included positive frequencies (neutrinos, no antineutrinos) so that this normalization coincides with our definition that µ is normalized to n νe .

Oscillation eigenmodes
In order to find unstable modes we seek solutions of our linearized EoM of the form S t,k,ω,v = Q Ω,k,ω,v e −iΩt , leading to an EoM in frequency space of the form Eigenvalues Ω = γ + iκ with a positive imaginary part represent unstable modes with the growth rate κ.

Monochromatic and isotropic neutrino distribution
In our explicit examples we will always consider monochromatic neutrinos with some fixed energy, implying a spectrum of two oscillation frequencies ω = ±ω 0 . Assuming that the energy and velocity distribution factorize, we may write the spectrum in the form The monochromatic energy spectrum is meaning that we have α antineutrinos (frequency ω = −ω 0 ) for every neutrino (ω = ω 0 ). The spectral asymmetry is = 1 − α.
We will consider isotropic velocity distributions which, in addition, are uniform, corresponding to blackbody-like angular emission in the SN context. In this case, where Γ v is the volume of the velocity phase space. The eigenvalue equation (2.26) finally simplifies to the form in which we will use it, We now consider systematically different cases of velocity distributions.

One-dimensional system
As a first case study we consider a 1D system, i.e., the toy model of "colliding beams" that has been used in the recent literature as a simple case where one can easily see the impact of spontaneous spatial symmetry breaking [46,[51][52][53]. We go beyond previous studies in that we include the multi-angle matter effect and study the "footprint" of the various instabilities in the two-dimensional parameter space −∞ < µ < +∞ and −∞ < λ < +∞. This schematic study already leads to the conclusion that essentially the largest-scale instabilities are "most dangerous" in the context of SN neutrino flavor conversion.

Eigenvalue equation
We begin with 1D systems, i.e., colliding beams of neutrinos and antineutrinos with different velocity distributions. The first case is what we call "single angle," a nomenclature which refers to the zenith-angle distribution of SN neutrinos. As we have explained, in our way of writing the equations, "single angle" means that the neutrino velocity distribution has |v| = 1. In our first 1D case this means we consider two colliding beams with v = ±1.
Matter effects can be rotated away. The eigenfunction Q Ω,k,ω,v now consists of four discrete components. We denote these four amplitudes with the complex numbers R for right-moving (v = +1) neutrinos (ω = +ω 0 ),R for right-moving antineutrinos, and analogous L andL for left movers. Our master equation (2.29) then reads corresponding to the equivalent result of Duan and Shalgar [52]. The eigenvalues Ω are found from equating the determinant of the matrix in square brackets with zero. This condition can be written in the form This expression depends only on µ 2 and therefore yields identical eigenvalues for positive and negative µ, i.e., for both neutrino mass hierarchies, as noted by Duan and Shalgar. It is also even under k → −k as it must because the system was set up isotropically, so the eigenvalues cannot depend on the orientation of k. For the homogeneous mode, k = 0, this eigenvalue equation simplifies considerably. We already know that we have the same solution for positive and negative µ, where the latter is the left-right symmetry breaking solution discovered in reference [46]. We here limit ourselves to µ > 0 and need to solve the quadratic equation It has the solutions Unstable solutions exist for 828. The maximum growth rate occurs at the interaction strength We show the growth rate normalized to its maximum in figure 2 as a function of µ/µ κmax . For this normalization, the unstable range is If we use α = 1 − and expand to lowest order in , this range is 2 /8 < µ/µ κmax < 2 − 2 /8. Therefore, even if is not very small ( = 1/2 in figure 2), the unstable range is close to its maximum range from 0 to 2.

Inhomogeneous modes (k > 0)
The quartic eigenvalue equation (3.2) is not easy to disentangle. However, for large k it simplifies and can be solved. We may guess that, for large k, the real part of Ω is approximately −k and, without loss of generality, we may go to a rotating frame such that Ω =Ω − k. Moreover, based on numerical studies, Duan and Shalgar [52] have conjectured that for large k, the unstable µ-range scales with √ k. This observation motivates us to write µ = m √ ω 0 k and, without loss of generality, the eigenvalue equation reads Now we can take the limit k → ∞ and approximate the second bracket as (1 − α)/2. The resulting quadratic equation is now easily solved and yields One can try the same exercise with the opposite rotating frame Ω =Ω + k and finds a similar-looking result, where however the argument of the square-root is always positive, i.e., there is no unstable mode with this property. From equation (3.2) one finds that the maximum growth rate κ max is the same as in the homogeneous case of equation (3.6), i.e., for large k the maximum growth rate does not depend on k and is the same as for k = 0. However, it now occurs at the interaction strength The unstable range is the same as given in equation (3.8) if we substitute µ/µ κmax with (µ/µ κmax ) 2 . The growth rate as a function of µ is the same as shown in figure 2 if we interpret the horizontal axis as (µ/µ κmax ) 2 with our new µ κmax .
For intermediate values of k, the maximum growth rate deviates slightly from the two extreme cases. It is somewhat surprising that the structure of the eigenvalue equation is such that the maximum possible growth rate depends only on the vacuum oscillation frequency ω 0 and α, but not on the potentially large scale k.

Eigenvalue equation
We may now study the multi-angle impact of matter, in the spirit of the SN system, in our 1D model by extending the velocity integration over the entire interval −1 ≤ v ≤ 1 instead of using only the modes v = ±1. One approach is to introduce n v discrete velocities on the interval 0 < v ≤ 1, i.e., 2n v modes on the interval −1 ≤ v ≤ 1, leading to a total of -14 -

JCAP01(2016)028
4n v discrete equations. One may then proceed to find numerically the eigenmodes. This approach is marred by the appearance of spurious instabilities and one may need a large number of modes to obtain physical results [43].
Therefore, we usually avoid discrete velocities and represent the eigenfunctions Q Ω,k,ω,v as continuous functions of their variables. Equation (2.29) reads for our specific case (3.12) Since our system was prepared "isotropic" in the sense of left-right symmetry, the instabilities cannot depend on the sign of k so that it is enough to consider 0 ≤ k < ∞. The r.h.s. of equation (3.12), as a function of v, has the form and v 2 are linearly independent functions on the interval −1 ≤ v ≤ +1, the l.h.s. must be of that form as well and we may use the ansatz for the eigenfunctions. Inserting this form on both sides yields This equation really consists of three linearly independent equations for the parts proportional to different powers of v, so we get three equations linear in A 0 , A 1 and A 2 . This set of linear equations is compactly written as where where we have dropped the prime from the integration variables. Notice that I n is odd under k → −k if n is odd, and it is even if n is even. Moreover, in the determinant of the matrix in square brackets in equation (3.15), every term involving I n with an odd power of n involves another factor I m with m odd. Therefore, the determinant is even under k → −k, in agreement with our earlier statement that without loss of generality we may assume k ≥ 0.

Homogeneous mode (k = 0) without matter effects (λ = 0)
As a first simple example we consider homogeneous solutions (k = 0) in the absence of matter effects (λ = 0). The latter assumption requires an exact cancellationλ = λ + µ = 0 between the matter effects caused by the background medium and by neutrinos themselves. We consider this case only for mathematical convenience without physical motivation. The velocity integrals vanish for odd powers of v. For even powers, and using the monochromatic frequency spectrum of equation (2.28), we find Before searching for solutions, we may diagonalize the 3×3 matrix. This leads to three independent equations of the form of equation (3.3), where we need to substitute Therefore, we get three instabilities: one for µ > 0, the usual bimodal instability (IH), and two negative-µ solutions (NH). The maximum growth rate is the same in every case as the one that was found in section 3.1.2 and was given in equation (3.6). The exact unstable µ-ranges have changed according to the µ scaling provided by equation (3.19). For our usual example with α = 1/2, the instability ranges for both hierarchies in the colliding-beam example were 0.68 < |µ/ω 0 | < 23.31. After v-integration they become explicitly We conclude that integrating over the velocity interval −1 ≤ v ≤ +1 modifies the unstable µranges, breaks the symmetry between normal and inverted hierarchy, and introduces another normal-hierarchy instability. Qualitatively, these results are analogous to the three types of instability discovered in the study of axial-symmetry breaking in the SN context [45]. The one inverted-hierarchy solution appearing in all cases is the bimodal instability and corresponds to the original flavor pendulum [21]. The first normal-hierarchy instability is what was termed the multi-azimuth angle (MAA) instability, although in our 1D system we have only two "azimuth angles," i.e., the two beam directions. The third instability, appearing in normal hierarchy for a much larger µ-range, is what was termed the "multi zenith angle" (MZA) instability. It requires, in the SN terminology, a nontrivial range of zenith angles, corresponding here to a non-trivial v-range, i.e., anything beyond the trivial v = ±1 velocity distribution.
Notice that our 1D MAA instability corresponds to an eigenfunction which is antisymmetric in v (it breaks the left-right symmetry) and corresponds to the middle entry −2/3 in the matrix of equation (3.18) which decouples from the remaining 2×2 block. The latter yields left-right symmetric solutions (even under v → −v). In this sense, it is the MZA instability which corresponds, for normal hierarchy, to the bimodal solution. It exists only for at least three velocity modes, and of course requires the presence of two vacuum oscillation frequencies, here always chosen as ω = ±ω 0 , i.e., we need at least a total of six modes, making a simple visual interpretation more difficult.
The growth rates of all three intabilities as functions of µ are shown in figure 3 as blue, green, and orange lines, all of them having the same maximum. In the top panel, we overlay the two instability curves for the original colliding-beam example, where v = ±1. Therefore, this panel directly illustrates the effect of the velocity integration ("multi-angle effects" in  SN terminology) in that the velocity integration takes us from the two black curves to the three colored ones.
We have also studied the equations using a set of discrete velocities, where the v = ±1 case is the simplest example with n v = 1 bins. (We count the number of bins in the range 0 < v ≤ 1, i.e., there are equally many bins for negative velocities, and the total number doubles for our two frequencies ω = ±ω 0 .) Adding the intermediate values v = ±1/2 takes us to n v = 2, shown in the second panel of figure 3. It reveals that the symmetry between the hierarchies (µ → −µ symmetry) is broken as soon as the velocity range is non-trivial and that there are two normal-hierarchy solutions. Increasing n v eventually approximates a uniform v distribution. A fairly small number of velocity bins is enough to achieve good agreement. We will see shortly that, including non-zero k and/orλ, changes the picture because spurious instabilities appear.

Inhomogeneous modes (k > 0) without matter effects (λ = 0)
Non-vanishing matter effects (λ = 0) and non-vanishing inhomogeneities (k = 0) modify the eigenvalue equation in similar ways: the range of effective oscillation frequencies given by 1 2λ v 2 + kv + ω increases considerably if k/ω 0 1 and/orλ/ω 0 1. Roughly one would suspect that significant collective phenomena require a neutrino-neutrino coupling exceeding this range of frequencies, i.e., a µ range exceeding something like the rms spread of this range. In this sense, one would expect that the µ-range of unstable solutions would be shifted roughly linearly withλ and/or k.
We first test this picture withλ = 0 and k > 0. The v-integrals in equation (3.16) can be performed analytically (appendix A) and the ω-integration amounts to summing over two terms with ω = ±ω 0 . However, finding the eigenvalues of equation (3.15) requires numerical tools. We use Mathematica and show the result in figure 4 for 0 ≤ k/ω 0 ≤ 100 as indicated above the curves. In the left panels, we only show the MZA mode. For relatively small k, the function shifts left and deforms somewhat, whereas for larger k values, it shifts left nearly linearly with k. We have checked that this linear behavior obtains numerically even for very large k values -we have tested values up to 10 7 . In other words, for large k, the instability curves are very similar as a function of µ/k, although they narrow somewhat. In contrast to the earlier two-beam example, there does not seem to be quite a universal function for large k. Moreover, in the earlier case, the scaling was with µ/ √ ω 0 k, i.e., the nontrivial v range has qualitatively changed the results with regard to the k-scaling.

JCAP01(2016)028
In the right panels of figure 4, we show analogous results for the MAA mode (µ < 0) and the bimodal mode (µ > 0), which show analogous behavior.
One can study the same case by solving the eigenvalue equation in terms of velocity bins, in analogy to what one would do if one were to solve the EoMs numerically instead of performing only a linear stability analysis. For k = 25 ω 0 we show the growth rates for all solutions in figure 5 for different choices for the number n v of velocity bins. (We recall that in our convention, the total number of positive and negative v bins is 2n v .) The large number of spurious modes is a conspicuous feature of these plots, although for sufficiently large n v , the physical modes stick out.
With non-vanishing k and/orλ, the functional form of the eigenfunctions in equation (3.13) no longer factorizes as a function of v and one of ω. It is conceivable that spurious modes can be avoided, or their impact reduced, if one were to find a better way of discretizing the neutrino modes than by simple bins in velocity and frequency [44].
It is noteworthy that for all modes, spurious or physical, the growth rates are of order ω 0 , i.e., they do not inherit a larger frequency scale from k or, in later cases, fromλ. Even for huge values of k andλ, this conclusion does not change and agrees with our explicit result in the two-beam example.

Including matter (λ = 0)
Including matter in our "multi-zenith-angle" case has the effect of introducing bothλ and k in the denominator of the integrals of equation (3.16). They can still be done analytically (appendix A), but lead to transcendental functions. Of course, numerically one can find the eigenvalues without much problem. The parameterλ, like k, has the effect of broadening the effective range of oscillation frequencies and of shifting the unstable collective modes to larger values of |µ|. We study the homogeneous and inhomogeneous cases separately.

Homogeneous mode (k = 0)
For the homogeneous mode (k = 0), the eigenvalue equation (3.15) simplifies considerably because in this case I 1 = I 3 = 0 and we are left with We need to solve the two equations (I 2 − 1) 2 = I 0 I 4 and I 2 = −1/2 . As an overview, we show in figure 6 a contour plot of the growth rate κ in the two-dimensional parameter space of the interaction strength µ and the effective matter densityλ = λ + µ. As explained earlier, the first quadrant (µ > 0 andλ > 0) corresponds physically to inverted mass ordering (IH), whereas the third quadrant (µ < 0 andλ < 0) corresponds to normal mass ordering (NH). The other quadrants would be relevant, for example, for a background medium of antimatter. Mathematically, µ andλ are simply parameters which we leave unconstrained by physical considerations. As usual, we use α = 1/2 and therefore = 1 − α = 1/2 so that matter-free space (λ = 0) corresponds to the lineλ = µ/2. Forλ = 0, the growth rates as a function of µ were shown in figure 3 in the form of the colored curves. The effect of increasing |λ| is to shift the unstable regions to larger  . Growth rates for all modes that appear when we consider n v positive and n v negative-v bins, using k = 25 ω 0 ,λ = 0, and α = 1/2. For growing n v , more and more spurious modes appear, but their growth rates become smaller and the physical modes begin to stick out.

Inhomogeneous mode (k > 0)
We next determine numerically the same instability footprints for non-vanishing wave numbers k. We already know the impact of nonzero k forλ = 0, i.e., the instability is shifted to larger µ values. We do not expect a big difference forλ k relative to the k = 0 case. These expectations are borne out by our results shown in figure 7. Considering first the simpler µ > 0 half of the plot, we show the instability footprints for k = 0, 10 2 , 10 3 and 10 4 as indicated in the plot. We can now easily diagnose the impact of small-scale instabilities: essentially they fill in the entire space between the k = 0 footprints between the quadrant with positive and negativeλ, whereas the region above the k = 0 footprint in the upper quadrant, and the space below in the lower quadrant remains stable. The only small caveat is that in the upper quadrant, for k ∼λ, there are "noses" of the footprint sticking into the previously stable region. So there is a narrow sliver of parameters above the k = 0 footprint, the envelope of the noses, which becomes unstable due to small-scale instabilities.
In the left half of the plot (µ < 0) the situation is somewhat more complicated because of the presence of two instabilities. For large k, the footprint actually connects asymptotically to the k = 0 instabilities in a crossed-over way which we illustrate only by the k = 10 4 case. In other words, for large k, the MAA and MZA instabilites strongly mix with each other.

Two-dimensional system
We now turn to a 2D system, corresponding to the SN example where neutrinos propagate within the "expanding transverse sheet" moving outward. The radial motion is parameterized by our "time" variable, whereas "space" is represented by two transverse directions. This example corresponds, with properly scaled variables µ and λ, to the usual treatment of selfinduced flavor conversion in SNe, except that now we can include small-scale instabilities with nonvanishing wavenumber k.

Eigenvalue equation
We begin with the "single zenith angle case," meaning that the rescaled neutrino speed within the transverse sheet is |v| = 1 and the matter effect can be rotated away. Our velocity phase space is the unit circle, described by an angle variable ϕ which we can measure relative to k. As the system is initially prepared axially symmetric, all vectors k are equivalent -the eigenvalues depend only on k = |k|. The eigenvalue equation (2.29) therefore reads where c ϕ = cos ϕ and s ϕ = sin ϕ. Proceeding as in our earlier cases, we notice that the r.h.s. of equation (4.1) has the form A 1 + A c cos ϕ + A s sin ϕ, i.e., a superposition of three linearly independent functions on the interval −π ≤ ϕ ≤ +π. Therefore, the l.h.s. must be of that form as well and we may use the eigenfunction ansatz We may insert this form on both sides, leading to three linearly independent equations, corresponding to the coefficients of the three functions 1, cos ϕ and sin ϕ. We may write the three equations in compact form where

JCAP01(2016)028
Here, f 1 (ϕ) = 1, f c (ϕ) = cos ϕ, f cc (ϕ) = cos 2 ϕ, and f ss (ϕ) = sin 2 ϕ. We have used that such integrals vanish if they involve a single power of sin ϕ because this is anti-symmetric on the integration interval, explaining the zeroes in the matrix in equation (4.3). We also note that I ss = I 1 − I cc , so we need only three different integrals.

Homogeneous mode (k = 0)
We first consider homogeneous solutions (k = 0) where the angle integrals in equation (4.4) can be performed explicitly. Using the monochromatic frequency spectrum of equation (2.28), the eigenvalue equation becomes These are three independent quadratic equations of the now-familiar form of equation (3.3). The first line corresponds to the usual bimodal solution, the second and third line to two degenerate multi-azimuth-angle (MAA) solutions which are unstable for negative µ (normal hierarchy). For these modes, the instability range is a factor of 2 larger.

Inhomogeneous modes (k > 0)
For the inhomogeneous modes (k > 0), the required integrals entering the eigenvalue equation are of the form With our monochromatic spectrum equation (2.28) we arrive at We define the auxiliary function of a complex argument w which, for complex numbers, is in general not equal to √ w 2 − 1. We then find These expressions allow us to write the eigenvalue equations explicitly, involving only polynomials and square-root expressions. We now have three non-degenerate solutions, in contrast to the original 1D example of Duan and Shalgar [52], one for positive µ (inverted hierarchy) and two for negative µ. We show contour plots of the growth rate κ as a function of µ and k for our usual example α = 1/2 in figure 8. The two quasi-symmetric regions correspond to the two solutions which correspond to those of the 1D case, i.e., the left-right symmetric and anti-symmetric cases also shown in reference [52] in a similar plot. In addition we see a third solution which is genuinely a result of the spatial 2D geometry with non-vanishing wave-vector k. The This quartic equation provides the large-k solution, which is symmetric for both hierarchies, i.e., symmetric under m → −m. Because the second-largest power of k is only a power 1/4 smaller than the largest, the asymptotic solution requires extremely large k-values. Notice that in the 1D beam example, the unstable range scaled with k 1/2 , while here it is k 3/4 .

Eigenvalue equation
We finally turn to the "multi zenith angle" case of transverse velocities within the full 2D disk described by |v| ≤ 1, meaning that multi-angle matter effects are now included. We describe the velocity phase space by the speed v = |v| and an angle variable ϕ which we measure relative to k as before. Noting that ( where c ϕ = cos ϕ and s ϕ = sin ϕ. The r.h.s. as a function of v and ϕ is A 0 + A 2 v 2 + A c v c ϕ + A s v s ϕ , so the eigenfunctions are of the form (4.12) As usual, we insert this form on both sides, leading to four linearly independent equations, corresponding to the coefficients of the four functions 1, v 2 , v c ϕ and v s ϕ . We may write the equations in compact form where we have used that dϕ vanishes for integrals involving an odd power of sin ϕ . The integrals are now (4.14) Here, f 1 (ϕ) = 1, f c (ϕ) = cos ϕ, f cc (ϕ) = cos 2 ϕ, and f ss (ϕ) = sin 2 ϕ. We also note that I ss 3 = I 1 3 − I cc 3 .

Homogeneous mode (k = 0) with matter (λ > 0)
In the homogeneous case (k = 0), the function cos ϕ in the denominator disappears, all angle integrals can be performed analytically, and I c n = 0. Therefore, equation (4.13) becomes where the integral expressions after performing the dϕ integration are

JCAP01(2016)028
The velocity integrals are explicitly The previous ss and cc blocks in equation (4.13) have now decoupled from the rest and are degenerate, leading to the MAA instability. The remaining 2×2 block provides the bimodal and MZA instability. In other words, we now need to solve (I 3 − 1) 2 = I 1 I 5 and I 3 = −1 , (4.18) in analogy to reference [45] with a slightly different notation. This entire development is very similar to the 1D case.
In the limitλ → 0, K 1 , K 3 and K 5 approach 1/2, 1/4 and 1/6 times 1/(ω − Ω), which one can also find by settingλ = 0 before doing the v-integrations. The matrix can be diagonalized and we find three independent quadratic equations of the form of equation (3.3) where we need to substitute µ → −µ/4 and µ → µ (3 ± 2 √ 3)/12. Following the same steps as in section 3.2.2, the instability ranges for our usual example α = 1/2 are found to be These results are numerically similar to equation (3.20) for the corresponding 1D-case. In analogy to the butterfly diagram of the 1D case (figure 6) we show in figure 9 a contour plot of the growth rate κ in the two-dimensional parameter space of the interaction strength µ and the effective matter densityλ = λ + µ. The result looks qualitatively similar to the 1D case. Again, for µ > 0 (inverted mass ordering), we obtain the bimodal instability while for µ < 0 (normal ordering) we find the MZA and MAA instabilities.

Inhomogeneous mode (k > 0) without matter (λ = 0)
Next we consider the relatively simple case of k > 0 without matter. We may write the integrals of equation (4.14) in the form and we have introduced w = (ω − Ω)/k. We find explicitly Notice that K ss 3 + K cc 3 = K 1 3 . With the help of these analytic integrals it is relatively easy to solve the eigenvalue equation numerically. We show a contour plot of the growth rate κ in the µ-k-plane in Growth Rate κ Figure 10. Growth rates for the 2D case with 0 < |v| = 1 ("multi zenith angle"), using α = 1/2. This figure is analogous to the single-angle case shown in figure 8, but now we have three instabilities for µ < 0 (normal mass ordering) and the usual bimodal one for µ > 0 (inverted ordering). For all four instabilities, the unstable µ range scales linearly with k as discussed in the text. figure 10. The 3×3 block in equation (4.13) provides three different solutions, i.e., one for µ > 0 (the usual bimodal solution in inverted mass ordering) and two solutions for µ < 0 (normal ordering). The 1×1 block provides a further solution for µ < 0. Figure 10 corresponds to figure 8 in the single-angle case. In comparison, we have one more instability, now a total of four, of which three are for µ < 0 (normal mass ordering). One may also extract the large-k asymptotic behavior (appendix D). For the 3×3 block in equation (4.13) one finds that the system is unstable for i.e., the unstable µ region scales linearly with k, in contrast to the corresponsing 1D case. The values of the coefficients a i and m max,i are given in appendix D.
For the 1×1 block in equation (4.13) one finds an instability on a very narrow strip around µ = −6 k. However, the maximum growth rate decreases with k −1/2 so that in the limit k → ∞ this instability disappears and we are left with those arising from the 3×3 block.

Inhomogeneous mode (k > 0) with matter (λ > 0)
As a grand finale, we now turn to the most general 2D case with matter (λ > 0) and inhomogeneities (k > 0). We need to find the zeroes of the determinant in equation (4.14) and write the integrals in the form where q = k/λ and w = (ω −Ω)/λ. These integrals can be performed analytically; we provide our results in appendix A. We next determine numerically the instability footprints for non-vanishing wave numbers k and show the result in figure 11, which looks qualitatively similar to the corresponding -29 -

JCAP01(2016)028
Effective Matter Density 1D case that was shown in figure 7. For the simpler µ > 0 half of the plot, we show the instability footprints for k = 0, 10 2 and 10 3 as indicated in the plot. These k > 0 footprints fill the space between the k = 0 footprint and the horizontal axis. In addition, in the upper panel, there are small "noses" of the k > 0 footprints which slightly extend in the space above the k = 0 footprint, but this is a very small effect. For µ < 0 the situation is more complicated because there are three instabilities. As in the 1D case, for large k the footprint connects asymptotically to the k = 0 instabilities in a crossed-over way, which we illustrate for k = 10 3 . The main novelty of the 2D case is the appearance of another instability, which merges with one of the others whenλ k. In other words, one of the instabilities somewhat splits into two unstable ranges. We also recall that for very large k, the maximum growth rate of one of them decreases and vanishes for -30 -

JCAP01(2016)028
k → ∞, in which case we are back to a total of three instabilities. In the unphysical third quadrant, we notice somewhat pronounced "noses" of the k > 0 footprints.
In all cases the main message is the same as in the 1D case: the small-scale instabilities fill the space between the k = 0 instability and the horizontal axis, whereas the space between the k = 0 instability and the vertical axis remains stable.

Conclusions
Several recent papers [46,[51][52][53] have studied the phenomenon of spatial spontaneous symmetry breaking in the "colliding beam" model of neutrinos interacting with each other refractively. One important finding was that for any neutrino density (or in our nomenclature for any value of the neutrino-neutrino interaction energy µ) there is some range of spatial wave vectors k where the system is unstable with regard to self-induced flavor conversion. As a consequence, it seemed that an interacting neutrino gas would never be stable for any conditions, with potentially far-reaching consequences for SN physics.
We have studied similar models, but including the multi-angle matter effect. We concur with the previous results in that smaller-scale modes are unstable at larger values of µ for a given matter density. This means that on a "footprint plot" such as figure 11, modes with k > 0 fill the space between the traditional footprint for k = 0 and the horizontal axis, but not the space towards the vertical axis. If we show the instability footprint in a plot like figure 1, adapted to more physical SN parameters, the large-k modes extend the instability region in the direction of larger neutrino density for fixed matter density.
Therefore, if the instability footprint of the traditional homogeneous (or rather spherically symmetric) mode does not intersect with the SN density profile, the large-k modes are safe as well. In this sense, the traditional large-scale mode remains the most sensitive stability probe. On the other hand, if the physical SN profile of density and neutrino fluxes intersects any instability region, instabilities on a large range of spatial scales will occur and one would not expect any simple outcome of the flavor conversion process.
Our analysis is based on a linearized stability analysis of a model which we have developed in section 2. We have formulated the problem is such a way that it includes, as special cases, the "colliding beam" examples of the previous literature, allows the inclusion of multi-angle effects by a simple modification of phase-space integration, and formulates the SN case as a 2D system evolving in time, i.e., the non-trivial dynamical evolution is in the 2D expanding sheet of neutrinos as a function of SN distance. All cases of the previous literature are covered in a single simple formulation.
In particular, our homogeneous, multi-angle, 2D system corresponds to the usual scenario described in the previous literature [45,82], where both the zenith and azimuthal multi-angle effects are included. Note that the usual single zenith angle description of the early days of the collective oscillation discussions [22] cover only a small part of the parameter space, for example the bimodal instabilities shown in our figure 8. In addition, for the first time we have have considered a scenario where inhomogeneous modes are included in the description of SN neutrino evolution. We have found that the homogeneous mode remains the dominant source of instability.
Still, in many ways our investigation is only a mathematical case study that may or may not apply to a realistic SN. Our main simplification is that we assume the flavor content of the SN neutrino stream to be stationary and to evolve only as a function of distance. Moreover, we assume a uniform boundary condition at the neutrino sphere, i.e., global spherical symmetry JCAP01(2016)028 of neutrino emission. In other words, our case study still contains substantial and nontrivial simplifying assumptions to reduce the complexity of the full problem. Future work will have to go beyond some of these simplifications to develop a more realistic understanding of what really happens to neutrino flavor in the dense SN environment.

A Analytic integrals
When searching for the complex eigenvalues Ω we need various integrals that can be found easily with Wolfram's Mathematica. There can be issues about the validity of the analytic expressions in the complex plane, so we here give the integrals explicitly.
In the 1D case, for a non-vanishing k and in the absence of matter effects (λ = 0), we need integrals of the form where w is a complex number. We first define two auxiliary functions A(w) = 2 + w i π sign Im(w) − 2 arctanh(w) . (A.2b) The required integrals are found to be The actual argument will be of the form w = (ω − Ω)/k. For non-vanishing matter effects (λ = 0) and non-vanishing k, we need integrals of the form where w is a complex number and p is real. Again we define two auxiliary functions The required integrals are found to be g 0 (p, w) = 2 B(p, w) , (A.6a) The actual arguments are going to be p = 2k/λ and w = 2(ω − Ω)/λ. In the 2D case to solve equation (4.14) we need the integrals defined in equation (4.21), i.e., integrals of the form where f 1 (ϕ) = 1, f c (ϕ) = cos ϕ, f cc (ϕ) = cos 2 ϕ, and f ss (ϕ) = sin 2 ϕ. These integrals can be found analytically with the help of Mathematica. We first define two auxiliary functions Our desired integrals are then found to be where S w = −i sign(Im w).

B Frequently encountered eigenvalue equations
Based on our "monochromatic" neutrino spectrum equation (2.28) with vacuum oscillation frequencies ±ω 0 and α the number of antineutrinos relative to neutrinos, we constantly encounter eigenvalue equations of the form

JCAP01(2016)028
whereμ is an interaction energy. In the simplest case, F (x) = x, this is the traditional eigenvalue equation for the flavor pendulum if we notice that ourμ = µ(1 − α)/2, whereas µ is the traditional interaction energy. We also encounter F (x) = √ x and F (x) = log(x). To study these equations, we transform them to dimensionless variables by the substitutions These substitutions allow us to easily take the limit of a "symmetric" neutrino distribution with = 1 − α → 0.
It has a nonvanishing imaginary part in the range We denote with K 1 (α, m) the function which is the positive imaginary part of w. In the present case, it simply is As a function of m, it is a semi-circle with center at m = 1 and radius 2 √ α/(1 + α), i.e., To compare with section 3.1.2, notice thatμ = µ(1 − α)/2 = mω 0 (1 + α)/(1 − α). Using this relationship between µ and m as well as the one between Ω and w reproduces the previous results. In particular, for α → 1 we find the physical growth rate κ = (2µ − ω 0 )ω 0 which grows without limit for µ → ∞. It is interesting that in terms of the scaled variable m one obtains, as a limiting results for α → 1, a semi-circle for K(1, m) as a function of m. So in principle one could study all of our problems in the α → 1 limit, and yet obtain representative results for α < 1, i.e., one could essentially eliminate the annoying parameter α and still obtain meaningful results for the asymmetric case.
It can be transformed to a quartic equation, but the results are too cumbersome to deal with and not particularly informative. Again we can expand this equation in powers of 1 − α and find, for the limit α → 1, the cubic equation 2w 3 + m(1 + 2w) 2 = 0. The imaginary part of the solution is This function is shown in the second panel of figure 12 and looks almost like a semi-circle, but is not quite one. It is is nonzero for 0 < m < 27/16 = 1.6875 and takes on its maximum show K 1/2 (α, m) also for several other value of α. The curves always begin at m = 0, i.e., there is no lower threshold in this case. We finally turn to the logarithmic case, F (x) = log(x), natural logarithm always understood, for which equation (B.3) becomes There is no general analytic solution, but again we can consider the α → 1 expansion where we find 1 = w[1 + log(−w/  figure 12. Again, it looks deceivingly like a semi-circle, but is not exactly one. Its maximum occurs at m = 1.76684 and K max log = 0.724611. We usually consider α = 1/2 as our main example. In this case, the eigenvalue equation can be solved analytically with the solution The imaginary part is nonzero for 0 < m < 4e/3, has its maximum at m = 2e/3, and takes on a maximum value of 2/3.

C Asymptotic solutions for 1D andλ → ∞
We can derive analytic asymptotic solutions for the 1D case with matter, i.e., the large-λ continuation of the contour plot figure 6. We begin with the bimodal and MZA instability for λ > 0 and consider the first block of the eigenvalue equation (3.22). It is of the form 1+C 1 µ+ C 2 µ 2 = 0, where the coefficients C 1 and C 2 depend on α, ω 0 , Ω andλ. We have evaluated the integrals according to the explicit transcendental functions given in equation (A.6). Inspired by numerical solutions, we assume that both the real and imaginary parts of the solutions Ω remain of order ω 0 and do not become large asλ → ∞, an assumption that later bears out to be consistent with the solutions. Therefore, we may expand C 1 and C 2 in powers ofλ −1 and find that the dominant terms are C 1 ∝λ −1 and C 2 ∝λ −3/2 . In terms of a dimensionless interaction strengthμ of order unity we write where the exact coefficient was chosen for later convenience. The lowest-order term in C 2 µ 2 no longer depends onλ, whereas the lowest-order term in C 1 µ ∝λ −1/4 and slowly becomes small asλ → ∞. To lowest order inλ −1 , the eigenvalue equation is found to be This is an example of the type of equations that we always encounter in this context and which are discussed in appendix B. The asymptotic solution derives from the term quadratic in µ and thus remains unchanged under µ → −µ, i.e., it applies to both hierarchies. We show the asymptotic solution as a blue curve in figure 13 on a linear and logarithmic scale. We also show the growth rates forλ = 10 2 , 10 4 and 10 6 where the solution is not symmetric under µ → −µ because the linear term in µ kicks in. We have already noted that one needs very largeλ values to obtain the asymptotic solution because the second-largest term only scales withλ −1/4 relative to the dominant term. The asymptotic behavior is achieved for much smallerλ values if µ > 0. The growth rate vanishes completely above a certain |µ| value, but obtains nonzero values otherwise, i.e., there is no lowerμ threshold. However, forμ 0.5, the growth rate is a steep power-law ofμ and can be taken to be effectively zero.
For our usual example α = 1/2 we find numerically that the maximum growth rate occurs forμ = 1.494. Therefore, we find that µ = ± 4.911 ω gives us the locus of the maximum growth rate in the µ-λ plane for the bimodal and MAA solutions. The maximum value ofμ before the growth rate becomes zero is 1.7724. On the small-μ side, the growth rate drops below κ < 1/100, our usual criterion, atμ = 0.3478. Therefore, the footprint of the instability is the region between the lines µ = 1.143 λ 3/4 and 5.826 λ 3/4 , where both µ and λ are given in units of the vacuum oscillation frequency ω 0 . This footprint is shown in the first quadrant (upper right) of figure 15. The corresponding footprint in the second quadrant (upper left) is also shown.
We next turn to the MZA solution which exists only in inverted hierarchy (µ < 0) and we consider the second block in equation (3.22). If we use µ = −λ/[2(1−α)] the leading terms cancel, leaving us with a leading term of orderλ −1/2 . To obtain the lowest-order equation, we introduce another dimensionless parameterμ and write where, of course, the detailed coefficients in the second term are chosen for later convenience. One then finds a quadratic equation with solutions Notice that these solutions require 0 ≤μ ≤ 1 and we have always assumed 0 ≤ α ≤ 1. The imaginary part, as a function ofμ 2 , has the familiar semi-circular shape. In figure 14 we show it as a function ofμ (blue curve) and we also show the full solution forλ = 10 3 and 10 2 . The asymptotic solution is quite good for relatively smallλ values. In contrast to the other solutions, on a logarithmic scale the unstable range becomes very narrow asλ → ∞.
The maximum growth rate obtains forμ = 1/ √ 2. Therefore, for our usual example α = 1/2 we find that for the MZA solution, gives us the locus of the maximum growth rate in the µ-λ plane. The growth rate becomes exactly zero forμ ≤ 0 andμ ≥ 1, so the footprint (see upper-left quadrant in figure 14) is delimited by the curves µ = −λ and µ = −λ − π 3λω 0 . The width of the footprint scales with √λ , i.e., on a logarithmic scale it becomes very narrow for largeλ. Forλ < 0, the above approach does not lead to unstable solutions. Numerically we observe that forλ → −∞, the real part of the solutions approaches Re(Ω) →λ/2, i.e., a large negative number. Therefore, to be able to expand the equation, we express Ω =λ/2 + w ω 0 and seek self-consistent solutions with the dimensionless eigenvalue w of order unity. After expansion forλ → −∞, the asymptotic eigenvalue equations are where e is Euler's number and as always the logarithm is with base e. This equation is one example for the type discussed in appendix B. For our usual example α = 1/2, we can solve this equation analytically with the explicit result

JCAP01(2016)028
It has a nonzero imaginary part for −∞ < a < log(8) = 2.0794, although it becomes exponentially small for a −1. The maximum imaginary part obtains for a = log(4) = 1.3863 and the maximum is 2. Therefore, the maximum growth rate obtains for Therefore, we have altogether three solutions, corresponding to the three instabilities, with maximum growth rates on the locus in the µ-λ plane given by where the limiting behavior is understood for A → ∞. Becauseλ → −∞, the first solution corresponds to positive µ and thus to the bimodal solution, the second and third solutions are the MAA and MZA instabilities, respectively. To draw the footprints in the lower quadrants of figure 15, we notice that κ = 0 for a > log(8) and on the other side κ < 1/100 for a < log(4 − √ 39999/50) = −9.90348. Therefore, the asymptotic footprints are limited by a 1 = log(8) and a 2 = log(4 − √ 39999/50) (C. 13) from which the limiting curves are extracted by solving equation (C.8) forμ. Once more we note that two of the footprints are "wide" and nearly symmetric between µ → −µ, whereas the third instability has a very narrow footprint.
D Asymptotic solutions for 2D withλ = 0 and k → ∞ We are looking for the large-k solutions of the 2D case without matter (λ = 0). We need to find the zeroes of the determinant of the matrix in equation (4.13). We first look at the 3 × 3 block and calculate it according to the explicit integrals that we have found. Next we substitute the variables as ω = 1, α = 1/2, Ω = −k + x, and where a is a coefficient to be determined and overall the substitution for µ is an educated guess. Except for the choice α = 1/2, everything is still completely general. The unknown frequency to be found is x. Its imaginary part is the growth rate which we are looking for. The parameter m is an effective interaction strength because it gives us µ in this parameterised form. Interaction Strength μ Figure 15. Footprint of the 1D instabilities in the µ-λ plane for k = 0 (homogeneous mode) and α = 1/2 as explained in the text. The colored regions derive from a numerical solution, where the blue footprints correspond to the 2×2 block in equation (3.22), the red solutions to the 1×1 block. The grey regions show the asymptotic solutions in the large-λ limit derived in this appendix.
Next we expand the determinant as a power series for large k and find to lowest nontrivial order det(3×3 block) = 2880 − 480 a − 424 a 2 − 11 a 3 2880 For term proportional to 1/ √ k to dominate we demand the first term to vanish, giving us three possible values for a from the requirement 2880 − 480 a − 424 a 2 − 11 a 3 = 0. The In other words, we have three asymptotic solutions, where one is for positive µ and two for negative µ as expected.
If we now imagine that a is one of these solutions, the first term in the determinant vanishes and in the second term we can substitute a 3 = (2880 − 480 a − 424 a 2 )/11 to remove the a 3 term. In anticipation of the result we further introduce the quantity (D.11) The solution has an imaginary part for 0 < m < m max . Therefore, the large-k footprint of the three instabilities is limited by the lines µ = a i k and µ = a i k + m max,i √ k . (D.12) For µ-values between these lines, the system is unstable. Finally we turn to the 1×1 block in equation (4.13). We proceed with the same substitutions except for µ = a (k + b) , (D. 13) where for the moment we leave open what b is supposed to mean. Expanding the 1×1 block determinant in powers of large k, we here find det(1×1 block) = 6 + a 6 + a 6 (−9 + b + 3x) 1 k + i 2 √ 2 a 3 2(x − 1) 3/2 − (x + 1) 3/2 1 k 3/2 + O(1/k 2 ) . (D.14) Again we can get rid of the first term, this time by setting a = −6, i.e., the footprint of this instability is for negative µ. The leading term does not provide an imaginary solution. In other words, for very large k we do not have an instability. If we keep both the leading and next to leading term, we finally need to solve the equation 16) Solving this equation actually leads to an asymptotic solution where the growth rate exists for a range of b-values. However, the maximum growth rate decreases with 1/ √ k. Therefore, we have overall four instabilities, but for k → ∞ the one from the single block disappears. E Asymptotic solutions for 2D with k = 0 andλ → ∞ We can derive asymptotic solutions for the 2D case with matter, i.e., the large-λ solutions of the eigenvalue equation (4.15), corresponding to the two equations (4.18). We begin with the 2×2 block andλ → +∞. As in the 1D case, we assume that Ω remains of order ω 0 , an assumption which is confirmed by the results. We express the interaction strength in terms of a dimensionless parameterμ in the form µ =μ 1 − αλ . (E.1) To lowest order inλ −1 the eigenvalue equation is, using w = Ω/ω 0 , This result is identical with equation (C.7), but with a different expression for a. To draw the asymptotic footprints we simply need to solve forμ using the limiting a-values given in equation (C.13). The result is shown in figure 16 as grey shaded regions in the upper panels, to be compared with the blue regions which derive from a numerical solution of the full eigenvalue equations. Next we turn to the 1×1 block for the limitλ → +∞ and express the interaction strength in the form µ = −λ +μω 0 log(λ/2ω 0 ) 1 − α .
(E.3) Withμ = 0 the eigenvalue equation is identically fulfilled to lowest order inλ −1 , i.e., to lowest order unstable solutions require µ = −λ/(1 − α). This simple behavior indeed corresponds to the very "thin" footprint shown in red in the upper left panel of figure 16. Includinĝ µ = 0 leads to an approximate eigenvalue equation which is not very simple and does not lead to simple asymptotic solutions. Expressing µ in terms ofμ as in equation (E.3) we can numerially find the growth rate κ as a function ofμ as shown in figure 17. It is clear that the instability footprint in the logarithmic figure 16 will be very narrow. We also notice that the maximum growth rate decreases with increasingλ. (For all of the other instabilities and for α = 1/2, the maximum growth rate κ max = 2 ω 0 .) For the next cases we turn to the limitλ → −∞. In this limit, we write Ω =λ/2 + w ω 0 in analogy to the 1D case. In theλ → −∞ limit, the eigenvalue is characterized by w values of order unity. We also write the interaction strength again in the form of equation (E. Interaction Strength μ Figure 16. Footprint of the 2D instabilities in the µ-λ plane for k = 0 (homogeneous mode) and α = 1/2 as explained in the text. The colored regions derive from a numerical solution, where the blue footprints correspond to the 2×2 block in equation (4.18), the red solutions to the 1×1 block. The grey regions show the asymptotic solutions in the large-λ limit derived in this appendix. where the first expression applies to the 2×2 block, the second to the 1×1 block of the eigenvalue matrix. As before, to draw the asymptotic footprints we solve forμ using the limiting a-values given in equation (C.13). The result is shown in figure 16 as grey shaded regions in the lower panels, to be compared with the blue and red regions which derive from a numerical solution of the full equations.