Neutrino Decoupling Is Altered by Flavor Conversion

The large neutrino density in the deep interior of core-collapse supernovae and compact binary merger remnants makes neutrino flavor evolution non-linear because of the coherent forward scattering of neutrinos among themselves. Under the assumption of spherical symmetry, we model neutrino decoupling from matter in an idealized setup and present the first non-linear simulation of flavor evolution in the presence of charged current and neutral current collisions, as well as neutrino advection. Within our framework, we find that flavor transformation occurs before neutrinos fully decouple from matter, dynamically affecting the flavor distributions of all neutrino species and shifting the location of the neutrino decoupling surfaces. Our results call for further work as they may have implications on the explosion mechanism of supernovae, the nucleosynthesis of the heavy elements, as well as the observable neutrino signal, all of which is yet to be assessed.


I. INTRODUCTION
Neutrinos are fundamental particles in the physics of core-collapse supernovae [1,2].Neutrinos are copiously emitted as matter accretes onto the proto-neutron star and transport energy and lepton number.According to the delayed neutrino-driven mechanism [3][4][5], neutrinos are essential to revive the stalled shock wave, powering the explosion.Similarly, in neutrino-cooled disks stemming from the merger of two neutron stars or the collapse of rotating massive stars, neutrinos radiate most of the heat stemming from the turbulent flow, the neutrino pair-annihilation rate may contribute to power the jet of short gamma-ray bursts, and the neutrino-driven outflow launched from the disk contributes, at least partially, to the synthesis of the heavy elements [6][7][8][9][10][11].
Despite swift progress occurring in the last decade, neutrinos are treated as radiation in hydrodynamic simulations of core-collapse supernovae and compact binary merger remnants [12][13][14][15][16].However, they should be modeled through a quantum kinetic approach describing the evolution of the distributions of all six neutrino species (three neutrino flavors and their antiparticles), including flavor mixing [17][18][19][20][21].This approximation is adopted because solving the full quantum kinetic transport of neutrinos is a formidable task, entailing the solution of a 7-dimensional non-linear problem (involving time, three spatial coordinates, energy, and two angular degrees of freedom), whose characteristic quantities rapidly vary by many orders of magnitude [22][23][24].To date, this problem is not solvable with available computational resources.
Flavor evolution depends on the neutrino mixing parameters, as well as on the coherent forward scattering of neutrinos onto each other and on matter [17,[25][26][27][28][29][30].Neutrino-neutrino interaction makes the equations of motion non-linear and correlates the flavor evolution of neutrinos of different momenta.A full solution of the quantum neutrino transport problem does not exist yet.
The omission of neutrino conversion in hydrodynamical simulations was not deemed worrisome in supernovae as it was expected to occur beyond the shock radius, and thus, not having an impact on the shock revival [31].However, this picture has been shaken by recent developments [32][33][34][35][36], hinting that neutrino conversion could already develop in the proximity of the neutrino decoupling region [37][38][39][40][41][42][43][44][45][46][47], with a characteristic time scale that is proportional to the number density of neutrinos and antineutrinos (hence dubbed "fast" to distinguish it from the "slower" conversion driven by the vacuum frequency).The impact of flavor conversion on the supernova physics remains to be understood.Additionally, in compact binary mergers, flavor conversion could affect the disk cooling rate as well as the neutron-to-proton ratio [48][49][50][51][52].
In the context of fast flavor conversion, an important role is played by the angular distributions of (anti)neutrinos: a crossing between the angular distributions of ν e and νe , or a change of sign in the electron neutrino lepton number (ELN) within a certain angular range, constitutes a necessary condition to trigger a flavor instability [36,53,54], in the limit of vanishing masssquared difference [24,32,34,36,38,[55][56][57][58][59][60].The amount of induced flavor conversion is, however, unknown [61].In the light of these developments, it becomes of paramount importance to ascertain whether the presence of ELN crossings is modified by flavor conversion and vice-versa.
In the supernova core, for example, neutrinos are trapped with a mean-free path of O(10) m at typical baryon densities of O(10 14 ) g/cm 3 [62].As the distance from the supernova increases and the baryon density decreases, neutrinos start decoupling from matter.However, not all flavors of neutrinos decouple simultaneously due to the flavor dependent cross section of neutrinos.Their angular distributions are initially almost isotropic and, as the density decreases, become forward peaked and ELN crossings can form, also because of collisions and advection [46,63].Recent work shows that, if flavor conversion should occur in the decoupling region, a non-trivial feedback between neutrino flavor transformations and collisions could take place [63][64][65][66][67][68][69][70].Similarly, neutrino advection and the spatial dimensionality of the problem could also dynamically modify the conditions for flavor conversion [56,71].
In this paper, we present the first investigation of flavor conversion as the angular distributions of neutrinos evolve in the collisional regime within an idealized toymodel shell and in the presence of neutrino advection.Our work is organized as follows.We introduce the mean field equations of motion in Sec.II, together with the modeling of the collision term.The decoupling of neutrinos in the absence and then in the presence of flavor conversion is discussed in Secs.III and IV, respectively.Section V focuses on the interplay between slow and fast conversion.The characteristic time scales of the problem are instead discussed in Sec.VI.Finally, our findings are summarized in Sec.VII.Appendix A provides additional material comparing simulations with our standard resolution with others obtained with a smaller number of bins in radius.

II. MEAN FIELD EQUATIONS OF NEUTRINOS
In this section, we introduce the kinetic equations of neutrinos.We also outline our modeling of the collision term, which mimics the hierarchy of decoupling of the different neutrino flavors, despite being heuristic.

A. Kinetic equations
For the sake of simplicity, we work in the two-flavor approximation, i.e. in the (ν e , ν x ) basis [57,59,72,73].The [anti]neutrino field is described through a 2 × 2 density matrix, ρ( r, p, t) [ρ( r, p, t)].The diagonal elements of the density matrix (ρ ii , with i = e, x) represent the occupation numbers of neutrinos of different species, while the off-diagonal terms (ρ ij ) encode flavor coherence.The evolution of the neutrino field at time t, location r and momentum p is governed by the following equation [74]: We work in spherical symmetry and assume only one energy mode E; hence, ρ( r, p, t) = ρ(r, E, θ, t), with θ being the local polar angle defined with respect to the radial direction and at a given r.On the left-hand side of Eq. 1, the advective term is, in agreement with e.g.Refs.[75,76].The commutator on the right hand side of Eq. 1 encodes the flavor evolution of neutrinos.The Hamiltonian, H = H vac + H matt + H νν , governs the coherent evolution of neutrinos and is composed of the vacuum and neutrino self-interaction term: where we have absorbed the constant factor coming from the integration over the azimuthal angle in µ 0 and suppressed the dependence on r and t for the sake of brevity.
The term, H matt , linked to matter effects should also appear in H.However, the effect of the matter potential is to suppress the vacuum mixing angle ϑ V [77], hence we set H matt = 0 and use ϑ V = 10 −3 instead 1 .The vacuum frequency is ω = ∆m 2 /2E, with ∆m 2 being the squared mass difference.For antineutrinos, ω → −ω in H vac , while all other terms of H remain unchanged.

B. The collision term
The collision term in Eq. 1, C ≡ C( r, E, t) = C emission + C absorb + C dir-ch , includes emission, absorption, and direction-changing collisions, respectively [78].It encapsulates all the reactions of neutrinos with the matter background [13,63,79,80].We assume that the collision term only affects the diagonal components of the density matrices 2 , as the time scales associated with flavor evolution are shorter than the collision ones [24].
For the sake of simplicity, given that we intend to explore the interplay between the neutrino decoupling and the flavor conversion for the very first time, we rely on a heuristic collision term, that allows for neutrinos to transition from the trapping regime to free streaming within a small spatial range (i.e., [15,30] km in our simulations, while providing decoupling regions and decoupling hierarchy among the different flavors inspired by Ref. [37]), since we are limited by the size of our simulation shelldue to technical complications induced by the treatment of flavor conversion.For simplicity, we retain the dependence of C on r and neglect its angular dependence because it is dominated by the interaction of neutrinos with nucleons that are heavier than the neutrino average 1 A smaller ϑ V would lead to a longer simulation time needed to achieve the quasi-steady state configuration.We have checked that the quasi-steady state configuration is virtually unaffected [for ϑ V = 10 −8 , the quasi-steady state is reached at 10 −4 s (as for the run with ϑ V = 10 −3 ), and we find a relative error of 1.7% averaged over the radius with respect to the case with ϑ V = 10 −3 ]. 2 We have tested that the inclusion of the off-diagonal components of the collision terms following the method adopted in Ref. [76] results in a negligible error on the final flavor configuration (the relative error averaged over the radius is 2.2%).energies: C νe,νe,νx,νx C νe,νe,νx,νx ρ ii (cos θ) Each of the above equations refers to all flavors as denoted by the superscripts, and the flavor-dependent length-scales, λ νi , are tabulated in Table I and shown in Fig. 1.Equations 5 and 6 determine the rate at which neutrinos traveling in the direction θ are created and absorbed, respectively.Equation 7 represents the probability that a neutrino traveling along θ will change to the θ direction.While Eqs. 5 and 6 do not conserve the number of neutrinos, Eq. 7 preserves the number of neutrinos.
In our simulations, neutrinos of all flavors are generated through collisions.The ratio between C emission and C absorb determines ρ ii for all flavors at r min = 15 km.Since only C emission and C absorb play a role in the trapping regime, the number density at r min can be determined analytically in the steady state configuration.At larger r, C dir-ch also comes into play in shaping the angular distributions and no analytical solution exists.In addition, in our simplified setup, we neglect Pauli blocking [81,82] and assume that the neutrino chemical potentials are negligible.
To investigate the decoupling of neutrinos, we focus on the simulation shell from 15 to 30 km.We use 7500 × 150 bins in r and cos θ, respectively.Additional technical details on the numerical implementation and convergence are provided in Secs.III and VI, as well as Appendix A and Ref. [83].
The radial profiles for the heuristic collision terms are shown in Fig. 1 (see also Table I).We can see that λ −1 dir−ch is the largest contribution for all flavors at small radii.The direction changing scattering is dominated by neutral current interactions with nucleons.Because of our single energy approximation, we adopt an heuristic approach to model the collision term such that ν e 's have the smallest mean free path, while non-electron type neutrinos have the longest mean free path, as expected [13,79,80].Moreover, the fact that the νe mean Inverse length-scales associated with emission, absorption, and direction-changing scattering for νe (in red), νe (in green), and νx (in blue) as functions of radius.The functional forms of these terms are reported in Table I, see also Eqs. 5-7.
free path is larger than the ν e one ensures that the νe 's decouple at a smaller radius than ν e 's [84,85].As the radius increases, all inverse length scales decrease, until they are negligible for r 25 km, when neutrinos decouple from matter as shown in the top panel of Fig. 3.
As neutrinos start to decouple from matter, the angular distribution for each neutrino flavor becomes forward peaked and the number density which is proportional to the diagonal components of the density matrix falls off.This is clearly visible from the radial profile of the number density for each flavor [which corresponds to 2, in the steady state configuration and in the absence of neutrino mixing.At radii larger than ∼ 25 km, the emission and absorbtion terms are almost zero.In this limit we find that the neutrino number density falls off like 1/r 2 as expected due to conservation of neutrino number for each flavor.Also, since there is almost no absorption or emission of neutrinos in the region of r 25 km, the angular distributions become forward peaked as the only neutrinos that are present at these radii are the ones which have escaped the core and travel along directions that are close to the radial one.

III. NEUTRINO DECOUPLING IN THE ABSENCE OF FLAVOR CONVERSION
For all the simulations presented in this paper, we use the central difference method for spatial derivatives, while an adaptive multi-step Adams-Bashforth-Moulton method is adopted for the temporal evolution with absolute and relative tolerances of 10 −12 .The central difference method adopted to calculate the spatial derivative and the derivatives with respect to cos θ cannot be applied to the edges of the simulation shell and hence to the related derivatives.As for the derivative with respect to cos θ, we use the forward and backward difference method at cos θ = −1 and cos θ = 1, respectively.For r = r min , the boundary condition is set by requiring a classical steady state solution.For r = r max , there are two distinct domains that need to be considered independently.For cos θ > 0, neutrinos travel outwards and, therefore, the boundary condition does not need to be specified.For cos θ < 0, the neutrino flux for all flavors is set to zero.
The kind of ELN crossing could affect the development of flavor conversion physics [24,42].We choose a representative case in this work.However, we investigate the flavor conversion physics for a set of different ELN crossings (i.e., different collision terms) in Ref. [83] and show that our main conclusions are not qualitatively affected.
The flavor configuration displayed in Fig. 3 represents the classical steady state configuration obtained in the absence of neutrino flavor transformation.Within our setup, the neutrino angular distributions are populated through collisions and advection across the simulation shell.While we rely on heuristic collision terms and focus on an idealized simulation shell, the decoupling regions and the hierarchy among the different flavors qualitatively reproduce the ones obtained from hydrodynamical simulations, e.g.Ref. [37], while also developing ELN crossings.
The decoupling of neutrinos from matter is driven by the flavor-dependent collision term along with neutrino advection [63,86].Solving Eq. 1 for H = 0 (i.e. in the absence of flavor conversion) and no neutrinos in the simulation shell at t = 0, the collision terms in Eqs.5-7 lead to a steady state configuration for each neutrino species after a sufficiently large time interval (t 10 −4 s).The top panel of Fig. 3 shows a contour plot of ρ ee in the plane spanned by cos θ and r.At small radii, the angular distribution (15 r 18 km) is uniform in cos θ; as r increases, the angular distribution becomes gradually forward peaked and its backward part (cos θ 0) is progressively emptied as full decoupling is reached.A similar behavior holds for all (anti)neutrino species.
As the presence of ELN crossings is crucial to the development of fast flavor conversion, the middle panel of Fig. 3 displays a contour of ρ ee − ρee in the plane spanned by cos θ and r.ELN crossings are not prominent at small r, but they develop at r 18 km and affect different angular regions as r increases.The dark blue band to the left of the dashed line in Fig. 3 is the region with an excess of ν e over νe due to the fact that the νe decoupling radius is smaller than the ν e one.The bottom panel of Fig. 3 displays snapshots of the angular distributions of ν e , νe , and ν x = νx at two representative radii.As C evolves as a function of r, the angular distribution of each flavor becomes forward peaked and ELN crossings are visible.
The dashed lines in the the top panel of Fig. 4 show the radial evolution of the angle-averaged occupation number for ν e and ν x .In order to determine the decoupling surface, we compute the flux factor: and approximately define the radius of decoupling as the radius at which the flux factor is F νi = 1/3 [37,49].The corresponding decoupling surfaces of ν e and ν x are shown in the top panel of Fig. 4.
Figure 4 shows fascinating insights into the effect of flavor conversion on neutrino decoupling.The solid lines in the top panel of Fig. 4 show the angle averaged occupation numbers of ν e and ν x .With respect to the case without flavor conversion (dashed lines in Fig. 4), we see that flavor conversion pushes the distributions of ν e and ν x towards each other, sensibly modifying them with respect to the case without flavor conversion.However, we stress that flavor equipartition is not a general outcome, but it is linked to the specific flavor setup adopted in this paper; we have found other flavor configurations that do not lead to flavor equipartition (results not shown here), see also Refs.[83,87].The impact of neutrino flavor evolution is evident in the shift of the radius at which the angle-averaged ρ ee and ρ xx start falling, with consequences on the decoupling radius.One of the consequences of neutrino flavor conversion is a shift in the radius of neutrino decoupling.The top panel of Fig. 4 shows the flavor dependent decoupling surfaces with (solid) and without (dashed) flavor conversion for comparison.It is also important to note that neutrinos start to change their flavor at radii much smaller than their decoupling surfaces; this is in stark contrast with the approximation of flavor-independent neutrinosphere commonly adopted in the literature [22,23,[26][27][28][29], under the naive expectation that flavor evolution near the neutrinosphere would not be present due to collisional damping [88][89][90].
The bottom panel of Fig. 4 represents the difference in the ν e occupation number in the absence and presence of neutrino mixing.As r increases, a region with a deficit of ν e forms due to flavor conversion (red band).This physics cannot be captured by any calculation that does not include advection and collision, together with flavor transformation-see Fig. 5 and the animations for comparison.
After obtaining a steady state configuration for all flavors (see Fig. 4), we switch on flavor conversion and explore its interplay with neutrino advection and collisions.In order to highlight the dynamical effects of neutrino advection and collisions on flavor conversion, we show in Fig. 5 the correspondent flavor outcome obtained by switching off collisions and advection once the steady state flavor configuration is achieved, i.e. by letting flavor transformation operate in the absence of collisions and advection.
From Fig. 5, one can see that collisions and neutrino advection not only affect the final angle-integrated flavor configuration, but also allow to spread the flavor instabilities across angular regions and spatial ranges [56,64,65].It is also worth noticing that the absence of neutrino advection makes the evolution of neutrino flavor much more oscillatory in nature.

V. INTERPLAY BETWEEN SLOW AND FAST FLAVOR INSTABILITIES
Interestingly, despite the presence of ELN crossings and because of the shape of our ELN angular distributions (see middle panel of Fig. 3) only a small radial range centered on r 21 km is affected by flavor instabilities in the limit of ω = 0 (i.e., when only fast conversion occurs).Because of ω = 0 and the system setup, slow collective neutrino transformations are also triggered at high densities, in contrast to what is commonly assumed in the literature.
In this section, we choose to carry out the linear stability analysis for the homogeneous mode only, as an illustrative example.However, the numerical simulations automatically take into account all modes.Hence, the interplay between slow and fast conversion is accounted for all modes consistently in the numerical simulations.
The investigation of the interplay between fast and slow instabilities for non-homogeneous modes would require the development of a novel analytical framework because of the importance of the scattering terms in our system setup; this task is left to future dedicated work.
The linear stability analysis consists of linearizing the equations of motion and calculating the growth rate of the flavor instability [91].By setting the collision and the advective terms to zero ( v = 0), expanding Eq. 1 in the off-diagonal component, and ignoring O(ρ 2 ex ) and higher order terms, the solutions are of the form: The collective nature of flavor evolution ensures that the eigenvalue Ω is the same for neutrinos and antineutrinos and also the same for all values of cos θ.When a flavor instability is present, Ω has a positive imaginary component, called growth rate and denoted by κ.It should be noted that, although the complex eigenvalues are always present in complex conjugate pairs in the absence of the collision term, this is not the case in the presence of collisions [92], which we do not consider here.However, the presence of a neutrino flavor instability does not necessarily imply significant flavor transformation [93].
Figure 6 shows the growth rate as a function of radius for ω = 0 (fast conversion) and ω = 0.32 km −1 adopted in our work.Interestingly, for most of the radial range, the neutrino flavor instability does not exist for ω = 0 for the homogeneous mode because of the shape of our ELN distributions.This implies that initially the neutrino flavor evolution is not driven by ELN crossings, but by the vacuum term for k = 0. Fast flavor instabilities are only present in a small radial range around ∼ 21 km.In addition, once triggered, fast flavor mixing is further affected by the vacuum frequency [55].We stress that these findings imply an interplay between fast and slow conversion and do not intend to suggest that the slow modes are dominant for any k.By comparing Fig. 6 with Figs. 4  and 5, we note that although the neutrino flavor instability exists for r 18 km, the magnitude of the neutrino flavor transformation is not directly correlated to κ, in agreement with the findings of Ref. [93] in the context of fast flavor transformation.

VI. CHARACTERISTIC SCALES OF THE PROBLEM
In the context of fast conversion, the time scale of flavor evolution can be as large as µ −1 in special cases (e.g., for ω = 0 and in the absence of collisions and advection).For more realistic cases the time scales are smaller by a few orders of magnitude.
One might naively confuse the time scale with the length scale over which spatial structure is present and in turn assume that this would infer the spatial resolution required to perform a reliable numerical study.However, as demonstrated in Ref.
[93], we stress that spatial and temporal scales are not identical when advection is taken into account; it is possible to get reliable results by coarse-graining over the small spatial structure, if the latter is at all present, and reproduce the average flavor ratio, see also Appendix A and Ref. [93].The typical length scale of our problem should not be considered to be O(µ −1 0 ), as in e.g.Refs.[71,87,[94][95][96]]-these works rely on a different system setup, with structures at length scales of O(µ −1 0 ) and periodic boundary conditions that may lead to cascades to smaller length scale with time [97].
In our case, the characteristic length scale of the problem is determined by the flavor conversion one, as well as the collision and advection scales.If neutrinos were assumed to be emitted by a single surface, as sketched in the left panel of Fig. 7, the typical length scale of the problem would be given by the product of the flavor conversion time scale and the speed of light (i.e., the advection velocity).This would be the length scale seen in a setup similar to the neutrino-bulb model, where the decoupling region has an infinitesimal extent.However, our decoupling region is much more extended than the oscillation length scale, see right panel of Fig. 7. Hence, at a given point, the resultant neutrino field is given by the neutrinos emitted from various locations in the decoupling region (i.e. the small-scale structures that one would see when considering one single emitting surface are smeared as a result of the superposition of oscillating patterns overlapping with each other, as sketched in the right panel of Fig. 7).
For a density profile that is rapidly decreasing, as ours, the extent of the decoupling region is roughly determined by the inverse of the collision term, ∼ C −1 .Although it is not possible to estimate a priori the exact length scale associated with the neutrino field in the presence of flavor transformation, it is clear that the resultant length scale should be between the oscillation and the collision length scales.In summary, µ −1 0 constitutes the minimum length scale expected for flavor conversion, the maximum being the collision length scale.As a consequence spatial resolution of O(µ −1 0 ) may not be needed for numerical convergence because of the interplay of flavor conversion with collisions and advection.).For a density profile that is rapidly decreasing, the decoupling region has a width that roughly scales as ∼ C −1 (right panel).In the first case, the typical length scale of the neutrino oscillation pattern is given by the product of the flavor conversion time scale ∆t and the advection velocity.When an extended decoupling region is considered, the resultant oscillation pattern (black curve) is given by the superposition of multiple oscillation patterns (gray curves), each separated from the other ones by O(C −1 ).

VII. OUTLOOK
For the first time, we follow the flavor evolution while neutrinos decouple in an idealized shell, in the presence of collisions and advection.We find that neutrino flavor conversion occurs well before all flavors decouple from the matter background and the neutrino decoupling surfaces are affected by flavor transformation.
While dealing with flavor transformation, advection and collisions simultaneously, the collision term adopted in this work relies on heuristic functions that allow us to transition from uniform to forward peaked neutrino distributions within a small spatial region because of the technical challenges induced by the treatment of flavor conversion in the stellar core.A further exacerbation of our findings could be determined by the inclusion of all six flavors, the energy dependence of the neutrino distributions, and the azimuthal emission angle of neutrinos [55,57,[71][72][73]98].Additional effects on the decoupling physics may arise from relaxing the assumption of spatial homogeneity [99].Not only can spatial inhomogeneity lead to additional flavor instabilities, neutrino advection causes dynamical effects which can only manifest in an inhomogeneous system [56].
Despite its caveats and the idealized simulation setup, this work clearly demonstrates that the neutrino decoupling physics is affected by flavor transformation, with possible implications on the neutrino energy deposition in the supernova gain layer and the physics of compact merger remnants that remain to be assessed.The fact that the dynamics of neutrino decoupling changes due to neutrino flavor transformation implies that the neutrino properties at the radius of decoupling may be different from the ones usually computed in hydrodynamic simula-tions treating neutrinos as radiation.However, a robust assessment of the radiated neutrino spectra in the presence of flavor conversion will have to take into account the hydrodynamic feedback on the thermodynamic properties and a more realistic implementation of the collision term.The occurrence of flavor conversion in the vicinity of the flavor-dependent neutrinospheres puts in perspective the existing literature naively assuming flavorindependent neutrinospheres and a larger radius for the onset of flavor conversion.and 10 2 km −1 with 150 spatial bins and compare the results with simulations performed with 7500 spatial bins adopted in the main text (i.e., 50 times higher spatial resolution).The choice of a smaller µ 0 allows us to perform numerical simulations with a spatial resolution that is accurate enough to resolve length scales of the order of µ −1 0 or smaller.Figure 8 shows results for µ 0 = 10 4 km −1 for our two simulations with 150 and 7500 radial bins (same as the one presented in the main text), with all other simulation inputs unchanged with respect to our benchmark model (including the number of angular bins kept fixed to 150).The top panel shows the difference in the electron neutrino number density, with and and without flavor conversion obtained using 150 spatial bins.Despite the appearance of small scale structure due to the coarser graining in radius,the results are same as the ones presented in Fig. 4 within reasonable numerical errors.This is evident from the middle panel of Fig. 8, where we show the comparison between the number densities averaged over angle for the simulations with 150 and 7500 spatial bins.The bottom panel of Fig. 8 displays the relative error between the two simulations.The error averaged over the radial range is 1.1%.This shows that 150 spatial bins are sufficient to capture the main features of the flavor evolution, despite the appearance of small scale structures.
Figure 9 shows a comparison between simulations using smaller values of µ 0 (100 and 1000 km −1 , left and right panels, respectively).From the bottom panels, we can see that the average relative error is 1.66 and 1.63% for µ 0 = 100 and 1000 km −1 , respectively.In all cases, the average relative error between the simulations with 150 and 7500 radial bins can reach up to 4-5% in a few radial bins but on average is less than 1-2%.The relative error has been computed by coarsegraining the results with 7500 radial bins and averaging over batches of 50 radial bins, to compare the results with the ones from the numerical simulation that used 150 radial bins.Since a quasi-steady state configuration is achieved, the relative errors quoted above overestimate the actual error due small time-dependent fluctuations that are present at every given point when the simulation is stopped.
It is easy to see that small time-dependent fluctuations are the main source of the error by investigating the error without averaging over the angle.Figure 10 shows the absolute error in ρ mix ee − ρ no−mix ee in the left panel and the relative error in the right panel for the simulations with µ 0 = 10 4 km −1 (see also Figs. 4 and 8).
The absolute error in the left panel of Fig. 10 provides an accurate representation of the error in our simulation, which is about 0.1-0.15throughout the simulation domain.The relative error is shown in the right panel of Fig. 10 for completeness.We plot the relative error only for the regions where ρ no−mix ee 0.1 (see red contour in the right panel of Fig. 10).To guide the eye, we also show an orange contour where ρ no−mix ee 0.5.The right panel shows that the relative error between the low-and high-resolution simulations is less than 5-10%.However, it can reach up to 40% for cos θ < 0 where the relative error is not an informative quantity since there are very few neutrinos in the tail of the angular distribution.We stress that the relative error is ill-defined in such regions and should be considered with extreme caution.This becomes clear by looking at Fig. 11 where we show the ν e angular distribution when the quasi-steady state is reached for r = 21 and r = 22 km with µ 0 = 10 4 km −1 .We plot ρ mix ee obtained from the simulation with 150 spatial bins, alongside with the correspondent angular distributions obtained from the simulation with 7500 spatial bins.The angular distributions obtained for the simulations with 150 and 7500 bins are in very good agreement, yet if one were to compute the relative error e.g. for cos θ = −0.9,one would obtain an error above 40%; this clearly shows the misleading information provided by the relative error in the right panel of Fig. 10 for some regions of the parameter space spanned by cos θ and r.
Figure 12 shows the radial evolution of ρ mix ee obtained from the simulations with 150 and 7500 spatial bins for cos θ = −0.5, 0, and 0.5 (from top to bottom panels, respectively).One can clearly see that the low-resolution curve tracks the average behavior of the high-resolution one.
It is worth highlighting that the value of µ 0 reported above is the self-interaction strength at r min ; as neutrinos start decoupling, the effective strength of self-interaction decreases as the fourth power of the radius, which is automatically taken into account in our simulations.12. Radial profile of ρee after flavor conversion, extracted at cos θ = −0.5, 0, and 0.5 (from top to bottom, respectively) for µ0 = 10 4 km −1 for the simulations with 150 (red) and 7500 (green) spatial bins.

FIG. 2 .
FIG.2.Neutrino number densities of νe (in red), νe (in green) and νx (in blue) as functions of radius in the steady state configuration in the absence of neutrino mixing.The number density has no radial dependence in the trapping regime (r 18 km), while it falls of like 1/r 2 at larger radii (r 25 km), as expected because the absorption and emission terms are very small.

FIG. 3 .
FIG. 3. Steady state neutrino flavor configuration in the absence of flavor mixing.These distributions (generated by imposing H = 0 in Eq. 1 and extracted at t = 10 −4 s) are the ones adopted as input to investigate the effects of flavor conversion.Top panel: Contour plot of ρee (proportional to the νe number density) in the plane spanned by cos θ and r.Middle panel: Same as in the top panel but for ρee − ρee (proportional to the ELN density).The dashed line marks the region where ELN crossings develop.Bottom panel: Angular distributions of νe, νe and νx at r = 19 km (solid) and 29 km (dotted).As r increases, the angular distributions become prominently forward peaked in a flavor-dependent fashion and ELN crossings develop.

FIG. 4 .
FIG.4.Steady state neutrino flavor configuration in the presence of flavor mixing (extracted 5 × 10 −5 s after the classical steady state configuration is reached in our simulation).Top: Angle averaged neutrino occupation numbers of νe (in red) and νx (in blue) with (solid lines) and without (dashed lines) neutrino mixing as functions of the radius.The vertical lines mark the radii of decoupling (approximately defined as the radius at which Fν i = 1/3).A similar trend holds for antineutrinos; however, note that flavor conversion induces a difference between νx and νx.Bottom: Contour plot of the difference between the νe occupation number without (when the classical steady state configuration is achieved) and with neutrino mixing in the plane spanned by cos θ and r.Due to the collective nature of the neutrino flavor evolution and flavor lepton number conservation, the corresponding heatmap for antineutrinos looks very similar and is not shown.

FIG. 5 .
FIG.5.Final flavor configuration obtained by switching on flavor mixing, after the classical steady state is reached, and switching off collisions and neutrino advection (i.e., H = 0, C = 0, and v = 0 in Eq. 1).The top and bottom panels are the analogous of the ones in Fig.2, which instead take into account the dynamical effects on flavor mixing due to collisions and advection (i.e., H = 0, C = 0, and v = 0 in Eq. 1).

FIG. 6 .
FIG.6.Growth rate of the flavor instability as a function of radius obtained by adopting the steady state flavor distributions.The orange line shows the growth rate for ω = 0 (fast conversion), whereas the blue line shows the growth for ω = 0.32 km −1 (slow conversion).Contrary to expectations, flavor mixing is induced by slow instabilities; fast flavor instabilities are only present around ∼ 21 km.

∆t∼ C − 1 FIG. 7 .
FIG.7.Schematic representation of neutrinos being emitted by a decoupling region with infinitesimal extent (left panel) and neutrinos emitted from an extended decoupling region composed of multiple emitting surfaces, each separated from the former and the following one by O(C −1 ).For a density profile that is rapidly decreasing, the decoupling region has a width that roughly scales as ∼ C −1 (right panel).In the first case, the typical length scale of the neutrino oscillation pattern is given by the product of the flavor conversion time scale ∆t and the advection velocity.When an extended decoupling region is considered, the resultant oscillation pattern (black curve) is given by the superposition of multiple oscillation patterns (gray curves), each separated from the other ones by O(C −1 ).

= 10 4 km − 1 FIG. 8 .
FIG. 8. Comparison between simulations with a different number of radial bins for µ0 = 10 4 km −1 .The top panel shows the difference between ρee with and without mixing for µ0 = 10 4 km −1 computed using 150 spatial bins.The middle panel shows the initial angle averaged ρee as red dashed line.The red solid line and green solid line shows are the angle averaged ρee calculated using 150 and 7500 spatial bins.The red solid line is difficult to see due to the overlap with the green solid line.The bottom panel shows the relative error between the two simulations as a function of r.The overall error averaged over radius is 1.1%.

FIG. 11 .
FIG. 11.Angular distributions of νe after flavor conversion, extracted at r = 21 and r = 22 km and for µ0 = 10 4 km −1 for the simulations with 150 (red) and 7500 (green) radial bins.The two curves are in very good agreement.