Fast neutrino flavor conversions in one-dimensional core-collapse supernova models with and without muon creation

In very dense environments, neutrinos can undergo fast flavor conversions on scales as short as a few centimeters provided that the angular distribution of the neutrino lepton number crosses zero. This work presents the first attempt to establish whether the non-negligible creation of muons in the core of supernovae can affect the occurrence of such crossings. For this purpose we employ state-of-the-art one-dimensional core-collapse supernova simulations, both in the presence and absence of muons creation. Although muon creation does not seem to modify significantly the conditions for the occurrence of fast modes, it allows for the existence of an interesting phenomenon, namely fast instabilities in the $\mu-\tau$ sector. We also show that crossings below the supernova shock are a relatively generic feature of one-dimensional simulations, which contrasts with the previous reports in the literature. Our results highlight the importance of multi-dimensional simulations with muon creation, where our results must be tested in the future.


I. INTRODUCTION
At the end of their lives, massive stars (M 8M ) may undergo a violent explosion, labeled as core-collapse supernova (CCSN). During such an event, about 10 53 ergs of energy is released into neutrinos of all flavors, corresponding to about 99% of the total released gravitational energy. Such neutrinos are thought to be a crucial ingredient of the astrophysical processes leading to the final explosion. Indeed, due to the extremely high density below the stalled shock, neutrinos might deposit enough energy to revive the shock, thus leading to the explosion. This is also known as the neutrino heating mechanism [1][2][3]. Neutrino interactions also affect the proton fraction Y e , which may indicate the presence of favorable conditions for the nucleosynthesis of heavy nuclei through r-processes, if its value is less than 0.5. Consequently, SN neutrinos represent an invaluable source of information of the processes occurring in the exploding star.
Neutrinos can experience flavor conversions during their propagation in the SN environment. In particular, in the highest density SN regions the neutrino selfinteraction potential becomes large enough to induce the so-called collective conversions [4][5][6][7], where neutrinos transform their flavors in a coherent fashion, despite having different oscillation frequencies. Such phenomena have a deeply non-linear nature and currently are not completely understood. This represents one of the most important unsolved problems in neutrino physics to date, since it impacts what we can learn from the next SN neutrino burst.
Lately, a lot of attention has been devoted to the so called fast flavor conversions [8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25], which are triggered by a crossing between the angular distributions of neutrinos and antineutrinos, i.e., a crossing in the angular distribution of the neutrino electron lepton number (ELN). In other words, such crossings occur when there is an excess of neutrinos in a given range of propagation directions, whereas antineutrinos dominate in another one. If such a condition is satisfied just above the neutrino trapping region, fast conversions will develop on time scales as short as a few nanoseconds and might lead to a significant modification of the flavor content of the neutrino fluxes. Recent studies have looked for crossings in a plethora of 1D and multi-dimensional numerical simulations outputs [26][27][28][29][30] and found them in the following regions: inside the proto-neutron star (PNS) [28][29][30], in the neutrino decoupling region [27] and in the pre-shock region [31]. Each of these zones rely on different microphysics for crossing generations.
The previous findings are all based on simulations adopting the standard assumption that the number densities and energy spectra of ν µ ,ν µ , ν τ ,ν τ are equal. This is an effective approach where only three different neutrino species exist: ν e ,ν e and ν x , where ν x = ν µ = ν τ = ν µ =ν τ . This assumption may only be taken in the absence of muon creation in the core of a SN. However, the temperature reached inside the core is usually enough for producing a non-negligible population of muons [32][33][34][35]. When this is taken into account, one finds that both ν µ andν µ production are enhanced, especially for the latter. This happens because the highly degenerate Fermi sea of e − is partially converted to µ − . Because the muon lepton number is conserved and initially equal to zero, this leads to an increasedν µ production. Furthermore, the creation of µ − and and µ + effectively softens the equation of state of the matter inside the core by converting a fraction of the thermal and degeneracy energy of electrons into rest-mass energy of muons. This leads to an accelerated shrinking of the PNS and therefore, to higher average energies of electron neutrinos, with a consequent increase of neutrino heating. The net result is a higher chance for a successful explosion compared to the case without any muon creation [33,36].
The goal of our work is to assess the potential impact of muon creation on the presence of ELN crossings. More specifically, we use two sets of data: • D1, provided by one-dimensional simulations discussed in Ref. [34] for a 18.6 and 20M progenitors (adopting the SFHo equation of state [37,38]), both with and without muon creation.
In these two data sets we look for the presence of ELN crossings through the method proposed in Ref. [40], which is based on the moments of the neutrino angular distributions. In D1 we find that ELN crossings are relatively common in one-dimensional SN simulations and they occur both below (post-shock region) and above (pre-shock region) the SN shock. This happens regardless of the inclusion of muon creation, though its inclusion makes crossings more numerous, or equivalently easier to capture with the moments. Furthermore, the inclusion of muon creation interestingly allows for a number of ELN crossings in the µτ sector, whereas none are otherwise expected. In D2 we do not find crossings in the post-shock SN region. We believe that such a discrepancy with the results of D1 is due to non equal theoretical inputs in the corresponding numerical models.
In D2 the neutrino angular distributions are also available so one can look directly at these for a crossing search. Using this unique opportunity we confirm the results obtained for D1 with some additional insights.
The structure of the paper is as follows. In Section II we briefly describe the SN models adopted for D1 and D2. In Section III we review the mathematical framework proposed in Ref. [40] to assess the presence of crossings by just using the moments of the neutrino angular distributions. Here we highlight how this method is modified when including muon creation. In Section IV we apply the method previously described to both D1 and D2 and we provide a list of original findings. In Section V we provide our interpretation of the results obtained in D1, D2, from their comparison and we finally draw our conclusions. When appropriate, we also comment on the results obtained in a similar analysis performed in Ref. [31].

A. D1
This data contains simulations for two progenitor masses: 18.6 and 20M , already presented in Ref. [34] 1 .
The progenitor models are taken from [41] for the 18.6M case and from [42] for the 20M one. The simulations are performed in 1D with the PrometheusVertex code with general-relativistic corrections [43] and six-species neutrino transport, which is solved iteratively through two-moment equations for neutrino energy and momentum, where a Boltzmann closure [44] is adopted. The effect of redshift and time dilation are also included in the transport. The tangent-ray method assumes neutrinos to propagate along straight paths instead of curved geodesics, since the effects of ray bending are relatively small. We consider the following set of neutrino reactions (see Appendix A of Ref. [45] for more details): neutralcurrent and charged-current interactions with nucleons including nucleon recoils and thermal motions, nucleonnucleon correlations, weak magnetism, a reduced effective nucleon mass and quenching of the axial-vector coupling at high densities; nucleon-nucleon bremsstrahlung; νν scattering; and neutrino-antineutrino-pair conversions between different flavors [46]. In neutrino absorption on nucleons we also take into account the mean-field corrections for the nucleons consistently with the equation of state. Electron capture and inelastic neutrino scattering on nuclei are also included [47,48]. Furthermore, neutrino interactions involving muons are taken into account, as described in Ref. [33]. A mixing-length approach is employed to simulate PNS convection [5]. Explosions are artificially triggered for the 18.6M model, a few 100 ms after bounce, as described in Ref. [5]. On the other hand the 20M model is not exploding. For all progenitor masses simulations are performed with the SFHo [38] equation of state. This equation of state is fully compatible with all current constraints from nuclear theory and experiment [49][50][51] and astrophysics, including pulsar mass measurements [52][53][54] and the radius constraints deduced from gravitational-wave and Neutron Star Interior Composition Explorer measurements [55][56][57]. To provide a better comparison with D2 (see next subsection) the simulation of the 18.6M model is performed also with the LS220 [58] equation of state.

B. D2
This data contains models for three progenitor masses of 11.2, 15 and 25M 2 . The progenitor models are taken from [59] for the 15M case and from [41] for all other cases. The nuclear equation of state is LS220 from [58] with compressibility modulus K = 220 MeV [60]. The simulations were performed with the 1D version of the angular_moments/, where the angular moments can be extracted.
Radial profiles are available upon request at https://wwwmpa.mpa-garching.mpg.de/ccsnarchive/data/ Bollig2016_radial_profiles/. same (Prometheus-Vertex) code employed for D1. Therefore the method for solving the hydrodynamical evolution and neutrino transport is the same as for model set D1, as well as the set of neutrino interactions accounted for 3 . The main differences with D1 are the following: here only three-neutrino species are considered (ν e ,ν e and ν x ) and neither muon production nor proto-NS convection are included. Moreover, no artificial explosion is triggered.

III. ANALYSIS METHOD
Under the assumption of ν x =ν x and azimuthal symmetry in the neutrino gas, an ELN crossing occurs when changes its sign at a given γ 0 . Here, γ = cos θ with θ being the neutrino propagation angle with respect to the radial direction and f να is the neutrino occupation number for a given flavor α. The quantity G e physically represents the angular distribution of the ELN. However, the inclusion of muon creation inevitably introduces a non-negligible difference between ν µ andν µ . In this case, the relevant quantities to be taken into account are [22] where G µ and G τ are defined equivalently to G e in Eq. (1), but for ν µ and ν τ , respectively. Analogously to G e , G αβ represents the angular distribution of the difference between the lepton numbers for flavor α and β. Considering that in general f νµ < fν µ , some differences with the standard calculations in the absence of muons should be expected. Nevertheless, whether the quantities in Eqs.
(2)-(4) can change sign has never been studied so far, mainly because the differences between ν µ andν µ are usually thought to be not comparable to those between ν e andν e . Moreover, until very recently no CCSN simulations has included muon creation. The first two-dimensional simulations in this direction have been presented in Ref. [33]. Usually, multidimensional simulations of CCSNe do not provide full angular information of G αβ (γ), since the equations for neutrino transport are solved for a number of the moments of G αβ (γ), which are defined as 3 In the charged-current neutrino-nucleon interactions, the simulations in D2 did not yet consider the nucleon mean-field potentials.
In Ref. [40], a simple method was proposed to assess the presence of crossings in the angular distribution of G αβ (θ). This method attempts to find a positive function F(γ) > 0, such that the quantity satisfies the following relation If such a function exists, then G αβ (γ) has a crossing. In order to write I αβ F in terms of the available moments I αβ n , defined in Eq. (5), we choose F(γ) to be a polynomial of γ where N is the highest moment available from numerical simulations, and usually is not larger than 3, and a i 's are an arbitrary set of coefficients for which F(γ) > 0. According to this choice, I αβ F becomes

IV. ANALYSIS RESULTS
Apart from observing crossings in the µ − τ sectors, we find a number of crossings below the SN shock in the region close to the shock front. This is indicated in Fig. 1, where with colored dots the post bounce times (x-axis) and radial distances from the center of the SN (y-axis) are shown for which we find crossings in D1, assuming F(γ) to be a quadratic polynomial. The top (bottom) row refers to the 20 (18.6)M model. The left (right) column refers to the case with (without) muon creation. The red, blue, green and orange dots refer to crossings for G eµ , G eτ , G µτ and G e , respectively. Note that, in order to have a clearer separation between different dots, the red, green and orange ones are slightly shifted horizontally from their correspondent post bounce time. The black line indicates the evolution of the shock in time and space. Obviously, the 20M model is not exploding, whereas the 18.6M case is exploding, where the explosion was triggered artificially around t pb ≈ 0.26 s. Fig. 2 is equivalent to Fig. 1, but it refers to the (not exploding) 15M case of D2, which is not including muon creation. Note that only three time snapshots are available in the public version of the data. To have a fair comparison with Fig. 1, we present the crossings found by looking directly at angular distributions with black dots, whereas those obtained with the moments method proposed in Ref. [40] are given by orange dots. Fig. 3 displays the shape of F(γ) that allows to capture the crossings. In all cases the functional form of F(γ) is assumed to be a quadratic polynomial. For the sake of brevity, this figure only contains information regarding the 20M model with muon creation of D1 at t = 0.3 s and r = 100 km. Similar results are obtained in all other cases. Fig. 4 displays the shape of the angular distribution uniquely provided by D2, for three radial distances where crossings are present in the 15M model at 0.15 s. For the first two distances crossings are found by directly looking at changes in the sign of G e and not by using the moments method proposed in Ref. [40].
The most important original findings of our analysis are: 1. Starting from t 0.08 s, there are crossings occurring below the SN shock in D1.
2. At fixed post bounce time, above the SN shock location there is usually a gap in terms of radial distance, where no crossing is found. The width of such a gap is time dependent and it nearly disappears when the explosion sets in the 18.6M model of D1. The gap is absent when finding crossings using directly the angular distributions.
3. There exist post bounce times and radial distances for which only crossings in the µ − τ sector are present. This is the first time that such a possibility is presented in the literature. By construction, such crossings cannot exist without considering muon creation, as usually done in the literature, since in this case G µ (γ) = G τ (γ).
4. All the captured crossings are in the backward directions. This means that the functional form of F(γ) satisfying Eq. (7) is always peaked at γ close to −1. Moreover, the crossings become more shallow at larger radial distances, since the value of F(γ) at γ = −1 increases with the distance.
5. Crossings are present for both quadratic and cubic polynomial forms of F(γ) (not shown), though they are unavoidably more numerous in the latter case. Using only quadratic polynomials for F(γ) removes crossings at lower radial distances. 6. A larger number of crossings are present in the eµ and eτ sectors in the presence of muon creation (left panels of Fig. 1). Crossings below the SN shock do not disappear in the absence of muon creation in the same model. 7. We find crossings in D2, whereas none were reported in Ref. [39]. Indeed the focus of that work was to look for deep crossings in the forward directions. Fig. 4 show that, if present, crossings come at least in pairs and that the region in γ where G e becomes negative gets larger with radial distance.

V. IINTERPRETATION OF THE RESULTS
Until now, no crossings has ever been reported in the post-shock region of one-dimensional CCSN simulations [31,39], though muon creation has never been taken into account. In particular in Ref. [39], no crossings was found, whereas in Ref. [31] they were found only in the SN pre-shock region, and associated with an enhanced coherent scattering of the more energeticν e on heavy nuclei. Our results, both with and without muon creation, confirm the findings of Ref. [31]. Indeed, in Figs. 1 and 2, most of the crossings are in the pre-shock region and they are all induced by scattering, since they are in the backward direction. This is confirmed by Fig. 3, where the shape of F (γ) has a peak at γ = −1. This is also supported by our reanalysis of D2 [39] with some further insight. As Fig. 4 shows, the crossings above the SN shock (panels for 118 km and 140 km) come in pairs: one in the very backward direction with γ −1 and the other at a larger γ, which can be even positive. The widths of these two crossings becomes larger with radial distance. Among the two, the crossing occurring at γ −1 is more difficult to explain and it has never been highlighted in the literature. Considering only the preshock region, a possible origin could be electron captures on free protons in the infalling matter. The electron neutrinos are dominantly radiated in the backward direction (or equivalently radially inward) due to the high infall velocity of the pre-shock matter. Such ν e emission might dominate theν e in the very backward direction. In principle,ν e can be similarly produced in positron captures on neutrons. Nevertheless, even if there could also be free neutrons ahead of the shock, only very few positrons are present, because the pre-shock matter is cool and degenerate. The same explanation, however, cannot be applied in the post-shock region, where the infalling velocity is much smaller.
Another interesting observation is the presence of crossings in the µ − τ sector, which can be explained as follows. Both the number densities and the average energies of ν µ are smaller than those ofν µ , but the neutralcurrent cross section for the former is larger assuming the same energy. This is due to weak magnetism corrections [61] and this is enough to compensate the not so large energy difference between ν µ andν µ (usually not larger than 0.5 MeV). This means thatν µ should be dominating over ν µ in the forward direction, whereas the opposite is happening in the backward one. In principle, the τ flavor can also play a role. But given the fact that the number densities of ν τ andν τ are very similar, the contribution from G τ (γ) should be negligible in G µτ .
One should also note that the larger the radial distance the more shallow the crossings. This is clearly visible in the third panels of each row in Fig. 3, where the values of F(γ) at γ = −1 become much larger compared to the other panels. The main reason for the increased shallowness is the geometrical reduction of the neutrino fluxes. This means that F(γ) needs to be very large in the back-  FIG. 2. Same as Fig. 1 but for the crossings in the 15M of D2. Note that the crossings can be captured here by looking directly at the ELN angular distribution (labeled as "with Ge(γ)"). In order to get insight, we have also found crossings using the moments method. Here F(µ) is assumed to be a quadratic polynomial, but its shape does not change significantly when assuming a cubic one. ward direction to be capable of capturing crossings. An intriguing aspect of our results is that some crossings are found with the moments method [40] in the post-shock region. None of these were found in Ref. [39], whereas in Ref. [31] they were found only by looking directly at the angular distributions and not at the moments. This means that probably such post-shock crossings are indeed present in all models, but they are more pronounced and easier to find in D1. This aspect is not related to the absence of muon creation in the previous works. Indeed, in D1 we find crossings below the SN shock both in the pre-and post-shock regions, regardless of the inclusion of muon creation. The reason for this difference compared to the previous studies might be related to non-equal set of theoretical inputs in the numerical models. A first possibility is the different equation of state of the PNS. Indeed, this has an influence on neutrino luminosities, neutrinosphere temperatures, mean neutrino energies and the PNS radius. For instance D1 and D2 are obtained with the following equation of states: SFHo for the former and LS220 [58] for the latter. We have checked the data from a version of D1 with LS220, with all the other inputs unchanged, and found that the pattern of crossings remains basically the same. This is shown in Fig. 5. A second possibility might be connected to PNS convection, which is included in all D1 models through a mixing-length approximation, but not considered in D2. Convection has the following consequences: it increases the PNS radius, it increases the ν e , ν µ and ν τ luminosities but reduces the electronν e luminosity relative to ν e [62,63]. It also changes the ν e andν e spectra relative to each other. Another possibility consists in different inputs in terms of microphysics. Indeed, the models in D1 and D2 use different neutrino reactions than those employed in Ref. [31].
The last aspect worth mentioning is the presence of the gap just above the SN shock, i.e. a range of radial distances for which no crossings is found with the moments method. Most likely here crossings are present, but they are harder to find since inaccessible higher order moments are required in order to capture them. Indeed, if we look directly at angular distributions in D2 we find a continuous set of crossings starting from a radial distance in the post-shock region: the smaller the distance the narrower the crossing. Another proof is that if we pass from cubic to quadratic polynomials for I α,β F most of the crossings disappearing are those at the low radial distances above the SN shock. Possible explanations of why crossings are harder to find in the gaps are the following. First, in the pre-shock region the fraction of heavy nuclei is increasing with the distance and consequently also the size of the crossings in the backward directions does. This possible connection is physically related to the fact that the temperatures in the pre-shock flow at small radii become large enough to dissociate iron-group nuclei to alpha particles, thus reducing coherent scattering cross sections. Moreover, the gap width is smaller in the 18.6M model of D1 especially after the explo-sion has set in. Under these circumstances, since the SN shock propagate outwards, the temperature in the preshock region is not high enough to produce a significant dissociation of heavy nuclei.

VI. CONCLUSIONS
A full understanding of the outcome of neutrino flavor conversions in CCSN is currently missing. Such an achievement would represent a milestone in both particle physics and astrophysics, since the next SN neutrino burst can be correctly interpreted only when the flavor composition at the source is precisely reconstructed. A step towards this ambitious goal consists in assessing whether the condition for fast conversions, i.e. a crossing in the angular distribution of the neutrino lepton number, is satisfied in realistic SN models. In this paper we have considered state-of-the-art one-dimensional CCSN simulations where muon creation is taken into account. In this context, a full three-flavor approach is mandatory, since ν µ and ν τ behave in different ways. We have concluded that muon creation has only a mild impact on the generation of crossings, though it seems to make the crossings more pronounced.
We confirm that backward crossings are relatively common above the SN shock, as first pointed out in Ref. [31]. However, we emphasize that they are also most likely common in the post-shock region, but they are not generically observed due to their dependence on different inputs used in CCSN simulations. In particular, PNS convection is taken into account in D1 [5] but not in D2 [39]. The consequences of PNS convection on the neutrino emission properties seem to enhance the depth of crossings, thus making them less difficult to be observed with the moments method.
Another original outcome of our analysis is the observation of crossings in the µτ sector, whereas none are usually expected because of the equivalent treatment of neutrinos of the muon and tau flavors. Intriguingly, there are locations in time and space where such crossings are found, where none are observed in the eµ and eτ sector. In analogy to all other crossings, these are also in the backward direction and are related to residual neutrino interactions in the free-streaming regime.
Our findings support the need for performing a similar analysis in the context of other one-dimensional models, with different sets of theoretical ingredients. Such a systematic study can augment the understanding we tried to develop in this paper. In the end, this will lead to an improved understanding of the creation of angular crossings. Moreover, our results require a conclusive confirmation from multi-dimensional models with muon creation. Two-dimensional ones have been already presented in Ref. [33] and the results from first three-dimensional simulations including muons have recently been published in Ref. [36].
Even if one assumes that the results obtained in the present work hold in the multi-dimensional models, too, the question of what is the impact of backward crossings remains unanswered. The backward crossings found in this study seem to exist only slightly below the SN shock. Nevertheless, more crossings may be present even deeper in the star, though they would be significantly more shallow. Although the occurrence of this backward ELN crossings can lead to fast instabilities with significant growth rates, it is not yet clear whether they can lead to significant flavor conversions. Indeed, these instabilities mostly affect (in the linear regime) the backward traveling neutrinos (see Fig. 7 in Ref. [64]), which are only a very tiny fraction of the neutrino population. The ultimate assessment must come through numerical simulations. One can perform simplified simulations of neutrino flavor evolution, completely decoupled from the hydrodynamics of the SN. On the other hand, a more reliable answer must come from fully including flavor conversions in self-consistent simulations of CCSNe, or treating them with an effective approach in the same context, which is computationally more affordable. We hope our work will trigger future endeavors in this direction.