Fast Neutrino Flavor Instability in the Neutron-star Convection Layer of Three-dimensional Supernova Models

Neutrinos from a supernova (SN) might undergo fast flavor conversions near the collapsed stellar core. We perform a detailed study of this intriguing possibility, analyzing time-dependent state-of-the-art 3D SN models of 9 and 20 Msun. Both models were computed with multi-D three-flavor neutrino transport based on a two-moment solver, and both exhibit the presence of the lepton-number emission self-sustained asymmetry (LESA). The transport solution does not provide the angular distributions of the neutrino fluxes, which are crucial to track the fast flavor instability. To overcome this limitation, we use a recently proposed approach based on the angular moments of the energy-integrated electron lepton-number distribution. With this method we find the possibility of fast neutrino flavor instability at radii<~20 km, which is well interior to the neutrinosphere. Our results confirm recent observations in a 2D SN model and in 2D/3D models with fixed matter background, which were computed with Boltzmann neutrino transport. However, the flavor unstable locations are not isolated points as discussed previously, but thin skins surrounding volumes where electron antineutrinos are more abundant than electron neutrinos. These volumes grow with time and appear first in the convective layer of the proto-neutron star (PNS), where a decreasing electron fraction (Ye) and high temperatures favor the occurrence of regions with negative neutrino chemical potential. Since Ye remains higher in the LESA dipole direction, where convective lepton-number transport out from the nonconvective PNS core slows down the deleptonization, flavor unstable conditions become more widespread in the opposite hemisphere. This interesting phenomenon deserves further investigation, since its impact on SN modeling and possible consequences for SN dynamics and neutrino observations are presently unclear. (abridged)

It has been speculated for a long time that neutrinos from a supernova (SN) may undergo fast flavor conversions near the collapsed stellar core. We perform a detailed study of this intriguing possibility, for the first time analyzing two time-dependent state-of-the-art three-dimensional (3D) SN models of 9 M and 20 M from recent papers of Glas et al. Both models were computed with multi-dimensional three-flavor neutrino transport based on a two-moment solver, and both exhibit the presence of the so-called lepton-number emission self-sustained asymmetry (LESA). The transport solution does not provide the angular distributions of the flavor-dependent neutrino fluxes, which are crucial to track the fast flavor instability. To overcome this limitation, we use a recently proposed approach based on the angular moments of the energy-integrated electron lepton-number distribution up to second order, i.e., angle-energy integrals of the difference between νe andνe phasespace distributions multiplied by corresponding powers of the unit vector of the neutrino velocity. With this method we find the possibility of fast neutrino flavor instability at radii smaller than ∼20 km, which is well interior to the neutrinosphere where neutrinos are still in the diffusive and near-equilibrium regime. Our results confirm recent observations in a 2D (axisymmetric) SN model and in 2D and 3D models with fixed matter background, which were computed with Boltzmann neutrino transport. However, the flavor unstable locations are not isolated points as discussed previously, but thin skins surrounding volumes whereνe are more abundant than νe. These volumes grow with time and appear first in the convective layer of the proto-neutron star (PNS), where a decreasing electron fraction and high temperatures favor the occurrence of regions with negative neutrino chemical potential. Since the electron fraction remains higher in the LESA dipole direction, where convective lepton-number transport out from the nonconvective PNS core slows down the deleptonization, flavor unstable conditions become more widespread in the opposite hemisphere. This interesting phenomenon deserves further investigation, since its impact on SN modeling and possible consequences for SN dynamics and neutrino observations are presently unclear.

I. INTRODUCTION
The deepest supernova (SN) regions provide a unique laboratory to probe neutrino flavor conversions in a nonlinear regime, where the neutrino evolution is determined mainly by their mutual interactions. Indeed, at distances r O(10 2 ) km from the centre of the SN, the neutrino density n ν is so high that it dominates the flavor evolution, leading to self-induced neutrino flavor conversions [1][2][3]. These have been a topic of intense investigation for over a decade [4][5][6][7][8][9][10]. See refs. [11][12][13][14] for recent reviews.
This possibility of fast flavor conversions and its potential effects on SN dynamics and nucleosynthesis has stimulated several studies to assess the occurrence of fast instabilities in different SN models. A first study in this direction was performed in [29], where a dedicated analysis of the angular distributions of the neutrino radiation field for several spherically symmetric (1D) SN simulations has not found any crossing in the ELN near the neutrinosphere. More generally, 2D or 3D models can exhibit a large-scale dipole in the ELN emission, termed Lepton-Emission Self-sustained Asymmetry (LESA) [30], which also makes a crossing more likely to occur. In arXiv:1912.00274v2 [astro-ph.HE] 17 Jan 2020 this context, the first analyses of fast instabilities in multidimensional SN models have been recently performed in [31,32]. In [31] the authors extracted three snapshots from numerical data in 2D and 3D SN simulations and looked for ELN crossings in the angular distributions of ν e andν e . They found favorable conditions in extended regions with the radius of 50-70 km. Then, by a linear stability analysis of the neutrino equations of motion, they identified the strength of this instability for a representative point. However, their neutrino distributions were obtained from neutrino transport calculations done in a post-processing step, dropping the time dependence of both hydrodynamical and neutrino quantities and ignoring matter motions entirely. Hence they are not fully self-consistent. Conversely, in [32] the authors used fully self-consistent simulations in 2D, computing neutrino transport with a multi-angle Boltzmann solver coupled to hydrodynamics. Applying linear stability analysis near the neutrinosphere they found no positive signatures of conversion at least for the spatial points and times studied in ther particular model. Furthermore, in [33] the claim was made that in the pre-shock region, at r O(100) km, the residual coherent neutrino-nucleus scatterings could produce a tiny crossing in the ELN, whose presence has been confirmed by the inspection of various numerical simulations. Despite the smallness of the crossing, according to a stability analysis it would be enough to trigger significant fast conversions. However, for a cautioning argument against overinterpreting results of stability analyses, see [34].
Recently, two publications, based again on the SN models considered in [31,32], reported positive detections of locations of ELN crossings deep inside the protoneutron star (PNS) when investigating the selfconsistent 2D core-collapse simulation of an 11.2 M star computed with Boltzmann neutrino transport [35] and 3D Boltzmann neutrino-transport results for fixed matter background at some instants during the post-bounce evolution of 11.2 M and 27 M progenitors [36]. Because the diffusive conditions for neutrinos in the deep PNS interior imply that the angular distributions of both ν e and ν e are nearly isotropic, ELN crossings were found only in regions where the "asymmetry parameter" Γ = nν e /n νe , i.e. the ratio of the number densities of both neutrino types, is close to unity (see also [26]). Consequently and naturally in the equilibrium diffusion regime, the chemical potential of ν e nearly vanishes in these regions. The authors of [35] speculated that the appearance of light nuclei (among them α particles as the dominant species) is causal for the development of such instability conditions. Just as Delfan Azari et al. [35], also Abbar et al. [36] diagnosed ELN crossings in deep regions inside the PNS only in a small number of isolated points at the analysed post-bounce moments. They also correlated their occurrence with locations where the chemical potential of electron neutrinos nearly vanishes and pointed out that the electron fraction Y e is relatively low there and the temperature is close to maximal values.
These interesting results motivate the need to extend the search for fast instabilities to other state-of-theart and fully self-consistent multidimensional SN models. However, most multi-D SN simulations [30,[37][38][39][40][41][42][43][44][45][46] evolve only the (energy-dependent) angular integrals ("moments") of the neutrino phase-space distributions with time, and not the fully angle-dependent distributions. Reference [39], used in [32,35], is a welcome exception. The lack of detailed angular information seems to preclude a linear stability analysis that requires knowing these distributions. To overcome this limitation, some of us have recently proposed an alternative method to diagnose the possibility of fast instabilities in the absence of detailed knowledge of the ELN distributions [47]. This recipe is based on identifying a specific Fourier mode of the flavor instability field called the "zero mode", which has an easily calculable growth rate depending only on the angular moments of the ELN up to second order. It has been shown with numerical examples that the growth rate of this mode, calculated from the stability analysis, nicely approximates the growth of flavor conversions for the same mode in detailed numerical calculations.
The purpose of the present work is to use this new method of analysis to scan the different regions in selfconsistent, state-of-the-art 3D SN models with fully 3D two-moment neutrino transport for the possibility of fast flavor conversions therein. Specifically, we employ two time-dependent stellar core-collapse (and explosion) simulations for 9 M and 20 M progenitors recently published by the Garching group [45,46]. Both of these simulations exhibit the LESA phenomenon. We will demonstrate that the direct evaluation of discretized numerical data provided by computational models leads to the identification of only few, isolated points of ELN crossings. We will argue that instead of being such point-like locations, the regions of fast flavor instability are thin 2D layers that first appear around small 3D volumes in the convective layer of the PNS, and which grow with time as the convective and diffusive transport of electron-lepton number drives a decrease of Y e in the PNS convection layer. We also observe a strong hemispheric asymmetry of the thin layer of flavor instability correlated with the asymmetry of PNS convection leading to the LESA phenomenon. The regions of ELN crossings are much more extended in the hemisphere opposite to the LESA dipole direction, where PNS convection is weaker and Y e is lower. Our analysis thus shows that the locations of ELN crossings are not dot-like and fluctuating because of stochastic hydrodynamical variations, but they are longlasting and large-scale structures (for possibly important implications of this fact, see [34]).
In Sec. II, we describe our method for diagnosing instabilities based on the angular moments of the neutrino phase-space distributions. In Sec. III we present the results of our search for fast flavor instability in the two investigated 3D SN models, first by directly using the discretized output of the numerical simulations in the instability condition (Sec. III A), which leads to the iden-tification of only few isolated points of flavor instability located interior to the neutrinospheres in the convective shell of the PNS. In Sec. III B we discuss the conditions for ELN crossings. We argue that this direct analysis on the discrete numerical mesh fails to correctly identify the regions of flavor instability, which are actually thin 2D layers in 3D space. We propose an alternative, better strategy to find these layers containing the locations of flavor instability. In Sec. IV we present the time evolution of the instability layers in our two model runs (Sec. IV A) and explain the reason for the development of the relevant physical conditions (Sec. IV B). We also discuss their correlation with the dipole of the leptonnumber emission and the asymmetry of the electron distribution in the PNS convection layer connected with the LESA phenomenon. Finally, in Sec. V, we conclude with a brief summary.

A. Instability equation
In order to track the existence of the fast neutrino flavor instability one has to perform a linear stability analysis of the neutrino equations of motion. We refer the interested reader to [22] for a detailed discussion of this analysis and for the related publications. Here, we simply mention that in the growing literature about fast flavor conversions, a novel approach to study these effects was recently proposed in [17]. This is based on the dispersion relation for the frequency and wavenumber (ω, k) in the mean field of ν e ν x coherence, which is essentially the offdiagonal element of the neutrino density matrix (p, x, t) that we will call S in the following. One looks for solutions of the linearized equations for the flavor evolution in the form Typically such a solution may exist only if ω and k are related by an appropriate equation, called the dispersion relation. Loosely speaking, if either k or ω develop imaginary parts leading to positive real arguments in the exponential, the solution is expected to grow in space or time, thus signaling an "instability". Unfortunately, identifying the instabilities rests upon a more complicated analysis that requires a full characterization of the complex analytic structure of the dispersion relation [18].
In [47] some of us proposed a simpler analytical tool to diagnose the fast neutrino instability. Our proposal is based on identifying a specific Fourier mode of the flavor instability field that we call the "zero mode" in a corotating frame, in which one can gauge-away the matter term from the neutrino equations of motion. We labeled M = 9 M t = 300 ms R [12,15]  this zero-mode as k = 0. 1 This is motivated by the fact that the calculation of ω for this mode is significantly simpler than a full characterization of the dispersion relations, D(ω, k) [18]. In fact, for this mode the dispersion relation becomes with η µν = diag(+1, −1, −1, −1), i.e., D(ω, 0) is a polynomial in ω. The specific model of SN neutrino populations and their angular distributions, encoded in the ELN, only enters the equation through the tensor V µν (with µ, ν = 0, 1, 2, 3) that contains the angular moments of the neutrino distributions up to second order in the neutrino velocity, namely, where v µ = (1, p/E), i.e. the zeroth component of the velocity four-vector is 1 and the spatial components are given by the unit vector v = p/E (with E = |p|). The function is the difference of the phase-space occupation functions integrated over energy space, i.e., the angular distribution of the electron-neutrino lepton number (ELN) [15].
Here we assume that ν x andν x have identical distributions. Therefore, V µν depends on the angular moments of the neutrino-species dependent phase space distributions up to second order in the neutrino velocities, i.e., where the notation . . . να refers to These (energy-integrated) moments are related to the difference of the number densities (n ν ) of ν e andν e , the corresponding difference of the number-flux densities (F r ν ), and the difference of the second angular moments of the ν e andν e number distributions (P rr ν ) as follows: 1 A remark is in order. Our method does not completely exclude the presence of an instability for those points where we find Im(ω) = 0. Indeed, as our analysis is based only on the zero mode k = 0, there might be a k = 0 that is unstable. Therefore, one cannot exclude larger instability regions than those that we will show. The reader is referred to [48] for a comparison of instability regions obtained with different instability criteria.
where for later quantitative evaluation we have reintroduced the factors c in the expressions in the rhs of these relations.
We remind the reader that the two-moment "M1" neutrino-transport scheme used in [45,46] evolves the "00" and the "0i" components (with i ∈ {r, θ, φ} being radius r, polar angle θ, and azimuthal angle φ of the polar coordinate system) of the moments in Eq. (6) in time for all neutrino species. The system of moment equations of the transport solver is closed by an algebraic relation for non-evolved moments (i.e., for the "ij" components of the tensor with i ≥ 1 and j ≥ 1), which depend on the "00" and the "0i" components (see e.g. [49]).

A. Direct analysis of discretized numerical results
In this work we employ energy-integrated angular moments of the neutrino number distribution provided by the neutrino transport solver used in the considered 3D SN simulations of Refs. [45,46], Models s9.0 FMD H and s20 FMD H there. Appropriate normalization constants as specified in Eqs. (7)-(9) are applied in order to match the quantities in Eq. (3). We choose to work mainly in the comoving frame of the stellar fluid (fluid frame), where SN neutrino transport usually provides its output quantities. Equivalent results in terms of flavor instabilities are obtained in the laboratory frame, i.e. the rest frame of the stellar center. We will briefly demonstrate this later.
We mainly focus on our progenitor with 9 M and discuss similarities as well as differences compared to the 20 M simulation. Figure 1 displays Aitoff projections of the 9 M model for post-bounce times of t = 300 ms and 500 ms (upper two panels) and of the 20 M star for t = 200 ms after bounce (bottom panel). Different radial directions in the 3D simulations correspond to discretized values of the zenith angle θ and azimuthal angle φ. Points (φ, θ) are marked by green color if the solutions of Eq. (2) yield Im(ω) > 0 for at least one value of the discretized radial coordinate r in the range 10 km ≤ r ≤ 30 km. If the color is blue then Im(ω) = 0 for all radii, i.e. there is no instability.
In the 9 M model at 300 ms we find instability points only in the interval of [12,15] km and for very specific directions, i.e., a few isolated and unconnected pairs of values (φ, θ). As specified in the plots, Im(ω) ∼ O(10 − 100) m −1 at unstable locations, which means the flavor instability can develop over a time scale even shorter than a nanosecond. In the middle and lower panels of Figure 1 we witness a larger number of points of instability than in the upper panel.  for νe andνe as well as their differences ∆nν , ∆F r ν , and ∆P rr ν along a radial direction at (θ, φ) = (107 • , 57.4 • ) in our 9 M model at 300 ms after bounce. The four panels on the lower right display the corresponding flavor-instability functional F and, for comparison, the first term of it, 16 9 (∆nν ) 2 , as labeled in the panels.
Nonzero values of Im(ω) in the 9 M model occur only later than ∼300 ms after bounce, whereas in the 20 M the condition for flavor instability shows up earlier and more points with positive imaginary part of ω are present already at 200 ms.
We note in passing that we have not detected any points of flavor instability between the neutrinodecoupling region and the SN shock (at r ≥ 30-50 km, depending on the post-bounce time). Some authors have speculated about the possibility that ELN crossings might occur in this region in multi-dimensional models and, indeed, Ref. [31] as well as Refs. [36,50] have reported such locations in 2D and 3D simulations of an 11.2 M star. In this context, however, it is important to keep in mind that our moment-based criterion employs differences of angle and energy integrals of the neutrino distribution functions. Such an integral criterion is not necessarily sensitive enough to diagnose low-level cross-  ings in the angular distributions of ν e andν e , i.e., small differences in the distribution functions leading to a reversed ordering of the distributions in some narrow region or in a sparsely populated part of the angle space, as, e.g., spotted in the preshock domain as a consequence of a small fraction of neutrinos that are backscattered in collisions with infalling nuclei [33].
Our findings are reminiscent of the individual points of ELN crossings that were identified deep inside the PNS by [35] in the 2D simulation of this 11.2 M model and by [36] in the 3D simulations of the same 11.2 M progenitor and of a 27 M model. In the following section we will argue that the regions of fast-flavor instability form thin 2D layers in 3D space rather than isolated points, and the identification of individual points in our analysis is an artifact connected with the discretization of the physical variables in the numerical treatment. This conclusion does probably also apply to the results of the previous investigations in Refs. [35,36].

B. Conditions for ELN crossings
Similar to Refs. [35] and [36] our points of flavor instability are located deep inside the PNS, i.e. below the neutrinospheres in a region where neutrinos diffuse and their phase-space distributions are very close to those of local chemical equilibrium.
Before we discuss in detail why we -and probably also the references mentioned above-have found only isolated points instead of extended regions of flavor instability, we will introduce a simplified and approximative criterion for the instability, which makes the underlying physics more transparent. In order to achieve this, we take advantage of the fact that the nonradial moments are much smaller than the radial ones in the diffusive core where we find instability. Moreover, in this region V rr = V θθ = V φφ holds because the neutrino phasespace distributions are nearly isotropic. In this limit, which effectively corresponds to the 1D case, Eq. (2) is explicitly quadratic in ω , namely with the solution With the additional relation that V rr = 1 3 V 00 , which is valid to high accuracy in the diffusion regime, the condition for instability becomes In terms of the differences of number densities ∆n ν and radial fluxes ∆F r ν of ν e andν e (see Eqs. 7 and 8) this instability condition reads The presence of flavor instability is thus mostly dependent on V 00 and V 0r , i.e., on the differences of the neutrino number densities and number-flux densities of ν e andν e . We explicitly checked that Eq. (11) gives almost identical results to what is shown in Fig. 1.
A subtlety concerns the evaluation of the instability condition either with the moments in the reference frame comoving with the stellar fluid, where the neutrino quantities (i.e., the angular moments) are computed by the numerical transport code, or in the laboratory frame. In a full-angle treatment of the neutrino transport this corresponds to the question whether ELN crossings shall be searched for with the angular distributions of ν e andν e in the comoving frame or in the lab frame. We will show here that at the level of the angular moments, which we use for our analysis, the results are basically independent of the specific frame where the neutrino moments are evaluated. For fluid velocitiesṽ c, the lab-frame . The first small raisin-like volumes withνe excess signalled by negative ∆nν can be found at locations of particularly low Ye (intense blue hues for negative Ye fluctuations relative to the average value). The spatial variations of Ye are connected to lepton-rich convective updrafts, which carry electron-lepton number from the convectively stable PNS core outward, and more lepton-poor convective downdrafts. Larger red patches that mark buoyantly rising plasma in the bottom panels are therefore correlated with bigger orange regions in the panels of the third row, and more extended blue inward flows in the bottom panels coincide with the deepest-blue areas in the panels of the third row. and comoving frame moments are related to lowest order inṽ/c through In the frame transformations of n ν (Eq. 14) and P rr ν (Eq. 16) we omit terms such as 1 c 2ṽi F i ν and 1 c 2ṽ r F r ν , respectively. Sinceṽ c holds and in the diffusion regime interior to the neutrinospheres also 1 c |F i ν | n ν r = 14 km applies 2 , the disregarded terms are many orders of magnitude smaller than the leading ones that we retain in 2 This can be easily verified by comparing the left panels in the top and third rows of Fig. 2 for r 30 km and taking into account that |F θ ν | ≈ |F φ ν | |F r ν |.
Sinceṽ r c and, as we shall argue below, the condition can be fulfilled only when ∆n ν ≈ ∆n lab ν ≈ 0, the relation in Eq. (17) is basically identical with which is identical to the instability condition of Eq. (13).  , the corresponding second angular moments P rr ν and their difference (four upper right panels), the radial neutrino-flux densities F r ν and their difference (four lower left panels) and the "flavorinstability functional" F of Eqs. (13) and (18) (four lower right panels). The angular direction (θ, φ) for the radial ray was chosen such that one of the instability points visible in the top plot of Fig. 1 was crossed. This can be seen in the four panels on the lower right of Fig. 2, where at r ≈ 14 km the flavor-instability condition is fulfilled. The four upper left panels demonstrate that at this location n νe and nν e are approximately equal. Lab-frame and comoving-frame quantities exhibit exactly the same behavior.
A comparison of the four upper left and four upper right panels shows that the same conclusion can be drawn from inspection of P rr ν , because in the diffusion region P rr ν = 1 3 n ν is very well fulfilled. This relation does not hold any longer when neutrinos begin to decouple from the stellar medium near the neutrinosphere and undergo the transition to free streaming outside. In this case P rr ν → n ν asymptotically for r → ∞, and therefore our flavor-instability conditions of Eqs. (13) and (18) are not valid any more. In the displayed model this is the case for radii r 30 km, for which reason the negative values of the flavor-instability functional for r 40 km do not signal flavor instability in this region exterior to the PNS.
The two lower right bottom panels of Fig. 2 also display the term 16 9 (∆n ν ) 2 as part of the flavor-instability functional for comparison with the full expression. One can see that this term usually dominates the second one, 4 c 2 (∆F r ν ) 2 , by several (typically by 2-3) orders of magnitude. This can also be directly verified by comparing ∆n ν in the upper left panels with 1 c ∆F r ν displayed in the lower left panels. We remark in passing that strongly negative values of theν e flux in the comoving frame occur because of a local temperature maximum that drives the diffusion flux ofν e inward while the more degeneracydriven diffusion flux of ν e can still be outward directed. Although the lab-frame and comoving-frame fluxes are considerably different (because the advective component v r n ν can dominate the diffusive component in the convection layer of the PNS), the radial profiles of ∆F r ν are more similar for lab-frame and comoving-frame fluxes, and the instability functional F in the lower right panels does not exhibit any visible frame dependence.
There are severe consequences of this huge imbalance between the first and the second term in the flavorinstability functional F when searching for ELN crossing points by evaluating the functional with discretized numerical results. In order to detect such points, i.e. in order to find grid locations where F < 0, the term 16 9 (∆n ν ) 2 must be very close to zero at exactly such grid positions, because only then the small second term can lead to a negative value of F. If, however, the discrete grid points are too far away from the root of F, the val-ues of 16 9 (∆n ν ) 2 at these points may be so large that the second term 4 c 2 (∆F r ν ) 2 does not achieve to produce negative values of F. This, in fact, is likely to happen in the far majority of all cases where the physical conditions enable flavor instability, and only in a minor fraction of such locations the discretized spatial points of the computational grid coincide incidentally with locations where the combination of terms can yield F < 0, and thus signal the presence of flavor-unstable conditions.
An example of such a missed point of instability can also be spotted in Fig. 2. The plots of F in the lower right panels show two local minima between 10 km and 20 km. Only in the case of the left one the minimum of F reaches a negative value, but for the right one the minimum is still on the positive side. Is the condition of flavor instability also fulfilled at this position and just not detected by the numerical analysis? Indeed, this is the situation as visualized in detail by a close-up of the region of the two local minima of F in Fig. 3 with crosses marking the positions of all radial mesh points of the computational grid. The orange line in the two panels corresponds to the radial direction chosen for the profiles in Fig. 2. It is obvious that ∆n ν (upper panel of Fig. 3) has two zero crossings and becomes negative between these two roots. If a mesh point happens to be close to the root, the small, second term in Eqs. (13) and (18) achieves to drive F (lower panel) to the negative side. Such a situation occurs for the left one of the two roots along the radial direction at (θ, φ) = (107 • , 57.4 • ) and for the right root in the case of (θ, φ) = (111 • , 57.4 • ). If, however, the mesh points are too far away from the root, then F remains positive at all discrete points of the grid. This is the situation for the direction corresponding to (θ, φ) = (109 • , 57.4 • ) displayed in Fig. 3, although also in this case ∆n ν possesses two roots. The fourth selected case in this figure for (θ, φ) = (105 • , 57.4 • ) does not exhibit any change in the sign of ∆n ν .
Discretization effects are therefore the reason why only very few points with ELN crossings could be identified by the analysis so far. Consequently, only individual, isolated points of instability appeared on the three panels of Fig. 1. This problem would become even more severe if the resolution of the computational grid used in our SN simulations had been coarser.
Our case is only an example how numerical discretization effects may impede the possibility to detect flavorunstable conditions that occur in narrowly delimited spatial regions. A specific condition such as, e.g., our instability criteria of Eqs. (2), (13), and (18) or any other analytical criterion relating physical quantities that are available on a discrete numerical mesh, may fail to identify the spatial locations of instability. We speculate that also the investigations in Refs. [35] and [36], determining ELN crossings by using angular distributions, suffered from the finite numerical resolution of the underlying SN simulations and therefore failed to find more spatial points of flavor-unstable conditions.
Actually, the problems encountered in our analysis with discretized physical variables can easily be circumvented. Evaluating the relations of Eqs. (13) or (18) to search for spatial locations where F < 0 (or for roots of F) on a mesh of discrete points is not a promising strategy. Instead, it is preferable to look for regions where ∆n ν changes its sign, which we understood as a sufficient condition to obtain roots of F. 3 When ∆n ν changes its sign between two grid points, there must be a root of this quantity between the two points. Close to this root F will become negative, unless ∆F r ν vanishes in this region. In such a pathological and very rare situation our approximate flavor-instability criterion based on a few angular moments of the neutrino phase-space distributions does not provide conclusive information.

A. Time evolution
With the approximative but numerically robust criterion of sign changes of ∆n ν , we have evaluated our 9 M and 20 M models in time to track the evolution of the volume of ELN crossings interior to the PNS in our 9 M and 20 M simulations.
For the 9 M model we show ∆n ν and the normalized fluctuations of the electron fraction, (Y e − Y e )/ Y e (the angle brackets indicate averages over zenith and azimuthal angles), in full-sphere Aitoff projections at a radius of r = 14 km as well as cross-sectional cuts in the x-z, y-z and x-y planes at post-bounce times of 300, 400, 500, and 600 ms in Figs. 4, 5, 6, and 7, respectively.
The fluctuations of Y e are connected to convective updrafts and downdrafts in the convection layer of the PNS (see Refs. [30,45]). Convective updrafts carry electronlepton number from the convectively stable high-density core of the PNS outward and therefore exhibit higher values of Y e than the angle average. In contrast, convective downdrafts are more lepton-poor, which also means that they contain more neutrons, which makes them specifically heavier so that they sink inward. The pattern of Y e fluctuations mirrors the familiar cell pattern of convection in spherical shells.
The zero-crossings of ∆n ν , and thus the locations very close to the flavor instability, are highlighted by yellow lines surrounding the volumes of negative values colored 3 Strictly speaking, a sign change of ∆nν is not necessary to get F < 0, but this condition for F can also be fulfilled when ∆nν dips nearly to zero while still remaining positive. However, because ∆nν and 1 c ∆F r ν are orders of magnitude different in the diffusion regime (compare left panels in the second and fourth rows of Fig. 2), such a situation is highly fine-tuned and not as common in hot PNSs as sign changes of ∆nν . This is obvious from our analysis of the time evolution of PNSs in two 3D SN simulations. . The development of a pronounced LESA dipole at t > 300 ms is obvious, the dipole direction is indicated by black asterisks and is always close to the +y axis of the computational polar grid (see Model s9.0 FMD H in [45]).
in blue. The physical thickness of these boundary layers of the ∆n ν < 0 volumes, i.e., the "skins" in which the flavor-instability condition is fulfilled, can be roughly estimated from Eqs. (13) or (18) by making use of the diffusion approximation to express the lepton-number flux, ∆F r ν = F r νe − F r νe through F r νi = −D νi ∂n νi /∂r, where D νi = 1 3 cλ νi is the diffusion coefficient and λ νi the (energy-averaged) mean free path. Introducing a mean free pathλ suitably averaged between ν e andν e , we can write for the effective neutrino-lepton number flux in the diffusion regime: The requirement for flavor instability, F < 0 (Eqs. 13 and 18), implies that |∆F r ν | > 2 3 |∆n ν | c, or, using Eq. (19), that the radial scale height of changes of the leptonnumber density must be smaller than 1 2λ : This means that the thickness of the skins of flavor instability is only a fraction of the ν e -ν e -averaged mean free pathλ, which is on the order of meters in the PNS region of relevance. This skin thickness is one to two orders of magnitude below the spatial resolution of the best 2D and 3D simulations, which explains why a direct evaluation of the flavor-instability criterion on the discrete grid points of the computational mesh can find, basically incidentally, only very few locations of ELN crossings. Around 300 ms first spots of flavor-instability can be seen near the radius of 14 km. However, while some moments earlier a small number of individual, isolated points may have fulfilled the instability condition F < 0, at 300 ms these points have already grown to 2D surfaces enclosing noticeable volumes where ∆n ν < 0. Such raisin-like inclusions are concentrated around regions where Y e is 10-15% lower than the average. At 400 ms the ∆n ν < 0 volumes have considerably grown and partly merged, enveloped by a coherent surface, besides still existing smaller droplets. This trend continues until our last displayed snapshot at 600 ms. It is obvious that by this time the quadrupole-dominated pattern of concentrations of instability regions that characterizes the situation at 400 ms and 500 ms has evolved to a distribution that is clearly dominated by a prominent dipolar asymmetry. While in one hemisphere there are extended blue 3D regions with ∆n ν < 0, the opposite hemisphere still exhibits only scattered spots where this condition is fulfilled. As time goes on these blue regions grow along with a decreasing electron fraction and rising density because the PNS gradually deleptonizes and contracts. It is obvious that the large-scale quadrupolar and dipolar asymmetries of these regions correlate with such asym-  [45]); its vector direction is marked by white asterisks in the Aitoff projections and black arrows in the cuts. Within only 50 ms initial "raisins" with ∆nν < 0 (blue) grow to large "pancakes" within a shell around r = 14 km. At 270 ms nearly the whole shell is included with a few remaining "holes" of ∆nν > 0 (red) remaining in the hemisphere of the LESA vector direction. At 370 ms the shell has become wider with no remaining holes. Yellow lines indicate the flavor-unstable locations at the boundaries between regions of ∆nν < 0 and ∆nν > 0. metries in the relative variations of Y e . Figure 8 visualizes by 3D volume rendering the situation in the 9 M model at the post-bounce time of t = 300 ms, when the first scattered "raisins" with ∆n ν < 0 inside and flavor-unstable conditions in their skins have grown. This is compared to the situation at 600 ms when a whole 3D shell between ∼10 km and ∼14 km withν e excess over ν e has developed, more widespread on one side of the PNS than on the other one. This hemispheric asymmetry is connected to the LESA phenomenon and establishes in correlation with the growing dipole amplitude of the LESA, which can be seen in the Aitoff projections of the normalized lepton-number flux, ∆F ν / ∆F ν = (F νe − Fν e )/ ∆F ν (evaluated far outside the PNS at r = 400 km; Fig. 9), and of the Y e distribution inside the PNS at r = 14 km (Figs. 4-7).
The black asterisks in the Aitoff plots in all of these figures and the black arrows in the cross-sectional cuts of Figs. 4-7 mark the dipole directions of the leptonnumber flux. These markers are missing in the plots for t = 300 ms, because at that early time the emission dipole is not well developed. But once it is present in a clear way, the dipole direction is extremely stable in the 9 M model, varying only very little around the +y-direction of the polar coordinate grid of the computation (see Model s9.0 FMD H in figure 3 of [45]). It is evident from all these figures that ELN crossings are less favored on the hemisphere in the direction of the LESA dipole vector. As we will discuss in more detail in Sec. IV B, the reason for this observation is the higher electron fraction in the PNS convective shell on this side. As discussed in [45] (see also [30]), convection inside the PNS is stronger in the hemisphere of the LESA dipole direction. This stronger convection transports electron-lepton number more efficiently out from the edge of the nonconvective central core of the PNS, thus raising the Y e in the convective shell as well as in the overlying near-surface layers of the PNS up to the neutrinospheres. Figure 10 presents the evolution of the flavorinstability regions in our 20 M model by color-coded Aitoff projections and cross-sectional cuts of ∆n ν in analogy to Figs. 4-7. The time sequence shows the same basic features as visible in the 9 M model, just accelerated because the PNS grows faster in mass due to the higher mass-infall rate in the more massive progenitor. Spots with flavor-unstable conditions are visible already at 170 ms after bounce near a radius of 14 km and have grown to large coherent structures at 220 ms. At 270 ms the regions of ∆n ν < 0 dominate and form a coherent shell around r = 14 km. The very few remaining volumes with ∆n ν > 0 in this shell have essentially disappeared until 370 ms after bounce. The direction of the LESA emission dipole is indicated by white asterisks in the Aitoff projections and black arrows in the cross-sectional cuts at the times when a strong dipole of the lepton-number emission exists. Again one can notice that flavor-unstable conditions are more widespread in the anti-LESA direction, where the electron fraction is lower than in the hemisphere the LESA dipole vector is pointing to. Accordingly, the last islands with ν e excess (i.e., ∆n ν > 0) can be found on this side of the PNS (see panels for t = 270 ms in Fig. 10).
In the following section we will explain the physical connection between the conditions of flavor instability and the evolution of electron fraction Y e , density ρ, and temperature inside the PNS in more detail.

B. Properties of ELN crossing points
Naturally, the relation ∆n ν = n νe − nν e ≈ 0, which we recognized in Sec. III B as a necessary condition for ELN crossings in a regime where neutrinos diffuse and are close to local chemical equilibrium, implies that µ νe = µ e + µ p − µ n ≈ 0 holds for the electron-neutrino chemical as a function of the chemical potentials of electrons, protons, and neutrons. Only then ν e andν e in local chemical equilibrium can have nearly equal number densities. This was pointed out already in Refs. [35] and [36].
But what is the reason why fast flavor unstable conditions develop in an increasing volume in the new-born NS? The relation of the chemical potentials that determines µ νe suggests that low Y e and correspondingly low µ e might allow µ νe to drop to zero and even to negative values. Figure 11 displays radial profiles of these and other quantities in the LESA and anti-LESA directions of the 9 M model at four post-bounce times. Inspecting the different panels shows that the locations where fast flavor instability is possible (∆n ν ≈ 0) lie within the convective shell in the interior of the PNS; the convective layer can be recognized by the convectively flattened gradient of the entropy per nucleon, s. However, neither Y e nor µ e reach a local minimum where the electron-neutrino chemical potential develops a trough near zero or below, and also the abundances of α particles or light elements up to 4 He are not particularly abundant in general in the regions where ELN crossings are favored, in contradiction to arguments in [35]. This can be directly concluded from a comparison of the LESA and anti-LESA directions in Fig. 11. In the hemisphere opposite to the LESA dipole vector, Y e as well as the abundances of α particles and light elements are lower than in the hemisphere the LESA vector points to. Nevertheless, flavor unstable locations are more widespread in the anti-LESA direction at late times (500, 600 ms after bounce) when the dipole is strong (as discussed in Sec. IV A). Lower values of Y e enhance the possibility of flavor instability, but they are not causal in the first place.
Instead of being linked to peculiarities in the Y e profile, the µ νe trough is found in a region whereμ ≡ µ n − µ p possesses a local maximum, because high values ofμ reduce µ νe = µ e −μ compared to µ e . The local maximum ofμ again correlates tightly with a local temperature maximum, which is a relic of shock heating by the initial propagation of the strong shock front formed at the moment of core bounce. The close connection between maxima ofμ and T can be easily understood from considering nondegenerate, nonrelativistic and noninteracting neutrons and protons as an ap-proximation of nucleons in the subnuclear regime. For such Boltzmann gases the chemical potential of particle species i = n, p with number density n i and rest-mass m i is given by µ i = m i c 2 + k B T ln(Λ 3 th,i n i /2), where The crosses indicate the conditions that are found inside the raisin-like volumes with ∆nν < 0 at these early times. Some of these conditions tend to be extreme compared to the angular averages, for which reason the crosses do not lie on the curves. Λ th,i = h (2πm i k B T ) −1/2 is the thermal wavelength, k B the Boltzmann constant, and h the Planck constant.
With Q ≡ (m n − m p )c 2 one therefore gets: This relation explains the direct dependence ofμ on the plasma temperature, and it also implies that lower Y e in the hemisphere opposite to the LESA vector reduce µ νe not only through lower values of µ e but also through higher ratios of n n /n p . 4 Non-negligible effects due to nucleon interactions in the density regime of interest between some 10 13 g cm −3 and ∼10 14 g cm −3 may lead to quantitative changes but do not qualitatively affect this argument. 4 We note that the logarithmic dependence ofμ on nn/np is weak, but also µe depends logarithmically on ne when electrons are nondegenerate, and µe is proportional to n 1/3 e when electrons are extremely degenerate.
The development of regions with negative electronneutrino chemical potential, µ νe < 0, inside the PNS does not only depend on temperature and Y e but also on density. This can be concluded from Fig. 12, where the evolution of the 9 M and 20 M models is compared. The upper panel displays the growth of the volume of regions with µ νe < 0. Vertical black lines mark the two earliest instants displayed for both simulations in Fig. 4 and Fig. 10, respectively. The solid lines represent angular averages of the quantities over the whole sphere at r = 14 km, and the crosses indicate the conditions inside the (still very small) volumes that fulfill ∆n ν < 0 at 300 ms (170 ms) in the 9 M (20 M ) model (the two instants marked by the vertical black lines). Since the conditions in the µ νe < 0 volumes are special and for some quantities tend to be extreme compared to the conditions on the rest of the sphere, the crosses do not lie on the curves of the angular averages. 5 The crosses in the second panel of Fig. 12 show that at the marked instants the electron-neutrino chemical potential begins to become negative in the raisin-like spots that are visible in the Aitoff projections of the two models (Figs. 4 and 10), whereas the angle-averaged values of µ νe at r = 14 km are still somewhat higher, with a monotonically decreasing trend with time. The remaining three panels contain the evolution of density, electron fraction, and temperature, again angle-averaged at r = 14 km. In both models flavor-unstable raisin-skins begin to occur when the density reaches roughly 5 × 10 13 g cm −3 , whereas Y e and temperature are different between the two models at this time.
From Fig. 13 it becomes clear where these differences originate from. The plot shows a set of isocontours for different, fixed values of the electron fraction (colorcoded) in the density-temperature plane. All contours are defined by the condition µ νe = 0. At each point (ρ, T ) this condition is fulfilled only for a single value of the electron fraction. Each µ νe = 0 contour encloses a ρ-T -domain where µ νe > 0 for the Y e value of the contour, whereas µ νe < 0 holds for the same Y e value outside of this µ νe = 0 contour (see also figures 11 and 12 in Ref. [52]).
In the inset of Fig. 13, the evolution tracks of both SN models in the ρ-T -Y e -space, with all quantities angleaveraged at r = 14 km, are superimposed on the contours of vanishing neutrino chemical potential. The tracks start at low values of temperature and density with a high value of the electron fraction, and evolve from the lower left to the upper right of the inset as the PNS contracts and heats up by compression and conversion of electron degeneracy energy to thermal energy (Joule heating). ELN crossings become possible when the electron fraction on the track matches the Y e value of one of the contours, i.e. when the color of the track and of a contour are the same. The symbols mark the (ρ, T ) locations when first noticeable flavor-unstable spots occur in the two SN runs at r = 14 km (corresponding to the two vertical lines in Fig. 12).
In the 20 M model the mass of the PNS grows faster than in the 9 M simulation because of the bigger mass-accretion rate in the more massive progenitor. Therefore the temperature in the PNS interior increases more rapidly and more steeply with density, and reaches greater values, allowing for a match of the electron fraction with one of the µ νe = 0 isocontours more quickly and at a higher value of Y e . As a consequence, fast flavor unstable conditions occur earlier in the 20 M model than in the 9 M case.
Compressional and Joule heating and continuous deleptonization are characteristic features of the neutrino-cooling evolution of new-born NSs. Our study includes a low-mass NS of a 9 M progenitor as well as a more massive NS in the 20 M model. Both of them develop flavor-unstable conditions in an increasing volume of the convective layer (Fig. 12), and the rising temperature naturally leads to a match of the decreasing Y e values along evolution tracks with the Y e of isocontours for µ νe = 0 (Fig. 13). We therefore expect that fast flavor unstable regions in the deep interior of PNSs are a generic phenomenon during PNS cooling.

V. CONCLUSIONS
We performed a detailed investigation of 3D state-ofthe-art SN models for the presence of fast neutrino flavor instability and to study the favorable conditions. These fast conversions are associated with crossings in the angular distribution of the ELN. However, the simulations we analyzed did not provide the detailed neutrino angular distributions. In order to overcome this limitation, we adopted a novel method for calculating the growth rate of the instability, which is based only on the angular moments of the ELN up to second order [47]. We applied our analysis to SN simulations of 9 M and 20 M progenitors, which were recently conducted by the Garching group [45,46]. In both models we found conditions for fast flavor instability deep inside the PNS in a radial range of 10 km r 20 km, where all flavors of neutrinos are in the diffusive regime and close to local chemical equilibrium.
We thus confirm similar recent detections of locations of ELN crossings in the PNS interior based on Boltzmann transport results in 2D time-dependent and 2D/3D fixed-background SN models in Refs. [35,36]. However, we showed that the direct evaluation of flavor-instability conditions with the discretized output of numerical simulations leads to the identification of only few individual, isolated points of ELN crossings. We argued that this finding on grounds of our method, and most likely also by the approaches used in the previous analyses, is a nu-merical artifact and misses the full phenomenology of the physical effects.
In reality the spatial locations of fast flavor instability are extended, narrow boundary layers surrounding volumes in 3D space where theν e number density exceeds the ν e number density. These surface layers have a thickness of roughly a neutrino mean free path and they grow from the skins of initially scattered, raisin-like inclusions to the envelopes of increasingly larger regions, until they finally form the inner and outer surfaces of a closed layer with nν e > n νe .
The region where this happens is located within the convective layer of the PNS where on the one hand convective transport of lepton number causes a rapid decline of the electron fraction, and on the other hand a local temperature maximum is further increased when the PNS contracts and compression as well as the conversion of electron-degeneracy energy to thermal energy heat the stellar plasma. Rising density and temperature combined with a decreasing electron fraction naturally drive this layer in the PNS towards conditions where the electron-neutrino chemical potential is µ νe ≤ 0, implying an excess ofν e relative to ν e . Since the lepton-number emission dipole associated with the LESA phenomenon is caused by stronger PNS convection in the hemisphere of the dipole direction, conditions for ELN crossings are more widespread in the opposite hemisphere, where PNS convection is weaker and therefore the electron fraction is not efficiently replenished by leptons carried outward from the edge of the nonconvective central core.
Our result backs the findings in [35,36] and points to ELN crossings in the PNS convection layer as a generic phenomenon during the cooling evolution of new-born NSs. In contrast to these previous works we have demonstrated that the regions of fast flavor unstable conditions are not fluctuating and point-like, but instead they are large-scale and long-lasting spatial structures.
This opens new directions of research. Presently, however, it is unclear whether fast flavor conversions in the deep interior of the PNS can have major consequences for the PNS cooling and/or SN evolution. The instability takes place in an extremely thin layer where µ νe is close to zero. Since these locations are deep inside the neutrinospheres of all neutrino species, neutrinos of all flavors are very close to chemical equilibrium. Therefore µ νe ≈ 0 implies that ν e andν e possess phase-space distributions that are very similar not only to each other but also to those of muon and tau neutrinos. At such conditions flavor exchange might have little impact on the overall conditions. Given the importance of potential effects, our work should stimulate further investigations. In view of the mutual support between our results and those obtained in Refs. [35,36], we conclude that our method based on angular moments is well suitable to analyze SN models computed with neutrino transport schemes that do not provide the detailed neutrino angular distributions.
The presence of fast conversions adds an additional layer of complexity to SN simulations. The current paradigm of flavor conversions is based on a separation from the treatment of neutrino interactions and transport. However, if fast conversions occur in the region of neutrino spectra formation or even deep inside the PNS, a self-consistent characterization of both of these phenomena will require new strategies for a simultaneous treatment of flavor oscillations and collisional effects. This will be a formidable task, since fast flavor conversions occur on time and length scales much shorter than usually resolved in global simulations of stellar collapse and explosions. Therefore, new methodical approaches will be needed. It is obvious that a huge amount of work still remains to be done to understand the role of neutrino flavor conversions in the cores of collapsing stars.