Turbulence patterns and neutrino flavor transitions in highresolution supernova models
Abstract
During the shockwave propagation in a corecollapse supernova (SN), matter turbulence may affect neutrino flavor conversion probabilities. Such effects have been usually studied by adding parametrized smallscale random fluctuations (with arbitrary amplitude) on top of coarse, spherically symmetric matter density profiles. Recently, however, twodimensional (2D) SN models have reached a space resolution high enough to directly trace anisotropic density profiles, down to scales smaller than the typical neutrino oscillation length. In this context, we analyze the statistical properties of a large set of SN matter density profiles obtained in a highresolution 2D simulation, focusing on a postbounce time (2 s) suited to study shockwave effects on neutrino propagation on scales as small as O(100) km and possibly below. We clearly find the imprint of a broken (KolmogorovKraichnan) powerlaw structure, as generically expected in 2D turbulence spectra. We then compute the flavor evolution of SN neutrinos along representative realizations of the turbulent matter density profiles, and observe no or modest damping of the neutrino crossing probabilities on their way through the shock wave. In order to check the effect of possibly unresolved fluctuations at scales below O(100) km, we also apply a randomization procedure anchored to the power spectrum calculated from the simulation, and find consistent results within fluctuations. These results show the importance of anchoring turbulence effects on SN neutrinos to realistic, finegrained SN models.
a]Enrico Borriello, b]Sovan Chakraborty, c]HansThomas Janka, d]Eligio Lisi, a]Alessandro Mirizzi
Prepared for submission to JCAP
Turbulence patterns and neutrino flavor transitions in highresolution supernova models

II. Institut für Theoretische Physik, Universität Hamburg
Luruper Chaussee 149, D22761 Hamburg, Germany 
MaxPlanckInstitut für Physik (WernerHeisenbergInstitut)
Föhringer Ring 6, D80805 München, Germany 
MaxPlanckInstitut für Astrophysik, KarlSchwarzschildStr. 1, 85748 Garching, Germany

Istituto Nazionale di Fisica Nucleare, Sezione di Bari, Via Orabona 4, 70126 Bari, Italy
Contents
1 Introduction
The neutrino flux emitted during a corecollapse supernova (SN) explosion represents an intriguing astrophysical case where flavor transformations exhibit a strong sensitivity both on the corecollapse dynamics and on unknown neutrino properties (such as the mass hierarchy). The relevant phenomena depend strongly on the background fermion density and thus on the radial distance from the SN center.
In the deepest SN regions [ km] the neutrino density is so high that neutrinoneutrino interactions become dominant [1, 2] and can induce a neutrino refractive index responsible for collective flavor conversions; see, e.g., [3, 4, 5, 6] and the recent review [7]. The main observable features of such flavor transitions consist in an exchange of the spectrum of the electron species () with the nonelectron ones () in certain energy intervals, giving rise to socalled spectral “swaps” and “splits” [8, 9, 10].
While neutrinos propagate towards larger radii, their collective conversions get suppressed as their distribution becomes increasingly radially beamed so that eventually the flavor conversions are dominated by the ordinary matter. In this context, many studies have analyzed possible signatures of the MikheyevSmirnovWolfenstein (MSW) matter effect [11, 12] in increasingly realistic and detailed SN density profiles. For example, since the neutrino conversion probabilities depend on the timedependent matter profile, a highstatistics SN neutrino observation may also reveal signatures of shockwave propagation at late times (i.e. s after the corebounce) [13, 14, 15, 16, 17]. The transient violation of the adiabaticity condition when neutrinos cross the shockfront would then emerge as an observable modulation of the neutrino signal. This signature would be particularly interesting to follow in “realtime” the shockwave propagation, as well as to probe the neutrino mass hierarchy.
A realistic characterization of matter effects during neutrino propagation across the SN shock wave must also take into account stochastic density fluctuations and inhomogeneities of various magnitudes and correlation lengths in the ejecta layer in the wake of the shock front. These anisotropies and perturbations are a result of hydrodynamic instabilities between the protoneutron star and the SN shock during the very early stages of the SN explosion. They lead to largescale explosion asymmetries and turbulence in a dense shell of shockaccelerated ejecta, which subsequently also seed secondary instabilities in the outer shells of the exploding star (e.g., Refs. [18, 19, 20, 21, 22, 23, 24]).
Neutrino flavor conversions in a stochastic matter background have been subject of intense investigations both in a general context [25, 26, 27, 28, 29, 30, 31] and specifically in relation to SN neutrinos [32, 33, 34, 36, 37, 38, 39, 40]. In general, one expects that stochastic matter fluctuations of sufficiently large amplitude may suppress flavor conversions and eventually lead to flavor equilibration between electron and nonelectron species. Therefore, the spectral properties of the fluctuations are very important for understanding the neutrino signal emerging from a corecollapse SN.
Ideally, the study of matter density effects on neutrinos propagating through a SN shock wave should be based on input matter profiles derived from advanced and preferably nonisotropic SN models, with space resolution smaller than typical oscillation lengths. Such models should be evolved to the late times (a few seconds) relevant for shockwave effects on neutrinos. In the absence of such ideal and rather demanding requirements, many phenomenological studies have characterized the turbulence in terms of “plausible” disturbances, with arbitrary or tunable amplitudes, superimposed to shockwave profiles derived by coarse, spherically symmetric SN models.
In this context, several relevant works aimed at characterizing and modeling the fluctuation spectra beyond earlier approximations [33, 34] of deltacorrelated noise. ^{1}^{1}1Neutrino propagation through non deltacorrelated turbulence was also studied in [35]. In particular, in [36] a general description of the turbulence was proposed, by modeling the density fluctuations via a Kolmogorovtype power spectrum. This model would allow to include matter fluctuations on all the scales between the shockfront radius and the cutoff scale by viscous dissipation. A “shadow effect” on the shockwave signature was estimated for sufficiently large fluctuations. In [37] it was described how to generate a fluctuation powerspectrum on top of an undisturbed density profile through the “randomization method” [38]. Parametric studies were recently performed to investigate the dependence of the damping effects on the amplitude of the fluctuations [39], and the impact on signatures related to the shock wave and to collective oscillations [40].
We build upon previous works on this topic, by anchoring the turbulence spectra to highresolution SN simulations which have become available in recent years, although still in twodimensional (2D) approximation. In particular, 2D simulations extending from core bounce to the late postbounce times relevant for neutrino propagation have become available by works of the Garching group [18, 19]. Very interestingly, the power spectra of radial density profiles from such simulations show the typical imprint expected from 2D turbulence.^{2}^{2}2A preliminary characterization of turbulence spectra for these 2D simulations was presented in [41]. The aim of the current study is thus to perform a selfconsistent study of the SN matter turbulence by using such 2D simulation inputs, eventually improved with randomized fluctuations at small scales below O(100) km, which are not necessarily captured by the simulations themselves at large radii.
We also mention that, although first threedimensional (3D) simulations covering the relevant SN stages some seconds after core bounce have been carried out more recently [21, 23, 24], the resolution of these models is still much coarser than those of the 2D calculations used for the analysis in our paper. Moreover, due to the lack of sufficient resolution and thus the overestimation of numerical dissipation, the expected Kolmogorov power law [42, 43] cannot be recovered by preliminary analyses of the corresponding fluctuation spectra. For these reasons we shall limit ourselves to the 2D simulations of [18, 19] in this work.
Our results are presented as follows. In Sec. 2 we describe the input matter profiles as taken from our reference 2D SN model. In Sec. 3 we perform a statistical study of the matter power spectra of the density fluctuations in the different radial directions. We find that these spectra can be described by broken powerlaws, with a “typical” structure which is actually the imprint of twodimensional turbulence. In Sec. 4 we use some representative turbulent matter profiles to calculate the neutrino crossing probabilities. We find that, by taking our reference SN profiles at face value, turbulence provides no or modest damping of flavor transitions. As a relevant consequence, the shockwave signature on the neutrino crossing probabilities should hopefully remain observable in the neutrino signal from a future galactic SN explosion. In Sec. 5 we investigate the effect of fluctuations on scales smaller than the typical resolution of the simulations, by generating them with a randomization method anchored to larger scales. We find that such small scale fluctuations may provide limited additional damping in some cases, but do not qualitatively alter the main picture emerging from the 2D simulations. Finally, in Sec. 6 we discuss and summarize our findings.
2 Our numerical SN model in 2D
Our investigations are based on results of the 2D SN explosion simulations of a nonrotating, 15 star (a blue supergiant progenitor model of Supernova 1987A from [44]) described in detail in Ref. [19]. The explosions via the neutrinoheating mechanism were launched by suitably chosen neutrino luminosities from the highdensity core of the nascent neutron star, and the evolution was continuously modelled in 2D from a few milliseconds after core bounce until 20,000 seconds later. The adaptive mesh refinement code AMRA was used, which allowed to temporarily achieve a radial resolution equivalent to covering the entire star (out to its surface radius of roughly cm) with 2.6 million equidistant radial zones, although only 3072 radial zones were effectively used on the finest resolution level of the adaptive numerical grid. A linear interpolation was used among the different zones. During the simulation all hydrodynamical quantities were interpolated with a higher order (PPM = piecewise parabolic) scheme within the grid cells. The effective number of angular zones on the finest grid level was 768. The exact use of the 2D grid was described in [18] to which we address the reader to that reference for further details. In the following we made use of the model run b18b of [19], whose explosion energy had a value of erg.
For the further analysis, the computational output in the relevant (dynamically changing) radial region was mapped to a uniform spatial polar grid of 3072 radial and 768 angular zones. From these data of the twodimensional SN model at different evolution times we extracted the matter density profiles on radial directions, having an angular separation of . On each direction the radial resolution is as fine as km at late times ( s).
In Fig. 1 we show the radial evolution of angleaveraged density profiles,
(2.1) 
at different postbounce times (1, 2 and 3 seconds). One recognizes that the SN shock wave at , while propagating outwards at supersonic speed, is followed by a dense shell delimited by two flow discontinuities, the forward shock front (at ) on the one side and an inner one, called reverse shock front at position with . The forward shock expands into the hydrostatic density profile of the progenitor star at . The reverse shock is a consequence of the fast neutrinodriven wind (a lowdensity, highentropy outflow blown off the nascent neutron star’s surface by continuous neutrino heating) colliding with the slower, shockaccelerated SN ejecta (a detailed discussion can be found in [45, 22]). The layer between both shocks contains the anisotropic and strongly perturbed neutrinoheated ejecta, which were subject to hydrodynamic instabilities before the onset of the explosion. The initial asymmetries continue to grow and trigger secondary mixing and fragmentation processes in this shell, leading to large inhomogeneities of the density distribution on all scales from the global one down to the resolution limit of the numerical simulation (see Fig. 2).
The MSW effect causes a resonant level crossing in neutrino oscillations when [46]
(2.2) 
where is the atmospheric masssquare difference, while the plus and minus sign at the righthandside refers to neutrinos and antineutrinos , respectively. ^{3}^{3}3We explicitly checked that the dependence of the SN resonant flavor conversions on the subleading “solar” squared mass gap and mixing angle is negligible for our SN model. Therefore, for normal mass hierarchy (NH, ) the resonance occurs in the channel, while for inverted mass hierarchy (IH, ) it happens for . The electron fraction in Eq. (2.2) is in the region relevant for the MSW effect.
With the horizontal band in Fig. 1 we indicate the density region corresponding to resonant neutrino oscillations for a typical SN neutrino energy range MeV and for eV (used hereafter), close to the current bestfit value [47]. It turns out that the case of sec in Fig. 1 represents the best time snapshot to have resonant flavor conversions for all the relevant neutrino energies, both across the sharp forward shock front and along its turbulent wake (back to the reverse shock front). The sec profile also happens to have the best radial resolution ( km) as compared to the others at lower or higher postbounce time (where km or higher). Therefore, in the following we will take this case as benchmark for our study.
The reason for showing first the angleaveraged radial profile is that, in our reference 2D simulations, the matter density in different directions exhibits large variations with respect to the average, which in some cases make it difficult to recognize the forward+reverseshock structure. In this context we present in Fig. 2 the twodimensional density plot for our reference SN simulation at s. The forward shock position in different directions (connected by the solid tick curve) can be rather different from the average one at km (dashed curve). Moreover, not at all angles the reverse shock is clearly visible in Fig. 2, although it emerges at km (inner dashed curve in Fig. 2) after angular averaging (see Fig. 1). Like the forward shock, the reverse shock is highly aspherical. It is strong at the head of slowly expanding, highdensity downflows, whereas it is weak in regions where highentropy, lowdensity bubbles rise faster by acceleration through buoyancy forces. The average position is determined by the massive, dense downdrafts. For later purpose, we shall take an average range of km as representative of the radial interval between the two shock waves, where turbulence is important.
In order to show the large density variations and differences radial profiles can exhibit, we plot the radial matterdensity profiles along three representative directions in Fig. 3. While the matter profiles in the direction (corresponding to ) in the left panel and at () in the right panel are rather smooth, the profile at () in the middle panel displays a ragged behavior downstream of the forward shock (at ), with many contact discontinuities due to intense turbulence. The horizontal bands represent the MSW resonant regions for the same energy range as in Fig. 1.^{4}^{4}4 We remark that the shock fronts from the SN simulations have been steepened “by hand,” in order to correct artifacts due to finite resolution. This correction is needed to properly account for strong violations of adiabaticity in neutrino flavor conversions across the shock discontinuity. We realize that, in general, multiple resonances can arise along the matterdensity profile behind the SN shock front (i.e., at ). Since the resonance pattern is quite different in these three cases, one might expect a different impact of turbulence on the conversion probabilities. However, as we shall see in Sec. 4, the impact of fluctuations is actually quite modest, although it may somewhat increase by adding further smallscale fluctuations via randomization (see Sec. 5).
3 Powerspectrum of the matter density fluctuations
While there is little to be learned from a mere collection of density profiles as in Fig. 3 –as a big variety of different cases can be observed– it turns out that the power spectra of these fluctuations show consistently the characteristic imprint of the twodimensional theory of turbulence, which predicts a brokenpowerlaw spectrum with exponents and on very general grounds (e.g., [48, 49]). In this section we thus perform a systematic statistical study of turbulence spectra along the 768 different directions of our reference SN model, which will also allow us to compare our findings with other assumptions most commonly found in the previous literature.
In order to compare the power spectra of the matter density profiles along different directions, it is useful to fix a common range over which the Fourier transforms will be evaluated. As commented in the context of Fig. 2, we assume km, corresponding to the distance between the radii of the mean forward and reverse shocks. This assumption is not crucial for our findings (as we have explicitly checked), and is adopted hereafter only for the sake of clarity and of homogeneity in the comparison of spectra at different angles. With this in mind, for each direction , we evaluated the power spectrum over the range , in which is the radius of the forward shock front along the th direction (solid curve in Fig. 2).
In the region where matter effects are relevant [see Eq. (2.2)], the neutrino oscillation length in matter is approximately given by [33]
(3.1) 
where in our numerical calculations we take a reference value , consistent with recent global data analyses [47] (and corresponding to ). Therefore, it turns out that the radial simulation resolution , the neutrino wavelength in matter , and the radial range chosen for the Fourier transform, satisfy the correct lengthscale hierarchy expected for a realistic representation of smallscale fluctuations,
(3.2) 
in the whole energy range relevant for SN neutrinos.
This hierarchy suggests that the radial resolution of our reference SN simulation is good enough to characterize the small scales relevant for neutrino oscillations in a turbulent environment, and (secondarily) that the shock wave extends over many oscillation cycles. Therefore, it makes sense to derive the spectral properties of the matter fluctuations directly from the matter density profiles, taken at face value along different directions. We anticipate, however, that the hierarchy in Eq. (3.1) is not necessarily guaranteed at large radii by the worsening of the transverse resolution; this issue will be separately discussed in Sec. 5.
For the spectral analysis, we define the density contrast
(3.3) 
where is the average density over the range , and expand it in a Fourier series:
(3.4) 
The Fourier transform
(3.5) 
is then squared to give the power spectrum:
(3.6) 
In the following, we shall refer only to wavenumber indices , since .
Figure 4 shows the values (dots) as a function of for the same three directions of Fig. 3 (corresponding to ). The points appear to be clustered along a rapidly decreasing curve, which, however, cannot be represented by a single slope in any of these three cases. Indeed, in such cases (as well as in all other directions, not shown), we find that the Fourier spectrum is quite reasonably represented by a twoslope line, i.e. by a broken powerlaw of the kind
(3.7) 
where is the overall amplitude of the powerspectrum and is the wavenumber index where the powerlaw spectrum changes the exponent from to .
For each of the 768 directions, we fit in terms of four independent parameters , , , and as in (3.7), see dashed thick lines in Fig. 4. The distributions of the fit parameters are shown in Fig. 5, where either linear or logarithmic scale was chosen in order to obtain more symmetric and “gaussian” distributions. The peak values and 1 range of the fit parameters are reported here:
(3.8) 
The spectral indices, and , of the two slopes of the broken power law are already expressed in terms of the values predicted for them by the KolmogorovKraichnan theory of twodimensional turbulence [48], namely and , respectively (see [49] for a recent review), which are indeed compatible with our results at .
Finally, Fig. 6 shows the synopsis of our fit results in terms of a broken power law with spectral indices and as in Eq. (3.7) (thick line) and 1 error ranges (dashed lines). The possible horizontal and vertical shift of the break is also shown, and takes into account the uncertainty at 1 C.L. on the determination of and of . Our findings show that 2D turbulence imprints clearly emerge with the expected spectral features from our reference SN simulation. We have found similar spectral results also for other postbounce times (not shown).
4 Neutrino crossing probabilities
In the previous section we have shown that the reference SN density profiles embed the expected turbulence patterns, and trace them with a space resolution appropriate to study flavor transitions in matter. We now analyze such transitions in detail.
The main theoretical quantity associated to neutrino flavor change in SNe is the socalled crossing probability , taking the same form for neutrinos and antineutrinos (see [14] and refs. therein)
(4.1) 
The crossing probability parameterizes the deviations with respect to a pure adiabatic neutrino flavor evolution in matter [46]. Indeed, violations of the adiabaticity would imply possible level crossings among the instantaneous neutrino mass eigenstates in matter. These violations could be relevant especially at the resonance points [Eq. (2.2)]. For a nonmonotonic matter density profile like the one encountered by neutrinos in SN, multiple level crossing can occur. In this case, one expects that the strongest violation of adiabaticity would occur when neutrinos propagate across the shockfront discontinuity. In this case, the crossing probability . Conversely a purely adiabatic propagation would correspond to . The survival probability of SN neutrinos at Earth (neglecting Earth matter crossing) is related to the crossing probability as follows:
(4.2) 
From the above equations, it appears that can modulate the (otherwise constant) survival probability of in normal hierarchy and of in inverted hierarchy, thus providing an important handle to solve the current hierarchy ambiguity.
In order to compute , we have numerically integrated flavor evolution equations along all the available matter density profiles, although we shall focus only on the representative ones shown in Fig. 3, for the sake of simplicity. We closely follow the approach of Ref. [14] to compute (the interested reader is referred to this reference for further details). In particular, since matter effects are dominant at large radii ( km) we neglect nonlinear effects associated with neutrinoneutrino interactions that occur at smaller radii.
In Fig. 7 we depict the crossing probability for the three chosen representative cases, as a function of the neutrino energy (gray curves). In all cases, exhibits the typical tophat structure widely discussed in previous literature (as quoted, e.g., in [14]). In particular, rapidly jumps from zero to for MeV. The energy range where corresponds to resonances occurring within the forward shock front, where extremely nonadiabatic transitions take place. Indeed, as shown in [14], for values of as large as the measured one, the crossing probabilities occurring before the shock front and along the static profile ahead of the shock are exponentially suppressed. Instead, the energy width of the tophat function depends on the height of the shock front relative to the resonance band in Fig. 3. In the left and central panels the resonance band intersects the forward shock front for all the range MeV. Therefore, in all this energy range the resonance occurring within the shock front gives , while the value of drops to zero outside this energy range. Conversely, in the right panel the energy range of is smaller, since the resonance occurs on the static profile already at MeV, where drops to zero.
We also realize that, while in the right panel the function has a smooth profile, in the other panels it presents an oscillatory behavior in the energy range around MeV. This oscillatory behavior is due to the phase effect caused by the interference between multiple narrowlyspaced resonances [50] corresponding to a rapidly and largely fluctuating turbulent profile (see Fig. 3). In practice, the observability of this oscillatory pattern is severely limited by the inability of the detectors to reconstruct the neutrino energy precisely [50], leading to a smeared profile. For illustration, the black curves in Fig. 7 are obtained by convoluting the gray ones with a “tophat” energy resolution function having MeV width, so as to enforce an effective phase averaging (smearing).
Independently of the details, we find that can be sizable in all the three panels of Fig. 7, and is not significantly suppressed in those cases where fluctuations are evident (left and middle panels) as compared with the case having a smoother matter profile (right panel). Moreover, although the three corresponding matter density profiles are associated with power spectra having a break at rather different (see Fig. 4), the global behavior of is rather similar, i.e., there is no evident link between the break in the power spectrum and the behavior of . Indeed, what is mostly relevant for the shape of is the presence of a crossing (resonance) within the forward shock. This fact is mostly related to the height of the shock front rather than to the presence of the turbulence. Indeed, as already noticed in [40], the “large” measured value of makes the resonances more adiabatic downstream of the forward shock front. Then, the main effect of turbulence is the phase effect which increases with the amplitude of the fluctuations. In addition, in none of the cases considered the turbulence has enough strength to strongly “damp” the profile of the . Large damping (with ) was found in [33] where “deltacorrelated” fluctuations were considered, with correlation scale of km and a fluctuation amplitude of a few . Such a hypothetical amplitude is actually much larger than what our simulation gives on the same length scale. Indeed, the density contrast on a scale km, is given by for our best fit values [see Eq. (3.8)].
In [40] a Kolmogorovlike fluctuation powerspectrum (characterized by a single exponent in 3D turbulence theory) was included through the “randomization method” [38] on top of an undisturbed density profile. A parametric study was performed for different values of the overall amplitude of the fluctuation powerspectrum, including cases with “small” values (of the same order of those emerging in our reference SN model). In particular, for a powerspectrum amplitude of , comparable with the amplitude from our simulations (see Fig. 3), no significant damping was found, while phase effects were still present. We qualitatively agree with such findings, but the intrinsic differences between our SN models (plus the fact that collective and matter effects are combined in [40]) prevent a more detailed comparison. In any case, we remark that our results are directly linked to a simulation rather than to a parameterization.
5 Reconstructing small scale fluctuations with the randomization method
In the previous Section the crossing probability has been evaluated by integrating the neutrino flavor evolution equations along same representative density profiles, taken directly from our reference 2D simulations. As discussed in Sec. 2, for the time snapshot we are considering ( s) these simulations have the best radial resolution km). However, given the angular separation of among two radial consecutive trajectories, it follows that the transverse resolution is given by
(5.1) 
Therefore, at the relatively large radii characterizing the shock front (average value km), where flavor conversions largely develop, the transverse resolution can be km, which is definitely worse than the radial one.
Since turbulence develops both in radial and transverse directions, this fact implies that the simulations might not really resolve length scales below km at typical shockfront radii. From the viewpoint of the power spectrum , this length scale corresponds to a multipoles number . In other words, the spectrum features for might not be correctly captured by the available simulations.
We stress that this limitation does not invalidate the statistical analysis of the power spectra in Sec. 3. Indeed, in the vast majority of analyzed cases the spectrum “knee” (i.e., the multipole ) is found below (see the upper right panel of Fig. 5). However, an overall resolution of O(100) km can be of the same order or larger than typical oscillation lengths, and may invalidate the hierarchy in Eq. (3.1), with possible alterations of the profiles discussed in Sec. 4.
In order to overcome this limitation and to investigate the effect of small scales not directly resolved by the simulations, it is necessary to rely on some procedure to “continue” the powerspectrum for multipoles larger than . At this regard, in the literature it has been proposed to generate random fluctuations with a given power spectrum by means of the “randomization method” [38]. In [37, 39, 40] this method was applied to generate fluctuations over smooth 1D density profiles, assuming a Kolmogorov spectrum with tunable amplitudes. Herein, we investigate a variant of this randomization procedure. We take the previous (Sec. 3) fluctuation spectra at face value for , and randomize them only for ; for definiteness we take , although the precise number chosen for is of little relevance. We also anchor the smallscale randomization to largescale spectral features as follows.
After truncating the power spectra at , we have studied the distribution of the amplitudes [see Eq. (3.6)] along each direction. In particular, we have considered their residuals with respect to the values expected from the bestfit, broken powerlaw spectrum (as represented, e.g., by the continuous lines in Fig. 4). As it can be seen at a glance in Fig. 4 (as well as for other profiles, not shown), the typical scatter of residuals is relatively similar in different profiles (in logarithmic scale). We found that the distribution of the fractional logresiduals, defined as
(5.2) 
is approximately gaussian in each of the 768 truncated spectra, with central value always close to zero, and variance quite similar for all the trajectories, within a factor less than two.
A conservative estimate for the overall variance of the residuals is , for . We then assume that the values for are also scattered around the bestfit KraichnanKolmogorov spectrum with this same overall variance, and with a random phase . More precisely, for we randomly generate amplitudes as
(5.3) 
where the are drawn from a gaussian distribution with null average and variance , and the phases are drawn from a flat distribution in . Then the density contrast in each direction [Eq. (3.4)] is derived as
(5.4) 
where for the coefficients are taken directly from the simulations , while they are randomized via Eq. (5.3) for . Of course, one can generate different realizations of the density profiles for (see below). Finally, we adopt a filter in order to suppress spurios localized fluctuations (spikes) in the reconstructed density profiles, which may emerge in singular points where the original profile is steplike.
The density profile in the central panel of Fig. 3 is particularly well suited to investigate the effects of the above randomization procedure, since it displays large fluctuations in the interesting energy region and for a large radial range. Conversely, in the left (right) panel of Fig. 3, level crossings occur only at small radii (energies), where randomization effects produce only a minor impact on (not shown).
Focusing on the profile in the central panel of Fig. 3, we have calculated the crossing probability for 100 realizations of the randomized density fluctuations above . We apply energy smearing as in Fig. 7 to remove (unobservable) phase effects in each realization. Figure 8 shows the average profile from the 100 realizations (dashed curve), together with its dispersion (gray band). For comparison, the profile in Fig. 7 is superposed (thick black curve). It appears that in the energy range MeV the effect of (randomized) small scale fluctuations can generate differences up to in the , but does not radically change the profile. In particular, the doublepeak shape of the is similar to the one obtained from the original density profile, taken at face value from the simulations. Moreover, the results from this profile are roughly in agreement within the dispersion band of the randomized profile. In general, by applying the randomization method described above also to other density profiles from our 2D simulations, we do not find dramatic changes (e.g., strong or complete depolarization) as in other cases considered in the previous literature [33, 37, 39, 40]. Although limited to 2D SN models so far, our results show the importance of anchoring studies of turbulence effect and randomization procedures, whenever possible, to highresolution SN simulations.
6 Summary and prospects
In this work we have carried out a study of the effects of matter density fluctuations on the conversion probabilities of SN neutrinos, characterizing the matter turbulence from a highresolution 2D supernova model. The latter was evolved to sufficiently late times to follow the shock dynamics in the density regime relevant for flavor transformations at energies typical of supernova neutrinos. We found that the power spectra of the matter density fluctuations along different radial directions show a “typical” structure which represents the imprint of generic 2D turbulence, namely a broken powerlaw spectrum with average exponents and . We evaluated the flavor evolution of SN neutrinos along representative density profiles, exhibiting powerspectra with breaks at different wave numbers. We found that the crossing probability exhibits, in all the different cases, a tophat structure where jumps to at energies where an extremely nonadiabatic resonance takes place across the shock front. We do not find significant damping of the crossing probabilities associated with the matter turbulence. Since the available simulations might not fully resolve details at length scales below O(100) km, we have also improved the treatment of smallscale fluctuations by means of a randomization procedure anchored to larger scales. We find that the randomized profiles generally preserve the main features of the profile, although with a slightly more pronounced damping in some cases. As a result, the shockwave imprint on the crossing probability seems to be a possible observable signature in the future corecollapse SN events, with no large suppression by turbulence. We emphasize that, if SN neutrino energy spectra of different flavor are not too similar at late times [54], shockwave effects would give an observable signature in the SN neutrino signal at large underground detectors [15, 16]. This would be an important tool to follow in realtime the SN explosion dynamics, as well as to determine the neutrino mass hierarchy. Therefore, we hope that the potential importance of shockwave signatures as a probe of neutrino physics and astrophysics will motivate further attention and progress in numerical simulations, in order to investigate the impact of matter turbulence in other independent SN models, possibly also in 3D.
As already commented, the transition from 2D to 3D is expected to bring turbulence spectra from a broken power law into a single power law, with a Kolmogorov exponent across all scales. This expectation (not yet confirmed by the available 3D results [42, 43]) would in principle bring more power to fluctuations at small scales, with a possible enhancement of the damping effect. In order to clarify this issue, a detailed study with highresolution threedimensional SN models is mandatory. In this context, another delicate issue is represented by the poorly known asymmetries in the progenitor stars prior to collapse. While all current SN simulations are started from spherical models, it is possible that the progenitor core at the onset of gravitational instability has turbulent density perturbations on different scales (e.g., [51] and references therein). These precollapse inhomogeneities and fluctuations could affect the SN explosion by seeding and enhancing the hydrodynamic instabilities in the postshock region [52], and their effects on MSW neutrino flavor conversions still need to be explored. Finally, it should be noted that, since the neutrinosphere source has a finite size, the neutrinos emitted from different points propagate through somewhat different (but correlated) density profiles. In this context, it might also be interesting to study the correlation of the transition probabilities along relatively close trajectories, as proposed in [53].
In conclusion, we believe that our study offers an interesting new perspective to the issues raised by SN neutrino turbulence, at least in the realm of 2D simulations. Although our approach has some limitations which prevent a wide generalizatios of the results, especially to the 3D case, it illustrates a methodology that can be usefully applied also to other SN models. When future, highresolution simulations will become available for detailed and comparative analyses, further studies will test the robustness of our results.
Acknowledgements
We warmly thank Konstantinos Kifonidis for providing the data of his 2D simulations and Timur Rashba, Lab Saha, Irene Tamborra, and Ricard Tomàs for useful discussions and help on various topics related to this work. The work of E.B. and A.M. was supported by the German Science Foundation (DFG) within the Collaborative Research Center 676 “Particles, Strings and the Early Universe”. The work of E.L. was partially supported by the Italian MIUR and Istituto Nazionale di Fisica Nucleare (INFN) through the “Theoretical Astroparticle Physics” project. H.T.J. acknowledges support by the Deutsche Forschungsgemeinschaft through the the Transregional Collaborative Research Center SFB/TR7 “Gravitational Wave Astronomy” and the Cluster of Excellence EXC 153 “Origin and Structure of the Universe.” S.C. acknowledges support from the European Union through a Marie Curie Fellowship, Grant No. PIIFGA2011299861, and through the ITN “Invisibles”, Grant No. PITNGA2011289442.
References
 [1] J. Pantaleone, “Neutrino oscillations at high densities,” Phys. Lett. B 287, 128 (1992).
 [2] Y. Z. Qian and G. M. Fuller, “Neutrinoneutrino scattering and matter enhanced neutrino flavor transformation in Supernovae,” Phys. Rev. D 51, 1479 (1995) [astroph/9406073].
 [3] H. Duan, G. M. Fuller and Y. Z. Qian, “Collective Neutrino Flavor Transformation In Supernovae,” Phys. Rev. D 74, 123004 (2006) [astroph/0511275].
 [4] H. Duan, G. M. Fuller, J. Carlson and Y. Z. Qian, “Simulation of coherent nonlinear neutrino flavor transformation in the supernova environment. I: Correlated neutrino trajectories,” Phys. Rev. D 74, 105014 (2006) [astroph/0606616].
 [5] S. Hannestad, G. G. Raffelt, G. Sigl and Y. Y. Y. Wong, “Selfinduced conversion in dense neutrino gases: Pendulum in flavour space,” Phys. Rev. D 74, 105010 (2006) [Erratumibid. D 76, 029901 (2007)] [astroph/0608695].
 [6] G. L. Fogli, E. Lisi, A. Marrone and A. Mirizzi, “Collective neutrino flavor transitions in supernovae and the role of trajectory averaging,” JCAP 0712, 010 (2007) [arXiv:0707.1998 [hepph]].
 [7] H. Duan, G. M. Fuller and Y. Z. Qian, “Collective Neutrino Oscillations,” Ann. Rev. Nucl. Part. Sci. 60, 569 (2010) [arXiv:1001.2799 [hepph]].
 [8] H. Duan, G. M. Fuller, J. Carlson and Y. Q. Zhong, “Neutrino Mass Hierarchy and Stepwise Spectral Swapping of Supernova Neutrino Flavors,” Phys. Rev. Lett. 99, 241802 (2007) [arXiv:0707.0290 [astroph]].
 [9] B. Dasgupta, A. Dighe, G. G. Raffelt and A. Y. .Smirnov, “Multiple Spectral Splits of Supernova Neutrinos,” Phys. Rev. Lett. 103, 051105 (2009) [arXiv:0904.3542 [hepph]].
 [10] A. Mirizzi, “Selfinduced spectral splits with multiazimuthalangle effects for different supernova neutrino fluxes,” arXiv:1308.5255 [hepph].
 [11] L. Wolfenstein, “Neutrino Oscillations In Matter,” Phys. Rev. D 17, 2369 (1978); S. P. Mikheev and A. Yu. Smirnov, “Resonance Enhancement Of Oscillations In Matter And Solar Neutrino Spectroscopy,” Yad. Fiz. 42, 1441 (1985) [Sov. J. Nucl. Phys. 42, 913 (1985)].
 [12] A. S. Dighe and A. Y. Smirnov, “Identifying the neutrino mass spectrum from the neutrino burst from a supernova,” Phys. Rev. D 62, 033007 (2000) [hepph/9907423].
 [13] R. C. Schirato and G. M. Fuller, “Connection between supernova shocks, flavor transformation, and the neutrino signal,” astroph/0205390.
 [14] G. L. Fogli, E. Lisi, D. Montanino and A. Mirizzi, “Analysis of energy and time dependence of supernova shock effects on neutrino crossing probabilities,” Phys. Rev. D 68, 033005 (2003) [hepph/0304056].
 [15] G. L. Fogli, E. Lisi, A. Mirizzi and D. Montanino, “Probing supernova shock waves and neutrino flavor transitions in nextgeneration waterCerenkov detectors,” JCAP 0504, 002 (2005) [hepph/0412046].
 [16] R. Tomàs, M. Kachelriess, G. Raffelt, A. Dighe, H. T. Janka and L. Scheck, “Neutrino signatures of supernova shock and reverse shock propagation,” JCAP 0409, 015 (2004) [astroph/0407132].
 [17] J. Gava, J. Kneller, C. Volpe and G. C. McLaughlin, “A Dynamical collective calculation of supernova neutrino signals,” Phys. Rev. Lett. 103, 071101 (2009) [arXiv:0902.0317 [hepph]].
 [18] K. Kifonidis, T. Plewa, H.T. Janka and E. Müller, “Nonspherical core collapse supernovae. I. Neutrinodriven convection, RayleighTaylor instabilities, and the formation and propagation of metal clumps,” Astron. Astrophys. 408, 621 (2003) [astroph/0302239].
 [19] K. Kifonidis, T. Plewa, L. Scheck, H.T. Janka and E. Müller, “Nonspherical core collapse supernovae. II. The latetime evolution of globally anisotropic neutrinodriven explosions and their implications for SN 1987 A,” Astron. Astrophys. 453, 661 (2006) [astroph/0511369].
 [20] L. Scheck, K. Kifonidis, H.T. Janka and E. Müller, “Multidimensional supernova simulations with approximative neutrino transport. I. Neutron star kicks and the anisotropy of neutrinodriven explosions in two spatial dimensions,” Astron. Astrophys. 457, 963 (2006) [astroph/0601302].
 [21] N. Hammer, H.T. Janka, and E. Müller, “Threedimensional simulations of mixing instabilities in supernova explosions,” Astrophys. J. 714, 1371 (2010) [astroph/0908.3474].
 [22] A. Arcones and H.T. Janka, “Nucleosynthesisrelevant conditions in neutrinodriven supernova outflows. II. The reverse shock in twodimensional simulations,” Astron. Astrophys. 526, 160 (2011) [astroph/1008.0882].
 [23] E. Müller, H.T. Janka and A. Wongwathanarat, “Parametrized 3D models of neutrinodriven supernova explosions. Neutrino emission asymmetries and gravitationalwave signals,” Astron. Astrophys. 537, 63 (2012) [astroph/1106.6301].
 [24] A. Wongwathanarat, H.T. Janka, and E. Müller, “Threedimensional neutrinodriven supernovae: Neutron star kicks, spins, and asymmetric ejection of nucleosynthesis products,” Astron. Astrophys. 552, 126 (2013) [astroph/1210.8148].
 [25] A. Schaefer and S. E. Koonin, “Influence of Density Fluctuations on Solar Neutrino Conversion,” Phys. Lett. B 185, 417 (1987).
 [26] R. F. Sawyer, “Neutrino oscillations in inhomogeneous matter,” Phys. Rev. D 42, 3908 (1990).
 [27] F. N. Loreti and A. B. Balantekin, “Neutrino oscillations in noisy media,” Phys. Rev. D 50, 4762 (1994) [nuclth/9406003].
 [28] H. Nunokawa, A. Rossi, V. B. Semikoz and J. W. F. Valle, “The Effect of random matter density perturbations on the MSW solution to the solar neutrino problem,” Nucl. Phys. B 472, 495 (1996) [hepph/9602307].
 [29] A. B. Balantekin, J. M. Fetter and F. N. Loreti, “The MSW effect in a fluctuating matter density,” Phys. Rev. D 54, 3941 (1996) [astroph/9604061].
 [30] C. P. Burgess and D. Michaud, “Neutrino propagation in a fluctuating sun,” Annals Phys. 256, 1 (1997) [hepph/9606295].
 [31] E. TorrenteLujan, “Finite dimensional systems with random external fields and neutrino propagation in fluctuating media,” Phys. Rev. D 59, 073001 (1999) [hepph/9807361].
 [32] F. N. Loreti, Y. Z. Qian, G. M. Fuller and A. B. Balantekin, “Effects of random density fluctuations on matter enhanced neutrino flavor transitions in supernovae and implications for supernova dynamics and nucleosynthesis,” Phys. Rev. D 52, 6664 (1995) [astroph/9508106].
 [33] G. L. Fogli, E. Lisi, A. Mirizzi and D. Montanino, “Damping of supernova neutrino transitions in stochastic shockwave density profiles,” JCAP 0606, 012 (2006) [hepph/0603033].
 [34] S. Choubey, N. P. Harries and G. G. Ross, “Turbulent supernova shock waves and the sterile neutrino signature in megaton water detectors,” Phys. Rev. D 76, 073013 (2007) [hepph/0703092 [HEPPH]].
 [35] F. Benatti and R. Floreanini, “Dissipative neutrino oscillations in randomly fluctuating matter,” Phys. Rev. D 71, 013003 (2005) [hepph/0412311].
 [36] A. Friedland and A. Gruzinov, “Neutrino signatures of supernova turbulence,” astroph/0607244.
 [37] J. P. Kneller and C. Volpe, “Turbulence effects on supernova neutrinos,” Phys. Rev. D 82, 123004 (2010) [arXiv:1006.0913 [hepph]].
 [38] A. J. Majda and P. R. Kramer, “Simplified models for turbulent diffusion: Theory, numerical modelling, and physical phenomena,” Phys. Rep. 314 237 (1999).
 [39] J. P. Kneller and A. W. Mauney, “The consequences of large for the turbulence signatures in supernova neutrinos,” arXiv:1302.3825 [hepph].
 [40] T. Lund and J. P. Kneller, “Combining collective, MSW, and turbulence effects in supernova neutrino flavor evolution,” Phys. Rev. D 88, 023008 (2013) [arXiv:1304.6372 [astroph.HE]].
 [41] T. Rashba, “Effects of random density fluctuations on Supernovae and solar neutrinos.” Talk at the Workshop on Next generation Nucleon decay and Neutrino detectors, NNN 2006, Seattle 21–23 September 2006. Slides available at http://neutrino.phys.washington.edu/nnn06/slides/ rashba_supernovae_density_fluctuations_v2.pdf.
 [42] J.C. Dolence, A. Burrows, J.W. Murphy and J. Nordhaus, “Dimensional dependence of the hydrodynamics of corecollapse supernovae,” Astrophys. J. 765, 110 (2013) [arXiv:1210.5241].
 [43] S.M. Couch and E.P. O’Connor, “Highresolution threedimensional simulations of corecollapse supernovae in multiple progenitors,” [arXiv:1310.5728].
 [44] S. E. Woosley, “SN1987A: After the peak,” Astrophys. J. 330, 218 (1988).
 [45] A. Arcones, H.T. Janka and L. Scheck, “Nucleosynthesisrelevant conditions in neutrinodriven supernova outflows. Spherically symmetric hydrodynamic simulations,” Astron. Astrophys. 467, 1227 (2007) [astroph/0612582].
 [46] T. K. Kuo and J. T. Pantaleone, “Neutrino Oscillations in Matter,” Rev. Mod. Phys. 61, 937 (1989).
 [47] G. L. Fogli, E. Lisi, A. Marrone, D. Montanino, A. Palazzo and A. M. Rotunno, “Global analysis of neutrino masses, mixings and phases: entering the era of leptonic CP violation searches,” Phys. Rev. D 86, 013012 (2012) [arXiv:1205.5254 [hepph]].
 [48] R. Kraichnan, “Inertial ranges in twodimensional turbulence,” Phys. Fluids 10, 1417 (1967); R. Kraichnan, “Inertialrange transfer in two and threedimensional turbulence,” Journ. Fluid Mech. 47, 525 (1971); R. Kraichnan and D. Montgomery, “Twodimensional turbulence,” Rep. Prog. Phys. 43, 547 (1980).
 [49] G. Boffetta and R. E. Ecke, “TwoDimensional Turbulence,” Annual Review of Fluid Mechanics 44, 427 (2012).
 [50] B. Dasgupta and A. Dighe, “Phase effects in neutrino conversions during a supernova shock wave,” Phys. Rev. D 75, 093002 (2007) [hepph/0510219].
 [51] W.D. Arnett and C. Meakin, “Toward realistic progenitors of corecollapse supernovae,” Astrophys. J. 733, 78 (2011) [arXiv:1101.5646].
 [52] S.M. Couch and C.D. Ott, “Revival of the stalled corecollapse supernova shock triggered by precollapse asphericity in the progenitor star,” [arXiv:1309.2632].
 [53] J. P. Kneller and A. W. Mauney, “Does the finite size of the protoneutron star preclude supernova neutrino flavor scintillation due to turbulence?,” Phys. Rev. D 88, no. 4, 045020 (2013).
 [54] T. Fischer, S. C. Whitehouse, A. Mezzacappa, F. K. Thielemann and M. Liebendorfer, “Protoneutron star evolution and the neutrino driven wind in general relativistic neutrino radiation hydrodynamics simulations,” Astron. Astrophys. 517, A80 (2010) [arXiv:0908.1871 [astroph.HE]].