Nonlinear flavor development of a two-dimensional neutrino gas

We present a numerical survey of the nonlinear flavor development of dense neutrino gases. This study is based on the stationary, two-dimensional ($x$ and $z$), two-beam, monochromatic neutrino line model with a periodic boundary condition along the $x$ direction. Similar to a previous work, we find that small-scale flavor structures can develop in a neutrino gas even if the physical conditions are nearly homogeneous along the $x$ axis initially. The power diffusion from the large-scale to small-scale structures increases with the neutrino density and helps to establish a semi-exponential dependence of the magnitudes of the Fourier moments on the corresponding wave numbers. The overall flavor conversion probabilities in the neutrino gases with small initial sinusoidal perturbations reach certain equilibrium values at large distances which are mainly determined by the neutrino-antineutrino asymmetry. Similar phenomena also exist in a neutrino gas with a localized initial perturbation, albeit only inside an expanding flavor conversion region. Our work suggests that a statistical treatment may be possible for the collective flavor oscillations of a dense neutrino gas in a multi-dimensional environment.

We present a numerical survey of the nonlinear flavor development of dense neutrino gases. This study is based on the stationary, two-dimensional (x and z), two-beam, monochromatic neutrino line model with a periodic boundary condition along the x direction. Similar to a previous work, we find that small-scale flavor structures can develop in a neutrino gas even if the physical conditions are nearly homogeneous along the x axis initially. The power diffusion from the large-scale to small-scale structures increases with the neutrino density and helps to establish a semi-exponential dependence of the magnitudes of the Fourier moments on the corresponding wave numbers. The overall flavor conversion probabilities in the neutrino gases with small initial sinusoidal perturbations reach certain equilibrium values at large distances which are mainly determined by the neutrinoantineutrino asymmetry. Similar phenomena also exist in a neutrino gas with a localized initial perturbation, albeit only inside an expanding flavor conversion region. Our work suggests that a statistical treatment may be possible for the collective flavor oscillations of a dense neutrino gas in a multi-dimensional environment.

I. INTRODUCTION
Due to the mismatch between its weak-interaction and vacuum (mass) states, a neutrino can experience flavor transformation or oscillations even in vacuum (see, e.g., Ref. [1] for a review). In the limit of coherent forward scattering (i.e., with no momentum exchange), a neutrino can experience a flavor dependent refraction in matter through the scattering by the charged leptons and nucleons in ordinary matter [2,3]. Similarly, the neutrinos in a dense neutrino medium can also experience flavor dependent refraction through the neutrino-neutrino scattering [4][5][6]. Unlike in ordinary matter, however, the so-called neutrino self-coupling potential can result in a feature-rich, collective flavor transformation which is coherent among different momentum modes of the neutrino medium (see, e.g., Ref. [7] for a review).
It is beyond the current means to solve the dynamic flavor evolution in an dense, anistropic, inhomogeneous neutrino medium. Instead, simplified models were studied with certain spatial and directional symmetries imposed. For example, the isotropic and homogeneous neutrino gas model has been employed to study the neutrino oscillations in the early universe (e.g., [8,9]), and the stationary, spherical neutrino bulb model was used to study the neutrino flavor transformation in supernovae (e.g., [10,11]). The studies of these restricted models have brought many valuable insights into this intriguing quantum phenomenon including the flavor pendulum model [12] and the analytical understanding of the spectral swaps/splits [13].
The real physical environments such as that inside a core-collapse supernova, however, can be highly anistropic and inhomogeneous. Even if certain directional and spatial symmetries are approximately held in * josh86@unm.edu a neutrino medium, they can be broken spontaneously during the collective flavor transformation [14][15][16][17][18][19][20][21][22] (see also Ref. [23] for a review). Consider, for example, the two-dimensional (2D) neutrino line model in which the neutrinos are emitted from the "neutrino line" along the x axis and propagate in the x-z plane. After linearizing the equation of motion one can show [18] that the amplitudes of the inhomogeneous oscillation modes (whose flavor compositions vary along the x axis) can grow exponentially as functions of z in the linear regime. Unlike the restricted models, it is numerically challenging to solve the neutrino oscillations without the spatial symmetries. As a result, there have been only a few published works devoted to this subject with very limited results [17,20,21].
In this work we presented a numerical survey on the inhomogeneous, nonlinear flavor transformation of dense neutrino gases. This study is based on the simplest 2D neutrino line model with two beams of monochromatic neutrinos and antineutrinos emitted from each point along the neutrino line. Through this work we wish to stimulate the interest of developing an analytic understanding of the nonlinear flavor transformation in the simplest 2D model which may bring further insights into the more realistic models. To keep the problem simple, we will ignore the temporal instabilities [24,25], fast flavor conversions [26,27], the collisional effects [28,29], and many other recent developments (see Ref. [30] for a review).
The rest of the paper is organized as follows. In Sec. II we will describe the two-beam, neutrino line model in detail and review the current analytic understanding of this model. In Sec. III we will briefly describe the numerical method we used to solve the line model and present a survey of the nonlinear flavor transformation of the neutrino gas in this model with various physical parameter choices. In Sec. IV we give our conclusions.

II. THE NEUTRINO LINE MODEL
In this section we first establish the formalism by explaining the neutrino line model proposed in Ref. [18]. We then briefly review the flavor pendulum model [12] which describes the flavor evolution in the neutrino line model when both the spatial and directional symmetries are preserved. We will also briefly review the linearized flavor stability analysis performed in Ref. [18] which describes the flavor evolution in the linear regime.

A. Equations of motion
We study the mixing between two (effective) neutrino flavors ν e and ν τ , with ν τ being a suitable linear combination of the physical ν µ and ν τ . We use the (flavor) polarization vector P p (t, r) in flavor space to represent the flavor quantum state of a neutrino of momentum p at time t and position r. We define the unit (flavor) basis vectors e i (i = 1, 2, 3) such that P 3 = P · e 3 gives the probability P νe of the neutrino being in the electron flavor through the relation and P 1 and P 2 describe the coherence between the two flavors. The polarization vectorP for the antineutrino is defined in a similar way. In absence of neutrino emission, absorption, and collision, the polarization vectors obey the following equations of motion [31] In the above equation,v = p/| p| and ω = δm 2 /2| p| are the velocity and vacuum oscillation frequency of the neutrino, respectively, B = sin(2θ v )e 1 − cos(2θ v )e 3 is a unit vector that denotes the neutrino mixing in vacuum, λ = √ 2G F n e is the matter potential with G F being the Fermi coupling constant and n e the net electron number density, and is the neutrino self-coupling potential with f ν ( p) (ν = ν e , ν τ ,ν e andν τ ) being the initial occupation numbers of the neutrinos of momentum p in the corresponding flavor states.
We adopt the convention that the neutrino masssquared difference δm 2 is always positive and that the normal and inverted neutrino mass hierarchies (NH and IH) are represented by the vacuum mixing angles θ v within the ranges of (0, π/4) and (π/4, π/2), respectively.
We restrict our study to the two-beam, monchromatic neutrino line model in which mono-energetic ν e 's and ν e 's are constantly emitted in two directions,v ± = (±u, 0, v z ), from each point on the x axis, where u = 1 − v 2 z . In this model, the convection operator becomes Because the presence of a large matter density does not suppress collective oscillations in this model other than reducing the effective neutrino mixing angle [12,32], we will work in the appropriate frame rotating about e 3 in flavor space in which where η = sgn(cos 2θ v ) is the signature of the neutrino mass hierarchy, and is the effective oscillation frequency of the neutrino. In the rest of the paper, we will measure the energies in terms of ω eff by setting Equation (2) now reduces to where the subscripts "±" denote the right-and left-going neutrino beams, respectively. In Eq. (8), the parameters and α = nν e n νe (10) measures the strength of the neutrino potential and the neutrino-antineutrino asymmetry, respectively, where n νe and nν e are the number densities of the corresponding neutrino species. We assume that both µ and α are constant in the whole space. From Eq. (8) it is easy to show that the average electron lepton number is constant along z.

B. The flavor pendulum
If the neutrino emission is homogeneous along the x axis and symmetric between the two directions, then In this limit, the two-beam neutrino line model reduces to the bipolar model [32,33], and its flavor evolution is equivalent to the motion of a pendulum in flavor space with a total angular momentum D = P−αP [12]. Defining Q = P + αP + ηe 3 /µ and using Eq. (8) one can show thatḊ whereḞ (z) = v z dF/dz for an arbitrary function F (z), and q = Q/|Q| and (η|Q|)e 3 are the "position" and "weight" of the pendulum bob, respectively. The spin σ = D · q of the bob is conserved, and is the orbital angular momentum of the bob. The pendulum model has been used to obtain many useful insights of the flavor evolution of the bipolar model [12,34]. For example, one can easily see that no (significant) flavor transformation can occur at µ 1 if the neutrino mass hierarchy is normal (η = +1) because the flavor pendulum is near its stable configuration initially. If η = −1 (IH), flavor oscillation will occur only if and P 3 has a minimum value No significant flavor oscillation occurs if µ > 2/(1 − √ α) 2 where the flavor pendulum behaves like a sleeping top [34].

C. Spontaneous symmetry breaking
It turns out that an initially homogeneous neutrino gas can lose its homogeneity during collective oscillations [18,20]. To see this, we assume the following periodic boundary conditions along the x axis: as in Refs. [18,20], where L is the size of the periodic box. Performing the Fourier transformation along the x axis one obtains and a similar expression forP (19b) Because of the coupling among the different Fourier moments, the "power" in the low Fourier moments with small wave numbers k m can dissipate into high moments and result in fine structures and the loss of large-scale, coherent flavor evolution. Before the translation symmetry is badly broken, however, the Fourier moments of different wave numbers evolve independently. This can be seen through the flavor stability analysis performed in Ref. [18] which is recapped below.
In the linear regime where the neutrino coherence fields are small, they obey the linearized equations of motion: The solution to the above equation is a linear superposition of the normal modes: where S In Fig. 1, we demonstrate the largest exponential growth rate κ(m) = max(Im[K (m) a ]) of the flavor coherence amplitude as a function of both the strength of the neutrino potential µ and the Fourier index m for the twobeam neutrino line with α = 0.6 and v z = 1/ √ 2. We also marked the boundaries of the regime where a homogeneous neutrino gas is unstable against flavor conversion [Eq. (15)].

III. NUMERICAL RESULTS
Following the pilot study in Ref. [20] which features a single parameter set (η = +1, α ≈ 0.77, and µ = 13), we have conducted a numerical survey of the two-beam neutrino line model for a wide range of the parameter space.We present a representative collection of our numerical results in this section. All the calculations assume the initial conditions of the following form where | ± (x)| 1. We choose the emission angles of both the left-and right-going beams to be π/4. A change to the neutrino emission angle is equivalent to a change of the neutrino density and a rescaling of the x and z axes simultaneously. We use a relatively large box with L = 30 compared to L ≈ 6.2 in Ref. [20].

A. Numerical method and validation
Unlike Ref. [20] which solved P (m) (19), we solve P ± (x, z) andP ± (x, z) in the direct space. We discretize the x axis into N equal intervals and solve the equations of motion by a finite difference method derived from the Lax-Wendroff algorithm [35]. Our code accurately recovers the pendulum-like flavor oscillation when + = − = const. In the linear regime where the translation symmetry along the x direction is slightly broken, our code produces the correct exponential growth of the coherence amplitudes as predicted by the linearized flavor stability analysis. We varied the number of discrete bins in the x direction and the error tolerance in the numerical integration along the z direction to test the numerical convergence of some representative calculations. We also varied the size of the periodic box L for the calculations with localized perturbations to make sure that the results are independent of the choice of L. We do not enforce the unitary condition, |P(x, z)| = |P(x, z)| = 1, in our code but we rather use it to check the correctness of the numerical solutions.

B. Sinusoidal initial perturbations
We first consider the neutrino gas with a sinusoidal initial perturbation: where ± m are constants. The results from this idealized initial condition demonstrates how the flavor evolution in the two-beam neutrino line model depends on the neutrino mass hierarchy, the neutrino-antineutrino asymmetry, and the neutrino number density.
In Fig. 2 we show P 3 (x, z) with µ = 25 and three values of α for both η = −1 (IH, top panels) and η = +1 (NH, middle panels), where · · · represents the average over the emission directions of the neutrino fluxes. In the bottom panels of the figure, we show the evolution of P 3 when averaged over both the neutrino emission directions and the x axis.
In the IH cases, we set ± m = 0 for all m except Although the left-and right-going neutrino beams have the same initial conditions at any emission point on the x axis, this left-right symmetry is lost at z > 0 when the x-translation symmetry is spontaneously broken. When averaged over the left-and right-going beams, however,  The initial conditions in the NH cases shown in the middle panels of Fig. 2 are the same as in the IH cases Because the homogeneous mode is perturbed symmetrically (i.e., + 0 = − 0 ) and because the flavor pendulum is stable with η = +1, the neutrino gas does not experience significant flavor conversion until the m = 1 mode becomes significant. Similar to the IH cases, the largescale flavor structures also break down into small-scale ones in the NH cases as z increases. But compared to their IH counterparts, P 3 (x, z) with the NH develop more prominent stream-like structures that are localized in the x direction and extended along the z direction. These stream-like flavor structures result in the "streamlines of the neutrino flux" observed in Ref. [20].
The results shown in the bottom panels of Fig. 2 suggest that, although the initial flavor evolution of the neutrino gas can be sensitive to the neutrino mass hierarchy and the initial condition at z = 0, the overall flavor con- version at large z appears to settle down on an equilibrium value which is almost independent of these parameters. In practice, we define P eq 3 andP eq 3 to be the mean values of P (0) 3 (z) and P (0) 3 (z) in the last distance unit of z after they are saturated.
To demonstrate the dependence of P eq 3 andP eq 3 on the neutrino-antineutrino asymmetry and the neutrino density, we plot P eq 3 andP eq 3 in Fig. 3 in terms of α (left panel) and µ (right panel), respectively, both with η = −1 and the initial condition defined in Eq. (28). From Fig. 3 one can see that the antineutrinos in the two-beam line model become almost depolarized in flavor (i.e.,P eq 3 ≈ 0) at large z. The depolarization is generally more complete with a larger value of µ or neutrino density. Because of the conserved electron lepton number L [Eq. (11)], the final mean flavor conversion of the neutrinos is related to that of the antineutrinos through As shown in Fig. 3, the values of P eq 3 in our calculations indeed demonstrate a nearly linear dependence on α but little dependence on µ.
Although the overall flavor conversion of the neutrino gas in the two-beam line model is largely independent of the neutrino self-coupling strength µ, a larger number of x bins is needed to achieve the same accuracy as µ increases (see the changing symbols in the right panel of Fig. 3). This is because more significant small-scale flavor structures are generated at smaller z as µ increases. To demonstrate this interesting feature, we show the flavor evolution of the neutrino gas with the same neutrinoantineutrino asymmetry (α = 0.6) but three different values of µ in Fig. 4. All three calculations assume the inverted neutrino mass hierarchy. The initial flavor perturbations of the neutrino gas in these calculations are the same as what is defined in Eq. (28) except with ± 0 = 10 −3 , ± 1 = 5 × 10 −5 and ± 200 = 5 × 10 −10 (28) for the case with µ = 50.
Comparing the first two panels of Fig. 4 in the first row, one can clearly see that the neutrino gas with a larger value of µ or a higher neutrino density develops finer flavor structures than the one with a lower density. This fact becomes obvious when one compares the magnitudes the Fourier moments of the neutrino polarization vector | P (m) | in these two cases, as shown in the middle panels of Fig. 4. It appears that | P (m) | has a semi-exponential dependence on the Fourier index m for a wide range of m and z: where the exponent β(z) decreases with increasing z. The comparison between the first two panels in the second row of Fig. 4 shows that β is smaller at the same z for the neutrino gas with a larger density which indicates more prominent small-scale flavor structures. The semi-exponential form of | P (m) | in the nonlinear regime is intriguing and cannot be explained by the linearized flavor stability analysis. In the middle panels of Fig. 4 we also plot the largest exponential growth rate κ(m) (dashed curves) predicted by the linearized flavor stability analysis. In the first two panels of this row, one observes that the high moments rise faster than the low moments in order to maintain the semi-exponential power spectrum, while κ(m) decreases with increasing m. In fact, the magnitudes of the Fourier moments with m 200 grow with z for µ = 10 even though they are predicted to be stable in the linear regime (left panel). These results indicate that the excitation of the high mo- ments in these calculations is not due to their own instabilities but because of the "power diffusion" from the low moments. This diffusion phenomenon is reminiscent of the development of the turbulence in an initially laminarly flowing fluid [20] and the kinematic decoherence in an initially isotropic neutrino gas [36].
In the bottom panels of Fig. 4 we also show the evolution of the strengths of a few Fourier moments over z. These panels show that different Fourier moments do not reach equilibrium at the same time. The low moments are saturated first, and the high moments later. This implies that small-scale flavor structures may still develop even after the overall flavor conversion of the neutrino gas has reached its equilibrium value.
The calculations shown in the last column of Fig. 4 is special because the flavor pendulum with α = 0.6 and µ = 50 is in the sleeping-top regime, and the Fourier modes with m 50 are stable (see Fig. 1). Therefore, we intentionally perturbed the Fourier modes with m = ±200 which exhibit the exponential growth in the linear regime (bottom panel). After these moments are excited, they couple the low moments to the high moments and, as a result, | P (m) | has semi-periodic dependence on m initially (middle panel). At large z, however, | P (m) | in this case has a much weaker dependence on m than those shown in the left two panels.

C. Localized perturbations
The sinusoidal initial perturbations discussed in the previous subsection are large-scale perturbations. Next we consider localized initial perturbations of the form where ± g , x 0 , and σ are constants. In Fig. 5 we show the flavor evolution of the neutrino gases with various values of α and µ for both the inverted (IH, top panels) and normal (NH, middle panels) neutrino mass hierarchies. For the initial perturbations, we take ± g = 2 × 10 −3 , x 0 = L/2, for the IH cases, and the same initial perturbations except with for the NH cases. The upper two rows of Fig. 5 show that, for an initial perturbation localized around x 0 , the flavor conversion first appears near x 0 and spreads out to the left and right. At large z the edges of the "envelope" of the flavor conversion region are parallel to the propagation directions of the two velocity modes of the neutrinos. The flavor conversion is coherent in space near these edges, but small-scale flavor structures develop deep inside the envelope. The flavor development inside the envelope shown in Fig. 5 is qualitatively the same as that shown in Figs. 2 and 4. A larger neutrino self-coupling strength or density induces finer flavor structures, and a smaller neutrino-antineutrino asymmetry produces a larger overall flavor conversion.
We intentionally limit our calculations to z < L/2 to avoid the unphysical consequences because of the artificial, periodic boundary condition. As a result, the flavor evolution of the neutrino gases with localized perturbations do not reach an equilibrium in the whole periodic box because unconverted neutrinos continue to stream into the envelope. Instead of averaging over the whole box, we average P 3 over the central region of the envelope and define with ∆L = L/3. In the bottom panels of Fig. 5 we compare the values of P 3 (z) C in both the IH and NH cases and P eq e obtained from the IH calculations with sinusoidal perturbations. Our calculations suggest that the overall flavor conversions in the central region of the envelope are approximately the same as those in the neutrino gases with the sinusoidal initial perturbations. However, we also notice the NH cases tend to develop flavor structures of the shapes of "streams" and "domains" as observed in the previous study [20]. Such regions can have interesting physical consequences if developed in a real astrophysical scenario.

IV. CONCLUSIONS
We have developed a numerical code to solve the 2D (xz), two-beam neutrino line model. We studied the flavor development of the neutrino gases with small sinusoidal and localized perturbations initially (at z = 0).
The neutrino gases with small sinusoidal perturbations behave like a flavor pendulum initially which is coherent in the x direction (along which the gas is almost homogeneous initially) when the neutrino density is below a critical value. This coherent flavor evolution breaks down into small flavor structures as the neutrinos propagate in the z direction. The Fourier analysis shows that the magnitudes of the flavor structures have a semi-exponential dependence on the wave number (in the x direction) and change with z. The overall flavor conversions in the neutrino gases (averaged over the x axis) eventually achieve constant equilibrium values at large z. In our calculations where there are fewer antineutrinos than neutrinos, the antineutrinos almost reach equipartition among different flavors, and the overall flavor conversions of the neutrinos are approximately determined by the neutrinoantineutrino asymmetries according the conservation law of the electron lepton number. This general behavior is largely independent of the neutrino densities or the neutrino mass hierarchy. However, the neutrino gas tend to develop prominent stream-and domain-like flavor structures when the neutrino mass hierarchy is normal. Also the increase of the neutrino density causes the development of more prominent fine flavor structures at smaller z which make the problem more difficult to solve.
The flavor evolution in a neutrino gas with an initially localized perturbation starts from the x coordinate where the perturbation is located and expands afterwards. The flavor development inside the region where the oscillations occur is very similar to that of a neutrino gas with an initial sinusoidal perturbation.
The qualitative results obtained in our numerical survey, such as the semi-exponential power spectrum of the flavor conversion and the final equilibrium conversion probabilities, are very intriguing. However, it remains to be seen whether these results will survive in the more sophisticated models (e.g., the line model with multiple neutrino beams [20] or the ring model with position dependent neutrino densities [21]). If they do, then these results suggest that some simple analytic understanding and statistical treatment may be possible for the flavor oscillations in a multi-dimensional neutrino gas. Such a treatment would be extremely useful in computing the neutrino signals and nucleosynthesis in astrophysical scenarios such as core-collapse supernovae and neutron star mergers.