Multipole expansion method for supernova neutrino oscillations

We demonstrate a multipole expansion method to calculate collective neutrino oscillations in supernovae using the neutrino bulb model. We show that it is much more efficient to solve multi-angle neutrino oscillations in multipole basis than in angle basis. The multipole expansion method also provides interesting insights into multi-angle calculations that were accomplished previously in angle basis.

Unlike the MSW effect, collective neutrino flavor transformation caused by neutrino selfinteraction involves the quantum states of neutrino themselves and, therefore, is nonlinear. Except for a few very simplistic models, the flavor evolution of a dense neutrino gas can be solved only through numerical methods. This non-linearity together with the inhomogeneous and anisotropic physical environment makes it a daunting task to solve collective neutrino oscillations in a (core-collapse) supernova even numerically. The most sophisticated calculations of this kind so far adopted the "neutrino bulb model" [7] which assumes the spherically symmetry for the supernova environment. The current approach of solving neutrino oscillations in the neutrino bulb model is to discretize the neutrino emission (zenith) angle as well as the neutrino energy and solve millions to tens of millions of coupled differential nonlinear equations simultaneously. Although this "angle-bin method" is straightforward to implement, a severe drawback of this approach is that a large number of angle bins ( 1000) are required to achieve numerical convergence even in the regime where no significant neutrino oscillation is observed [7,22].
Because of its spherical symmetry, the neutrino bulb model does not completely capture the complexity of the supernova environment revealed in recent multi-dimensional supernova simulations (e.g. [23][24][25]). Even if the supernova is approximately spherically symmetric, the spherical symmetry in neutrino oscillations can be broken spontaneously because of the vectorvector coupling nature of the neutrino self-interaction [19,26,27]. In the cases where the spherical symmetry is indeed maintained, the current technique is still insufficient for solving neutrino oscillations in the presence of the neutrino halo [17] and/or transition magnetic moments of Majorana neutrinos [20,21]. Although a new computing model for supernova neutrino oscillations is yet to be found, it seems clear that the angle-bin method needs to be replaced by a more efficient approach.
A multipole expansion approach using Legendre polynomials was employed in an earlier study with the focus on the "kinematical decoherence" of a homogeneous neutrino gas [28]. A different moment expansion method was later proposed for the neutrino bulb model [29]. However, the latter approach seems unnecessarily complicated and has not been widely accepted.
In this paper we propose a multipole expansion method similar to that in [28] but for the neutrino bulb model. The rest of the paper is organized as follows. In Section 2 we explain the formalism of the new multipole expansion method. In Section 3 we validate this method and show that it can be much more efficient than the traditional angle-bin method. We also discuss some physical insights behind the multipole expansion method and how it may be improved. In Section 4 we give our conclusions. Some of the approximations used in our approach are listed in Appendix A.

Formalism
Here we consider two-flavor neutrino oscillations (i.e. ν e ↔ ν x andν e ↔ν x ) in the neutrino bulb model. Generalization to three neutrino flavors is straightforward. The neutrino bulb model assumes the spherical symmetry about the center point of the supernova and the azimuthal symmetry about any radial direction from this center point. In this model the flavor density matrix ρ E,u (r) of the neutrino at radius r depends only on the energy E and trajectory u = cos(2ϑ 0 ) of the neutrino, where the (zenith) emission angle ϑ 0 is defined with respect to the normal of the neutrino sphere at radius R. We adopt the convention of the neutrino flavor isospin [6] so that the diagonal elements of ρ E,u (in flavor basis) with E < 0 are proportional to the negative number densities ofν e andν x with energy |E|. The flavor density matrices of neutrinos and anti-neutrinos are normalized in the following manner: where Φ ν and Φν are the total number fluxes of the neutrino and anti-neutrino, respectively. The equation of motion for density matrix ρ E,u is where ∂ r is differentiation with respect to radius r, is the radial component of the neutrino velocity, and is the Hamiltonian. In the above equation M 2 is the neutrino mass-squared matrix in flavor basis, G F is the Fermi coupling constant, and L is the matrix of net charged-lepton number densities. Here, for simplicity, we have assumed that neutrinos are emitted isotropically from the neutrino sphere. Solving eq. (2.2) directly usually involves numerical integration of millions to tens of millions of coupled, nonlinear differential equations. One of the reasons for the necessity of solving so many equations is that the solution to eq. (2.2) with insufficient number of discrete angle bins, say N a 1000, has spurious flavor instabilities at small radii and, therefore, are qualitatively different from the solution in the continuous limit N a → ∞ [7,22]. As a result, a large number of angle bins is necessary to ensure that no spurious oscillation occurs in the region where nothing really happens. This nuisance can be avoided by solving the equation of motion in multipole basis.
We define the n'th multipole/moment of the neutrino flavor matrix to be where P n (u) are the standard Legendre polynomials. The reason that we use u = cos(2θ 0 ) instead of v u (as in [28,29]) in the definition of the multipoles is that u has a fixed range which changes with r. Using the normalization condition Instead of transforming eq. (2.2) to multipole basis directly, we note that (multi-angle) neutrino oscillations in the neutrino bulb model usually occur at Therefore, we first expand eq. (2.2) in terms of z as in [12,30]. Using the appropriate lowest order approximation for each term in eq. (2.4) we obtain (see Appendix A) where with n e being the net electron number density, σ 3 is the third Pauli matrix, and In eq. (2.10), ω ≡ ∆m 2 /2E is the vacuum oscillation frequency of the neutrino with ∆m 2 being the mass squared difference of neutrino, and θ eff 1 is the effective mixing angle of neutrino in matter [6,8].
Using the above approximation and also identity we obtain where a n = n + 1 Solving eq. (2.15) can be much more efficient than solving eq. (2.2). In the region where no significant neutrino oscillation occurs, all multipoles with n > 0 are very small. Therefore, one can set all but the first two multipoles to 0 initially and solve eq. (2.15) for ρ E,0 and ρ E,1 only until ρ E,1 become significant. 1 (Solving eq. (2.15) for ρ E,0 only is equivalent a single angle approximation, which may lead to wrong results [16].) More moments can be added adaptively from this point on as necessary.

Validation and discussion
To validate and show the usefulness of the multipole expansion method we solve eq. (2.15) numerically with two sets of initial neutrino energy spectra. We assume that the matter density is not large enough to suppress collective neutrino oscillations, and we take λ = 0. As comparison we also solve eq. (2.2) numerically with replacement In both kinds of calculations we use mass-squared difference ∆m 2 = −3 × 10 −3 eV 2 (i.e. with the inverted neutrino mass hierarchy), effective mixing angle θ eff = 0.01, 100 energy bins per neutrino flavor, and the radius of the neutrino sphere R = 11 km.
The two sets of neutrino spectra used in our calculations are adapted from [7] and [13], respectively, which are known to produce a single swap/split (SS) and multiple swaps/splits (MS) in final neutrino fluxes. In both sets the neutrino spectra are described by the Fermi-Dirac distribution with the parameters listed in table 1.
In figure 1 we show the neutrino fluxes at 200 km computed in multipole basis with 25 mutipoles (i.e. with ρ E,n>25 = 0) and with the SS spectra. In the same figure we also plot the differences between these results and those computed in angle basis with 1200 angle bins. One can see that these two calculations agree with each other very well. We note that the SS spectra MS spectra L νe (10 51 ergs/sec)  small differences between the two results are partly due to the approximations we made in eq. (2.9) which is employed in the multipole expansion method. To understand why so few multipoles are needed, we consider the strengths of the multipoles defined as follows, outside the neutrino sphere, Tr(ρ E,u ) and, therefore, Tr(ρ E,n ), do not change with radius. Furthermore, because we do not consider the quantum decoherence of the neutrino states, Tr(ρ 2 E,u ) are also constant. Using eqs. (2.6) and (2.7) it is straightforward to show that d dr n S 2 E,n = 0.
In figure 2 we show how evolve as functions of radius. One can see that, right after collective neutrino oscillations begin, both S n andS n grow exponentially with radius. This result suggests that all multipoles become unstable simultaneously which one may expect from the ansatz of the linear stability analysis in angle basis [30]. However, this result seems to be contrary to a previous study of the homogeneous gas of neutrinos where higher multipoles are populated successively by diffusion from lower multipoles [28]. From figure 2 one can also see that S n andS n are never large for n > 0. This property together with eq. flavor depolarization in this particular case. This result is, of course, already known from figure 3 in [32], which motivated us to solve neutrino oscillations in multipole basis in the first place. From figure 2 one may conclude that only the first very few multipoles, far fewer than 25, are needed for this calculation. This is indeed true but only up to a certain radius after which spurious oscillations of big amplitudes would occur. This phenomenon implies that some kind of diffusion among multipoles (as suggested in [28]) may indeed exist while collective oscillations are active. Such spurious oscillations may be suppressed by implementing an appropriate closure scheme (as e.g. in [33]). We note that any closure that insures the accuracy of the first multipoles would be sufficient for most physical applications at r R. We also performed the calculations for the MS spectrum case. It was reported in [13] that one needed 15000 angle bins to achieve "apparent convergence" for this case. We found that true convergence were achieved with 50000 angle bins in our calculation. In comparison only 300 multipoles are required to achieve numerical convergence in multipole basis. In figure 3 we show the neutrino fluxes at 400 km in these calculations.
In figure 4 we show the multipole strengths S n andS n in the MS spectrum case as functions of radius. Although the multipoles in this calculation also show an overall behavior of exponential growth right after collective oscillations begin, the first several multipoles oscillate with radius, and their strengths are of similar magnitude. For this reason we plot S n andS n with n = 1, 5, 10 and 15 instead. Because the strengths of the multipoles do not decrease rapidly with increasing n, significant more multipoles are required to achieve numerical convergence in the MS spectrum case than in the SS spectrum case.

Conclusion
We have developed and validated a multipole expansion method for the multi-angle calculations of collective neutrino oscillations in supernovae. For the neutrino spectra that we tested this method requires solving far fewer number of equations (by approximately two orders of magnitude) than the traditional angle bin method does. In particular, one needs to solve only the equations of motion for the first two multipoles before collective neutrino oscillations begin. The multipole expansion method can help us understand the emergence of the angular structure in collective neutrino oscillations. We found that the strengths of all but the first multipoles of the neutrino fluxes grow approximately exponentially right after collective neutrino oscillations begin. We also demonstrated that the strengths of the multipoles decrease relatively slowly with increasing multipole index n when multiple swaps/splits appear in the final neutrino energy spectra. Such cases can require significant more multipoles to be solved to achieve numerical convergence than in the case where a single swap/split exits.
Although our calculations were performed for two flavor neutrino oscillations with the neutrino bulb model, this method can be easily generalized to three flavor calculations and can potentially be extended to more realistic supernova models. The efficiency of the multipole expansion method may be further improved if an appropriate closure scheme is employed.

A Approximations
Like in previous works (e.g. [12,30]) we expand the Hamiltonian in eq. (2.4) in terms of z = R 2 /4r 2 and keep the leading order terms for each part of the Hamiltonian. Higher order terms can be included for more accuracy. For vacuum Hamiltonian we have where ω = ∆m 2 /2E is the vacuum oscillation frequency with ∆m 2 being the neutrino masssquared difference, and θ v is the neutrino vacuum mixing angle. We have dropped the trace term which has no impact on neutrino oscillations. We note that the dispersion of H vac /v u in neutrino energy E, which is much larger than its dispersion in angle parameter u, plays an essential role in collective neutrino oscillations [6]. We also note the terms of order z or higher are unlikely to be important because H vac is smaller than H νν where collective oscillations are active.
For matter Hamiltonian we have where λ(r) = G F n e (r)z(r) and where we have again dropped the trace term. The terms on the right hand side of eq. (A.2) that do not depend on u have no effect on collective neutrino oscillations other than resetting the effective mixing angle θ v to a smaller value θ eff [6,8]. The angle dependent term in eq. (A.2) can lead to suppression of collective oscillations in the case of very high matter density, which can be important near the neutrino sphere and/or during the accretion phase of the explosion [12]. The terms of order z 2 or higher in eq. (A.2) are unlikely to be important because the matter density itself decreases rapidly with r. In the regime where the matter density is much larger than G −1 F |ω| for typical neutrino energies, we obtain The angle dependence of H νν can lead to the multi-angle suppression of collective oscillations [16,30].