Fast flavor instabilities and the search for neutrino angular crossings

With the recognition that fast flavor instabilities likely affect supernova and neutron-star-merger neutrinos, using simulation data to pin down when and where the instabilities occur has become a high priority. The effort faces an interesting problem. Fast instabilities are related to neutrino angular crossings, but simulations often employ moment methods, sacrificing momentum-space angular resolution in order to allocate resources elsewhere. How can limited angular information be used most productively? The main aims here are to sharpen this question and examine some of the available answers. A recently proposed method of searching for angular crossings is scrutinized, the limitations of moment closures are highlighted, and two ways of reconstructing angular distributions solely from the flux factors (based respectively on maximum-entropy and sharp-decoupling assumptions) are compared. In (semi)transparent regions, the standard closure prescriptions likely miss some crossings that should be there and introduce others that should not.


I. INTRODUCTION
Calculations predict that supernova and neutron-starmerger neutrinos experience fast flavor instabilities at some times and locations, possibly with important consequences for the dynamics, composition, and neutrino signal . Of the physical conditions that affect stability, the electron-lepton-number (ELN) angular distribution G(µ) is recognized as being particularly important: where g νe and gν e are the neutrino and antineutrino distribution functions, E ν is neutrino energy, and µ = cos θ is the angle of propagation away from the radial direction. (Axial symmetry is assumed and differences among the heavy-lepton flavors are ignored.) Angular crossings of G(µ)-points where the function passes through 0are associated with fast instability.
Since the presence of angular crossings is a simple criterion, searching for crossings is a useful alternative to running every single location through linear stability analysis. But while the criterion is straightforward to apply if the neutrino angular distributions are known, they usually are not. Because of the high complexity of supernova and merger physics, sacrifices simply have to be made in some aspects of the problem, even when using the most powerful supercomputers. Computing neutrino transport with high angular resolution in momentum spacethat is, finely resolving thep-dependence of the neutrino distributions-is expensive, and many radiationhydrodynamics simulations opt for evolving only the lowest moments of the angular distributions (the energy densities and fluxes, for example). Thus we face the question * NASA Einstein Fellow (ljohns@berkeley.edu) that occupies us here: How can this limited amount of information best be used in the search for angular crossings?
Recently a few publications have advocated and adopted a search procedure that involves applying polynomial weighting functions to G(µ) [42][43][44]. In Sec. II we examine this approach. We find that it never introduces spurious crossings but sometimes misses physical ones. In our view, more evidence is needed to establish whether this approach has any advantages over searching point-by-point for crossings in G(µ).
In Sec. III, inspired by Ref. [42]'s contrasting of the polynomial method with other criteria found in the literature, we offer our own remarks on the similarities and differences.
Regardless of how the search is conducted, it inevitably runs into problems in (semi)transparent regions if only the first few angular moments are employed. In Sec. IV we demonstrate this point, emphasizing that a closure prescription can be adequate for energy-momentum transport while being inadequate for angular-crossing searches. The paradigmatic example is a highly forwardpeaked distribution. This case is easy to describe from the transport perspective: particles travel like a beam, and pressure is nonzero only along the flux direction. But the same distribution flusters angular moments, as very many of them are required to specify an approximation that is not compromised by unphysical features like negative particle numbers.
In this paper we use completion to mean a prescription for reconstructing angular distributions from the first few moments (or, equivalently, generating all higher moments). We use this term in contrast with closure, by which we mean a prescription for generating a finite number of higher moments. In practice, closures usually only specify moments that are one or two ranks above those that are actually being evolved, because that suffices to close the hierarchy of coupled transport equations. Our goal, however, is not to approximate transport by trun-cating the hierarchy of equations: our goal is to fill in the details of the angular distributions. We would prefer to complete, rather than close, the angular distributions that are output by moment-based simulations.
Here we focus on two completions that depend only on the neutrino flux factor f = F/E, where F is the energy flux and E is the energy density. The maximum-entropy completion fills in the details by assuming that the angular distribution is of an entropy-maximizing form. The same distribution types have been used to motivate closures; here we are simply adopting them for a different (but obviously related) purpose. We compare this first prescription with the sharp-decoupling completion, which maps an angular distribution with flux factor f onto the angular distribution in the supernova bulb model with the same flux factor. These two completions are presented in Sec. V.
We also find in Sec. V that, given the same 0th and 1st moments, closed and completed distributions do not always agree on whether a crossing occurs. The disagreement can go either way, with closed distributions exhibiting a crossing not found with completed distributions or vice versa. In view of this fact, we urge caution when presenting or interpreting stability analyses that apply closures to obtain the angular distributions.
The maximum-entropy and sharp-decoupling completions have their own strengths and weaknesses, and we are not endorsing either one as a final answer to the question of how best to search for angular crossings using moment data. In Sec. VI we indicate the steps we will take elsewhere to provide a more actionable answer to this question.

II. POINTS VS. POLYNOMIALS
In the method proposed in Ref. [42] and utilized in Refs. [43,44], a weighting function is introduced, such that F(µ) > 0 for all µ ∈ [−1, 1]. The coefficients a n are arbitrary aside from having to satisfy the positivity constraint, and N labels the highest angular moment for which simulation data is available. A weighted version of the ELN angular distribution is then formed, and the search for angular crossings is conducted by searching for functions F(µ) that satisfy where The idea is that if there is an angular crossing, a wellchosen polynomial F(µ) will exaggerate the region where the crossing occurs, causing the integrated quantity I F to have a different sign than I 0 . We call this search procedure the polynomial method.
An alternative is to simply scan over G(µ), evaluating the function at different values of µ and checking if it changes sign (the point-by-point method ). Ref. [42] does not explain why (or whether) the polynomial method is superior to the point-by-point method-the latter procedure is not mentioned-but the paper does observe, just as F(µ) is being introduced (footnote 2), that an angularmoment expansion cannot be used to approximate G(µ) because higher moments are not necessarily negligiblethey simply are not provided by moment-based simulations. We agree with this point and develop it further in Sec. IV.
However, the polynomial method does not evade the issue of non-negligible higher moments. The product F(µ)G(µ) that appears in Eq. (3) is a kind of reshaped version of G(µ), and so it might be thought that the method in some way explores a fuller range of ELN angular distributions (because of the flexibility in F) while at the same time being informed by the simulation's output (because the starting point is G). This is not the case, unfortunately. Let G (µ) ≡ F(µ)G(µ) be the reshaped ELN angular distribution. Then and the criterion I F I 0 < 0 becomes I 0 I 0 < 0, which simply checks if G(µ) is capable of being reshaped in such a way that the overall lepton number changes sign. The polynomial method is thus in no way akin to searching for crossings in physically motivated modifications of the moment-reconstructed distribution. In fact the procedure is in quite the opposite spirit: the distribution is reshaped so as to artificially exaggerate crossing regions. The polynomial and point-by-point methods are equally limited by truncation of the angular-moment expansion. Moreover, the polynomial method sometimes misses angular crossings. (Because F(µ) > 0, spurious crossings are never created by the search procedure itself, although they may appear due to neglecting higher moments in the ELN angular distribution.) As an example of missed crossings, let N = 1, define by and consider the case with I 0 > 0. Since a crossing exists in the backward direction for > 0. The criterion I F I 0 < 0 can be rewritten as I F /I 0 < 0, giving [for crossing to be found] .  7)). The angular crossing at µ ∼ = −0.7 is not detected by the polynomial search.
Assume that a 0 > 0. Then Eq. (9) entails This says that if 0 < < 2/3, then F(µ = +1) must drop below 0 in order to satisfy I F I 0 < 0. If instead a 0 < 0, then F(µ = −1) < 0 for all > 0. Therefore many ELN angular distributions with crossings will not be identified as such in the polynomial search. The requisite polynomials are excluded by the positivity condition. The idea being demonstrated here is that because F(µ) is required to be nonnegative, it is limited in how much it can amplify certain regions relative to other ones. The ELN angular distribution itself is not limited in this way, and so small enough crossings can go unnoticed. This intuitive point is illustrated in Fig. 1. F(µ), the purple curve, weights the backward direction as much as possible relative to the forward direction, hoping to amplify the crossing region and make I F go negative, but F(µ)G(µ) can only be distorted so much. The area under the solid green curve remains positive.
Based on the considerations above, we are unable to see any advantages to the polynomial method as compared to scanning point-by-point over G(µ). The latter is in fact equivalent to using delta functions in place of F(µ), and because the philosophy of the polynomial method is that crossing regions should be amplified as much as possible, delta functions would intuitively seem to be the best way to accomplish this. The point-by-point method could also very well be more efficient, since the scan is directly over µ rather than over N th-degree polynomials with coefficients satisfying a constraint that must itself be satisfied at all µ ∈ [−1, 1].
The possibility certainly remains that we have overlooked something in our assessment. Since the polynomial method has now been used in multiple publications on fast instabilities, it would be a valuable contribution to the literature for any advantages we may have overlooked to be spelled out.

III. ZERO MODES, PENDULUMS, & POLYNOMIALS
Seeing as how several other instability criteria can be found in the literature, it may be worthwhile to clarify how they all relate to one another.
The most rigorous test of stability involves linearizing the equations of motion-more particularly, making them linear in the off-diagonal elements of the density matrices-and searching for normal-mode solutions of the form e −iΩt+iKx , where for simplicity we assume a single spatial dimension [47,48]. The result is a dispersion relation that involves integrals such as where w = w(Ω) and z = z(K). Generally the dispersion relation is therefore transcendental and must be solved numerically.
Ref. [10] pointed out that by picking K such that z = 0, the dispersion relation simplifies from a transcendental equation to a quadratic one, allowing one to solve for Ω(K 0 ) analytically. (K 0 , the zero mode, is the value of K for which z = 0.) The criterion Im Ω = 0 can then be written as Ref. [24] derived two different criteria that apply to homogeneous (K = 0) fast instabilities. The first of these comes from the demand that there be a trajectoryṽ ∈ [−1, 1] that satisfies a certain kind of resonance condition in the linear regime. If such a trajectory exists, the system is unstable. The criterion is simply The other criterion comes from the pendulumlike nature of fast flavor conversion and reads In Ref. [24] these criteria served the purpose of demonstrating the validity of the paper's analysis, showing in particular how the analysis allows for criteria to be formulated that improve (in the homogeneous setting considered) upon the angular-crossing criterion. They can be viewed as being complementary to the zero-mode test of Ref. [10], but one should be cautious in employing them on simulation data because, unlike the zero-mode test, they do not derive from exactly solving a particular mode's dispersion relation. One should also be cautious in applying the zero-mode test, of course, since it only assesses the stability of a single mode.
Ref. [42] claims that the polynomial search "offers an infinite number of inequalities similar to" Eq. (12) and suggests that this particular inequality corresponds to F ± (µ) = µ 2 ± 2µ + 1. While it is true that this expression does not generally show up in the search because What can be said in favor of the polynomial search is that if the zero mode is unstable, then it will register a crossing. If I F+ I F− < 0, then one of the two factors must have the opposite sign of I 0 . (Since F + (µ = 1) = 0, a nuance here involves whether the polynomials must be non-negative or strictly positive.) But suppose that the factors I F+ and I F− are both negative. Then, although the zero mode is stable, the search registers a crossing anyway. Testing each factor individually is not equivalent to evaluating their product. Thus the polynomial method neither generalizes nor even reproduces the zero-mode criterion (or for that matter the criteria in Eqs. (13) and (14)). The difference between the two is in fact a crucial distinction. Eq. (12) implies that the discriminant of the quadratic dispersion relation is negative, hence that Im Ω = 0 for a particular mode. The criterion I F I 0 < 0 lacks any such interpretation.
To summarize, the following instability tests are found in the literature: 1. Linear stability analysis for all (or a range of) Ω and/or K.  The presence of an angular crossing has been shown to be necessary and sufficient for the existence of an instability Im Ω = 0 for some wave vector [46]. Besides that equivalence, no two of these tests produce identical results in all situations. The point-by-point and polynomial methods are similar in that they adopt angular crossings as the instability criterion, differing in how they search G(µ) for crossings. But it is a separate question, which we now turn to, how G(µ) is constructed in the first place.

IV. CLOSURES & COMPLETIONS
We are drawing a distinction in this paper between closures (which supply a finite number of angular moments) and completions (which supply an infinite number) of the neutrino distribution functions. The distinction is helpful here because it reflects the differing needs of transport calculations and fast-instability searches.
Much attention has been devoted to closure in the context of radiative transfer (e.g., Refs. [49,50] and the studies they cite). Its importance can be seen immediately by taking angular moments of the Boltzmann equation (here assuming flat spacetime with spherical symmetry and a static background): and so on ad infinitum. (As before, E and F are the energy density and flux density. P is the radiative pressure.) Particle streaming alone couples the moments in an infinite hierarchy of equations. With a savvy implementation of closure, an accurate solution is obtained for the transport of energy and momentum. Many different closure methods have been proposed, implemented, and analyzed. Some of themvariable Eddington tensor methods-solve the Boltzmann equation under certain simplifications to all orders in the moment expansion and could be used to output approximate angular distributions. We are interested, however, in simulations that in practice only output the lowest moments, regardless of how they are obtained. These include simulations that adopt algebraic Eddington tensor methods.
In such cases, accuracy in E and F does not guarantee accuracy in the reconstruction of angular distributions using solely those moments. In other words, while it may not be necessary to track higher moments for the purposes of approximately solving the Boltzmann equation, higher moments must sometimes be supplied if the angular distributions are to be approximated well.
This point is illustrated in Fig. 2, which shows the angular distributions ψ that result from applying two different M1 closure prescriptions. Three different flux factors f are shown for each closure. In one spatial dimension, the radiative pressure is P = pE, where p is the Eddington factor. Shown in the figure are the Levermore closure [51] p = 3 + 4f 2 and the Minerbo closure [52] These prescriptions follow from distinct physical assumptions. In the Levermore case, closure is obtained by asserting that the pressure is isotropic in the frame in which the flux vanishes. In the Minerbo case, the pressure is calculated by assuming that the underlying distribution function is classical and entropy-maximizing. From Fig. 2 we see that the distribution functions do not radically differ for the plotted flux factors. In particular, they both exhibit the same unphysical features at f 0.5: negative particle densities in non-radial directions and positive but highly inflated particle densities in the backward direction.
Accurate modeling of forward-peaked distributions demands that higher moments be specified.

V. MAXIMUM-ENTROPY & SHARP-DECOUPLING DISTRIBUTIONS
To illustrate the latitude in modeling angular distributions-even with the densities and fluxes taken as given-we examine two completions in this section. For simplicity we consider monochromatic neutrinos, so that the angular distributions of number density and energy density have the same shape.
The first is the maximum-entropy distribution ψ (max) , which maximizes the entropy s = ψ log ψ at fixed E and F [52,53]. It has the form where k = 0. The distribution takes on k = +1 (Fermi-Dirac statistics) or k = −1 (Bose-Einstein) if the entropy functional is modified accordingly, but here we focus strictly on the classical limit as a representative case. For a distribution of this form, the flux factor is calculated to be which can be numerically inverted to find a(f ). The other parameter η sets the overall normalization of ψ (max) . Normalizing so that the 0th moment is unity, we have The maximum-entropy assumption is the origin of the Minerbo closure and its quantum-statistical variations . Flux factor f of the bulb distribution as a function of radius r in units of the neutrinosphere radius Rν . Note that f is never below 0.5. Compare to Fig. 6 of Ref. [57], which shows the flux factor vs. radius in simulations with Boltzmann neutrino transport. [54][55][56]. Although ψ (max) models astrophysical neutrino angular distributions imperfectly, it does capture some of the essential features. An alternative way to specify all higher moments is to use a sharp-decoupling completion. A distribution of this type is found in the bulb model of supernova neutrino emission, in which all neutrinos decouple at a single sharp surface (the neutrinosphere) at radius R ν [58]. The distribution is where r is the radial coordinate, θ is the step function, and A is a normalization factor. We choose which ensures as before that the 0th moment is unity. The flux factor in the bulb model is in one-to-one correspondence with the radial distance r: or, inverting, The flux factor f (bulb) is plotted in Fig. 3 as a function of r/R ν . It starts at the neutrinosphere with a value of 0.5 and by twice the emission radius already exceeds 0.9. Because neutrinos are radiated strictly outward, ψ (bulb) is intrinsically forward-peaked. These completions are both suitable for certain applications, but they differ markedly. The left panel of  Another way to inspect the angular information contained in these completions is by observing their angular moments These are shown in the right panel of Fig. 4. In order to write a recursion formula for ψ (max) , we note that a −a dx x n e x = a n e a + (−1) n e −a − n a −a dx x n−1 e x .
This then leads to Note that with the normalization of Eq. (22), ψ (max) n → 1 as a → ∞ (corresponding to f → 1) at finite n. In the other limit, as a → 0, ψ (max) n alternates between 0 (at odd n) and finite but declining values (at even n). This is because all even n have some contribution from the isotropic moment, whereas all odd n vanish by symmetry. If angular distributions are to be expressed using truncation, it is better to truncate an expansion in orthogonal functions (e.g., Legendre polynomials) than an expansion in µ n moments. Closure prescriptions make the distinction moot, but only in part, by ensuring agreement in the isotropic limit.
For the sharp-decoupling completion, using the nor-malization of Eq. (24), we have As the figure makes clear, forward-peaked distributions converge slowly in n and vary in exactly how they do so. In Sec. IV we saw that standard M1 closures become unphysical at flux factor f 0.6. By comparing the Minerbo closure and the max-entropy completion, Fig. 5 shows that the failure of closure to accurately describe forward-peaked angular distributions leads to closures being unreliable when used in the search for angular crossings.
The figure shows potential errors of two kinds. In the left panel, an angular crossing in the completed distributions is missing from the closed ones at µ just above 0.9. The reason for this crossing being absent is, again, that closures struggle with forward-peaking. In order for closure to place a lot of weight at µ ∼ 1, it must impose unrealistic negative fluxes in more transverse directions, as we saw above. The right panel of Fig. 5 shows that this undesired consequence can introduce spurious crossings. Even in cases where closure and completion both have crossings, a single crossing in the completion is typically turned into a double crossing in the closure.
While ψ (max) is not necessarily an optimal approximation, it does capture the main feature-concentrated flux at µ ∼ 1-that underlies the inadequacies of closure. We conclude that an angular-crossing search using closure is susceptible to these errors.

VI. DISCUSSION
Several recent publications have searched for fast instabilities, using angular crossings as a proxy. This manner of searching faces technical challenges associated with the limited angular resolution of simulations. The root of the problem is simple: as angular distributions become more forward-peaked, more angular moments are required for their accurate approximation.
Closure-which, in practice, usually means that the angular-moment expansion is truncated at the pressure or heat-flux tensor-is inadequate for forward-peaked distributions. Completions-where the available moments are matched to a physically motivated angular spectrum-are a superior alternative, but the question arises exactly which physics should be used as the motivation. As their names indicate, the sharp-decoupling and maximum-entropy distributions are based on assumptions regarding decoupling and entropy. They are easy to formulate and work with, but they are not necessarily the best choices.
Elsewhere we will present phenomenological fits to the neutrino angular distributions in a supernova simulation using Boltzmann transport [59]. From these fits, a completion can be constructed that better reflects the actual distributions found in supernovae.