On probing turbulence in core-collapse supernovae in upcoming neutrino detectors

Neutrino propagation through a turbulent medium can be highly non-adiabatic leading to distinct signatures in the survival probabilities. A core-collapse supernova can be host to a number of hydrodynamic instabilities which occur behind the shockfront. Such instabilities between the forward shock and a possible reverse shock can lead to cascades introducing turbulence in the associated matter profile, which can imprint itself in the neutrino signal. In this work, we consider realistic matter profiles and seed in the turbulence using a randomization scheme to study its effects on neutrino propagation in an effective two-flavor framework. We focus on the potential of upcoming neutrino detectors — DUNE and Hyper-Kamiokande to constrain the parameters characterizing turbulence in a supernova. We find that these experiments can effectively constrain the parameter space for the amplitude of the spectra, they will only have mild sensitivity to the spectral index, and cannot inform on deviations from the usual Kolmogorov 5/3 inverse power law. Furthermore, we also confirm that the double-dip feature, originally predicted in the neutrino spectra associated with forward and reverse shocks, can be completely washed away in the presence of turbulence, leading to total flavor depolarization.


I. INTRODUCTION
The detection of neutrinos from SN1987A can be regarded as one of the first instances of multi-messenger astronomy, where neutrino events were observed hours before visible light.Although limited in statistics -a total of O(30) neutrino events in the Kamiokande detector [1,2], the IMB detector [3] and the Baksan detector [4] were observed -this data has been used widely by the particle and astroparticle physics community to probe important questions of fundamental physics.
It is already well-known that neutrinos from a galactic core-collapse supernova (SN) will be our best bet to probe the supernova dynamics.These neutrinos originate from the deepest regions of the SN and propagate through the entire matter envelope before exiting the SN and arriving at the Earth.The flavor composition of the neutrinos is extremely sensitive to the background fermion density it encounters en route.Close to the neutrinosphere, the neutrino density is large enough, so neutrino forward scattering on a neutrino background is the dominant contribution to the propagation Hamiltonian.This can lead to collective flavor conversion -both "slow" [5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24] and ''fast"  (see [66][67][68][69] for detailed reviews on the subject)-which is the major mechanism of flavor transformation within radii of ∼ O(100) km from the neutrinosphere.
At larger radii, neutrino density falls off and matter effects, due to coherent forward scattering on background electrons, become important.Multiple studies have demonstrated the sensitivity of the neutrino spectra to the underlying SN density profile through the Mikheyev-Smirnov-Wolfenstein (MSW) resonances [70,71], including signatures of the shockwave propagation at late times [72][73][74][75][76][77][78][79].In the absence of the shockwave, neutrino propagation is shown to be adiabatic, governed only by the vacuum mixing angle and the matter density.However, violations of adiabaticity can arise due to an abrupt change in density at the shockfront.This non-adiabaticity leaves its imprint on the flavor transition and can be used to reconstruct a real-time propagation of the shockwave.
A distinctive signature of the non-adiabaticity introduced by the shockwave propagation can manifest in the form of dips in the energy-dependent survival probabilities of neutrinos.A few seconds post the SN core bounce, a possible reverse shockwave can arise if a supersonic neutrino-driven wind collides with the ejecta, as shown in Fig. 1.Note that this is very sensitive to whether the neutrino-driven wind is supersonic or subsonic, for a detailed discussion, see [80].In [77], it was demonstrated that for neutrinos passing through the two density discontinuities, characterised by the reverse and the forward shocks, a "double-dip" feature would be present in the neutrino survival probability as a function of energy.Such a signal can be present in electron neutrinos for normal mass ordering (NMO) and electron antineutrinos for inverted mass ordering (IMO) and can be used to trace the position of the two shockfronts.
However, this is not the complete story.A SN progenitor is host to a variety of hydrodynamic instabilities, which take place behind the shockfront.This leads to turbulence in the ejecta, causing stochastic fluctuations in the matter density.Turbulence is generally expected to be present in aspherical simulations of SNe, and can play an important role in seeding the explosion [81].These stochastic fluctuations in the matter density have a non-trivial impact on neutrino flavor evolution.
Neutrino flavor evolution in a turbulent medium has been studied extensively to determine its impact on the neutrino spectra from the next Galactic supernova [82][83][84][85][86][87][88][89][90][91]. .Generically, the effect of turbulence on flavor propagation can be broadly classified into two types: (i) 'stimulated' or non-adiabatic transitions between different propagation eigenstates if the stochastic density fluctuations are present near the MSW resonance region [84], and (ii) phase effects due to the presence of semi-adiabatic resonances [78].The latter effect is more subtle and can be suppressed if the amplitude of fluctuations is large.
It is numerically challenging to obtain realistic turbulent matter density profiles from aspherical (3D) hydrodynamic simulations.Theoretically, the effect of turbulence has been modelled by adding disturbances to a smooth matter profile obtained from coarse, spherically symmetric simulations.This can be done in a number of different ways -for example, by adding delta correlated noise [82], or by modelling the turbulent density fluctuations through a Kolmogorov-type power spectrum [83].In particular, recent 2D hydrodynamic simulations have shown the presence of a power-law dependence of the power spectrum on the wavelength of fluctuations.The long-wavelength fluctuation indeed shows a power-law dependence with the mean spectral index given by p = 1.85 +0.54  −0.77 , thus indicating that 2D turbulence might be closer to having a Kolmogorov-type power spectrum [85].However, simulations also indicate that the turbulence might not be fully developed, hence other values of the spectral indices cannot be ruled out.Moreover, it is well known that the cascading of turbulence in 3D is different from 2D [92,93].However, the lack of sufficient resolution makes it difficult to recover the Kolmogorov power law behaviour from 3D simulations.
A theoretical study of the impact of turbulence on the neutrino survival probability was performed in [83], where power-law fluctuations in the density profile were assumed.It was found that the electron neutrino survival probability exhibits a power law when the amplitude of turbulence is increased and a corresponding analytical estimate to explain the power law for the simple turbulence model was derived.This was followed by a series of works [88][89][90], where the authors used rotating wave approximation to predict the final survival probability amongst other techniques both numerical and analytical.The latter formalism was also used to predict the dependence of the neutrino flavor transition on the amplitude and the spectral index of turbulence.However, various unanswered questions remain.For example, does the "double-dip" feature pointed out in [77] survive in the presence of turbulence?Can turbulence inside a SN media follow a different power-law index than the Kolmogorov prediction, and with what amplitude?Can future-generation neutrino experiments designed to detect thousands of neutrinos from the next galactic supernova actually distinguish between different parameters characterising the turbulence in a supernova?
We aim to shed some light on these questions in this study through a detailed effective two-flavor study of neutrino propagation through a typical SN media profile, characteristic of the cooling phase.We add turbulence to the smooth density profile following the prescription in [86].We find that the addition of turbulence with an amplitude beyond 10% washes out the double-dip feature completely.We simulate neutrino events in the upcoming DUNE and Hyper-Kamiokande experiment and find that they can distinguish between a turbulent and a non-turbulent SN matter profile with high significance.On the other hand, they do not have enough sensitivity to probe the index of the power spectrum.
Our paper is organized as follows.We discuss turbulence in supernovae and the slab-approximation used as the numerical technique in Sec.II.The algorithm to seed the turbulence along with its impact on neutrino propagation is discussed in Sec.III.The effects of turbulence on the double-dip feature in the neutrino spectra are studied in Sec.III C. In Sec.IV, we discuss the event rates in DUNE and Hyper-Kamiokande along with their potential to constrain the physical parameters associated with turbulence for a galactic scale supernova.We conclude and discuss the implications of our work in Sec.V.The numerical algorithm developed for this work is summarized in Appendix.A.

II. TURBULENCE IN SUPERNOVAE
We work in the approximate framework of 2-neutrino oscillations, where we have the electron neutrinos (ν e ) and ν x where, x = {µ, τ }, along with each of their antineutrino components.The Equation of Motion (EoM) for a neutrino flavor |ν α ⟩ can be written as where H vac = diag(∆m 2 /(2E), −∆m 2 /(2E)) is the vacuum Hamiltonian and U is the 2 × 2 unitary mixing matrix.The matter potential V (r) has a contribution from the forward scattering of neutrinos off the background fermions and is given by V 0 (r) = √ 2G F n e (r)diag (1,0), where G F is the Fermi constant and n e (r) is the electron number density at a given distance r.This is extracted from a snapshot of the matter density at a certain time from a simulation.We further neglect the contributions from the non-electron flavor neutrinos, which cancel at the tree-level and arise only at the loop-level [94][95][96].Neutrinos close to the neutrinosphere also experience a forward scattering potential due to their scattering on other background neutrinos [97].This potential is given by V νν ≃ √ 2G F n ν , where n ν is the total background neutrino density.This neutrino-neutrino potential dominates deep inside the SN, and can lead to fascinating collective effects.However, in the region of our interest, close to the shockfront, this neutrino potential is subdominant, and collective oscillations can be safely ignored.
The turbulence is included in our model by adding a random field F (r) as The random field F (r) is defined in the following way where C * is the amplitude of the fluctuations, F rand (r) is a real-valued homogenous Gaussian scalar random field.The reverse shock position is given by r rs while the forward shock position is r fs .The tanh terms suppress the fluctuations close to the shocks and prevent discontinuities at r rs and r fs .The parameter λ is the length scale over which the fluctuations reach their extent.We choose λ = 100 km.The standard way of calculating the evolution of the neutrino wavefunctions involves numerically solving the Schrödinger equation directly.Unfortunately for irregular potentials (matter densities), this method is not efficient.In particular, the current scenario has turbulences making the matter potential highly noisy and the standard scenario for solving such problems is unsuitable.Hence, we use the slab approximation [98] to calculate the spatial evolution of the neutrino wavefunctions.
In this approximation, the density profile is approximated by a series of slabs of constant density.We discretize the spatial dimension of physical size L into n slabs, thus the width of each slab is, ∆x = L/n.In such a discretization scheme, the starting point is given by x 0 and the boundaries of the slabs are at x 1 , x 2 , . . ., x n .The mass eigenstates of the neutrinos propagate as plane waves in constant density regions with phases given by exp(±i∆m 2 M ∆x/4E), where the effective mass-squared difference in matter is given by with θ being the vacuum mixing angle.A CC is the effective effective charge current potential, All the information about the matter density profile is provided by this term.The wavefunction in the flavor basis is given by, The wavefunction at the boundaries is joined according to the following scheme Ψ e (x 0 ) , (7) where, • the unitary evolution operator in each slab diagonal basis U M (∆x) is given by, • U M is the unitary mixing matrix in matter given by where the effective mixing angle in matter • The notation [. . .] (i) indicates that all the matter-dependent quantities inside the brackets need to be evaluated with the matter density profile of the i-th slab, where the slab goes from x i−1 to x i .
The wavefunction in the matter basis can be easily obtained using ( 7) and ( 9) as The probability that a neutrino of flavor α is found in some flavor β after passing through the matter profile is given by Similarly, the corresponding survival probability of the matter eigenstate is given by The adiabaticity parameter, which is an important quantity to understand the changes in oscillation and transition probabilities related to neutrino propagation, is defined as [98] Armed with this, we can investigate the propagation of neutrinos in an irregular potential.

III. TURBULENCE: NUMERICAL ANSATZ
In this section, we discuss the Randomization Method used to implement stochastic fluctuations following a powerlaw on top of a smooth matter profile.A. The Randomization Method: Generating the turbulence In previous works, the authors of [84] used the Randomization Method [99] to generate and sample the turbulence.In [86], a slightly advanced version of the randomization method -'Variant C' was used.For this work, we will use the Variant C of the randomization method to generate the turbulence.It is important to note here that in this method the modes are distributed in a logarithmic scale which makes the distribution scale-invariant by ensuring that the realizations are uniform over the decades of length scales considered in the work.We summarize the method below.
Let the wavenumber space from which we draw the wavenumbers have a size ∆, such that, ∆ = [0, ∞).For numerical implementation, the space ∆ is divided into n non-intersecting intervals ∆ j , such that, ∆ = ∪ n j=1 ∆ j .We then select n 0 independent random wavenumbers from each of these intervals ∆ j , which can now be indexed as k jl , where, l = 1, 2, . . ., n 0 .The selection of the k jl 's are according to the probability distribution function p j (k) defined as where E(k) is the power spectrum, which we assume to be of the following form, In the above definition of the power spectrum, α is the spectral index, k * = 2π/λ * is the wavenumber cutoff, i.e., the cutoff that decides the longest non-zero turbulence wavelength (λ * ) and hence the smallest wavenumber.Thus, it can be referred to as the driving scale for our turbulence model [87].For this work, we treat the spectral index α as a free parameter that is varied and the corresponding sensitivity of the next-generation neutrino detectors is tested.We choose, λ * = 2(r fs − r rs ), which is twice the distance between the positions of the reverse and forward shocks.In Eq. ( 14), σ 2 j = 2 ∆j E(k)dk.We can now define F rand (r) (see Eqs. ( 2) and ( 3)) as where, ξ jl and η jl are standard Gaussian random variables, which are independent of each other and of k jl , have zero mean, and unit variance.Note that the indices j = 1, 2, . . ., n and l = 1, 2, . . ., n 0 .In this strategy of stratified sampling of the wavenumber bins, we choose the sampling bins such that they are equally spaced in the log-space, that is with respect to ln k.Thus we have, for a given index j, kj+1 = q kj , where q is some ratio parameter of the geometric distribution which defines the wavenumber intervals ∆ j = ( kj , kj+1 ].Since we have a minimum wavelength in our model, we have k1 = k * .The bounded domain for k also introduces the constraint, k ≤ k max , where we choose k max = 100, such that kn+1 = k max .For this work, we chose n = 4 and n 0 = 40.We have also verified that pushing n and n 0 to higher values of 5 and 50 respectively does not change the results.It is important to note here that this strategy of stratified sampling is efficient in way that allows us to work with a comparatively low value of total number of modes, O(100), to give us accurate and physical scale independent results, as compared to just sampling the modes trivially, thus, gaining some computational advantage.The density profile extends from r start = 0.105×10 4 km to r end = 11.35 × 10 4 km.This region is sampled in spatial steps of δx = 10 km.The turbulence is introduced to this density profile, starting at the reverse shock position which is at, r rs = 0.39 × 10 4 km and it is terminated at the position of the forward shock, r fs = 5.6 × 10 4 km.The region of turbulence has the same spatial discretization scale as the original density profile (= 10 km).As stated earlier, the primary parameters that we vary and examine in this work are the amplitude associated with the turbulence C * and the power spectral index α.In Fig. 2, we show the effective matter potential for a fixed α = −5/3 and various values of C * .The red thick line is the case when turbulence is absent, that is, C * = 0.This is essentially the same curve as shown in Fig. 1 (solid blue curve) rescaled to convert the matter density to effective potential using natural units ℏ and c.As the amplitude of turbulence increases, the fluctuations grow as expected.The effects of this variation in the amplitude on the survival probabilities for neutrinos are discussed in the following section.It is important to note here that we want to focus on the effects associated solely with the non-adiabaticity at the forward and reverse shocks, hence we artificially smoothed the contact discontinuity so as to avoid any non-adiabaticity effects arising from it.

B. Results
In this section, we discuss the results obtained by solving for neutrino flavor propagation through the matter profile as shown in Fig. 2. The two primary scenarios we are interested in are when turbulence is absent and when turbulence has a high amplitude.Thus, we compare the results for C * = 0 (which we label as without turbulence) and C * = 0.3 (which we label as with turbulence) respectively to implement the above scenarios.The choice of C * = 0.3 is representative of an extreme case where complete flavor depolarization occurs and the oscillation signatures are washed out.
The results are shown in Fig. 3 for the electron neutrino survival probability (P ee (r)), the survival probability for the first instantaneous matter eigenstate (P 11 (r)), the effective mixing angle in matter (θ M ) and the adiabaticity parameter (γ(r)) as ν e propagates through the effective matter profile from r start to r end .In Fig. 3a we show the survival probability (P ee (r)) for ν e as a function of distance for the cases of with (solid orange) and without (dotted purple) turbulence.We note that the survival probability initially starts from unity, corresponding to a pure ν e state.The position of reverse shock introduces a drop in the survival probability.Between the forward and the reverse shocks, the behaviour of P ee (r) is clearly distinct for the cases with and without turbulence.The propagation of the neutrinos in the presence of turbulence is extremely non-adiabatic due to the highly fluctuating background matter density.A better understanding can be achieved by studying P 11 (r) (in Fig. 3b), which is expected to track the changes in the density.The matter eigenstate survival probability shows a jump due to non-adiabatic evolution corresponding to the density jumps.While the case without turbulence has jumps only at the positions of the reverse and forward shocks, the case with turbulence has numerous jumps.
A similar trend can be seen in the effective mixing angle in the matter for the two cases, as shown in Fig. 3c.The case without turbulence exhibits a maximal change in θ M at the locations of the forward and reverse shocks, whereas, in the presence of turbulence, θ M has large jumps throughout the density profile.To demonstrate the non-adiabaticity associated with the neutrino flavor evolution, Fig. 3d shows the associated adiabaticity parameter, γ(r).A highly non-adiabatic evolution is characterized by γ(r) << 1.We see that for the case with turbulence between the reverse and the forward shocks, the adiabaticity parameter is much smaller than 1, implying a highly non-adiabatic evolution.The turbulence is only seeded between the forward and reverse shocks, so beyond the forward shock, the evolution once again becomes adiabatic in the absence of any further shocks or turbulence.This can be seen from Figs. 3b  and 3d, where in the former case, the transition probability does not show any jumps beyond the forward shock, whereas in the latter case, γ(r) >> 1 in the same region, implying adiabatic evolution.

C. On the double-dip associated with shocks
The presence of shocks characterized by rapid changes in the density profile can have imprints on the neutrino oscillation probabilities.The MSW resonance associated with neutrino oscillations is given by, ∆m (low density) and H-resonance (high density) respectively.The propagation of the shockwave through the densities corresponding to H-resonance regions can break adiabaticity leading to distinctive features in the neutrino energy spectra.
As discussed in Sec.I, numerical simulations hint towards the presence of both forward and reverse shocks in the matter density profile, resulting in a double-dip feature in the neutrino spectra [77].In Fig. 4 we show the average survival probability for ν e as a function of the neutrino energy (E ν ).Note the presence of the double-dip feature (solid red curve) in the neutrino spectra corresponding to the density profile (in Fig. 2) in the absence of turbulence (C * = 0).The first peak corresponds to E ν ∼ 18 MeV and is due to the forward shock.Following the first peak, the adiabatic transitions corresponding to the forward and reverse shocks partially cancel out to give the valley.Finally, a second peak appears for E ν ≳ 60 MeV which is due to the neutrinos crossing one of the forward or reverse shocks.
We now aim to study the effects of the presence of turbulence on the double-dip feature in the neutrino spectra.The results are shown in Fig. 4 for a range of turbulence amplitudes given a fixed value of α = −5/3.Interestingly, we find that the double-dip feature vanishes as the amplitude of turbulence increases leading to highly non-adiabatic flavor evolution.For C * > 0.01, the first peak disappears.As the value of C * is increased, a gradual wash-out of the feature is seen.In fact, for the maximum value of C * used in this work, i.e., C * = 0.3, the survival probability averages to 1/2, signifying complete flavor depolarization.Thus, the presence of turbulence can completely wash away the double-dip feature in the neutrino energy spectra associated with forward and reverse shocks.This can be used to probe the presence of hydrodynamic turbulence in core-collapse supernovae.

D. Revisiting previous analytical results for a simplified model of turbulence
In [83], the effects of a simple model for turbulence on neutrino propagation in supernovae were studied.In particular, analytical estimates corresponding to neutrino flavor depolarization in the presence of a Kolmogorov-type turbulence were derived.The authors argued that evidence from numerical simulations hints towards the development of Kolmogorov cascades due to various instabilities.In such a turbulence scenario, the Fourier transform of the velocity correlator has the following form, where α = −5/3.The matter potential, A(x), was assumed to be linear and turbulent random noise was introduced as where F is the amplitude or normalization factor for the noise, k is the mode number, β = α/2, and ϕ(k) is a random phase which is a random value between 0 and 2π for every mode k.The above equation can be compared to Eq. ( 2) used in this work.We reproduce the numerical results from [83] using the same parameters to check our numerical code.We show the result in Fig. 5a.The spatial coordinate goes from x = −100 to x = 100, ∆ = 50, and θ = 0.1.We choose ∆x = 0.1.The agreement between our results and the same obtained by the previous work is good.We find that for low turbulence amplitudes, the effect on survival probability is less pronounced, which is expected and matches with our own analysis discussed above.When the turbulence amplitude is increased from F ∼ 0.2 to higher values, a power law develops.The power law rise seen in the figure can be explained using an analytical estimate derived in [83], where n ′ 0 is related to the spatial derivative of the matter density.We note that for values of F ≳ 0.05, a complete flavor depolarization happens, implying the flavor constitution becomes equal for ν e and ν x , in the two-flavor oscillation scenario considered here.
We extend the analysis in [83] by studying the variation of the average electron neutrino survival probability (⟨P ee ⟩) in the β − F plane.The result, in the form of a density plot for the average probability1 , is shown in Fig. 5b.A rainbow color scheme is used to illustrate ⟨P ee ⟩, where the violet end of the spectrum implies ⟨P ee ⟩ ∼ 0 and the red end signifies ⟨P ee ⟩ ∼ 1.The area where complete flavor depolarization occurs is shown by a greenish-yellow shade corresponding to ⟨P ee ⟩ ∼ 0.5.For β ∼ 0, we find that ⟨P ee ⟩ ∼ 0.5 for any given value of the turbulence amplitude F .As β decreases to −1, higher amplitudes of turbulence lead to flavor depolarization.This generalizes the results of [83] and confirms that for large values of F , flavor depolarization holds for other values of β as well.
There are distinct differences between the analysis in [83] and our current work.We assume a more realistic matter density profile based on simulations in contrast to the linear profile assumed for simplicity in the former work.While the simple profile helps in getting some analytical estimates, the effects of a realistic profile need to be addressed numerically and lead to more realistic predictions for the upcoming detectors.Furthermore, the seeding of the turbulence in this work uses a variant of the Randomization method which is different from the simple random noise turbulence used in the previous work.

IV. EVENTS RATES AND STATISTICAL ANALYSIS
In previous sections, we discussed the effects of turbulence on the neutrino signatures and survival probabilities.This section focuses on the prospects and the capabilities of the upcoming neutrino detectors to constrain the parameters associated with turbulence, in particular, the amplitude (C * ) and the spectral index (α).The detectors we consider for this work are the Deep Underground Neutrino Experiment (DUNE) and Hyper-Kamiokande (HK).While the former is largely sensitive to a flux of ν e , the latter is better suited for detecting the νe flux.
The supernova neutrino spectrum at the source (the neutrinosphere) is assumed to be a pinched Fermi-Dirac spectrum [100], where N src ν is the number of neutrinos at the source for any given flavor with energy E ν , L ν is the neutrino luminosity for that particular flavor, ⟨E ν ⟩ is the mean neutrino energy, α ν is a pinching parameter associated with the given flavor.Considering the cooling phase, we choose L νe = L νe = 5 × 10 51 erg, L νx = L νx = 10 52 erg, ⟨E νe ⟩ = 12 MeV, ⟨E νe ⟩ = 15 MeV, ⟨E νx ⟩ = ⟨E νx ⟩ = 18 MeV, and α νe = α νe = α νx = α νx = 3.The neutrinos emitted from the neutrinosphere then propagate through the matter profile encountering the forward and reverse shocks, along with the turbulence which leads to neutrino oscillations.Let P ee be the survival probability of ν e while propagation.In this case, the number of ν e neutrinos at the end of the matter profile is given by where N SN νe is the number of ν e after the initial neutrino spectra from the neutrinosphere have propagated through the matter profile.Similarly one can obtain the spectra for νe .Once the neutrinos leave the supernova they propagate through the interstellar medium to reach Earth.This can be effectively treated as propagation in vacuum2 .The vacuum survival probability for ν e is given by where θ vac is the vacuum mixing angle and ∆m 2 vac is the vacuum mass-squared difference.We assume θ vac = 0.9967 rads and ∆m 2 vac = 2.6 × 10 −3 eV 2 .Thus the neutrino spectra at the Earth, neglecting Earth-matter effects, is given by A similar analysis can be repeated for obtaining the νe spectra on earth which will be relevant for HK.We do not consider earth matter effects for this work.The total number of neutrinos of a given species α (N να ) detected in a neutrino detector from a source at a distance R is given by where N tar is the number of targets for a given channel in the detector fiducial volume, σ να is the neutrino interaction cross-section with the targets for the relevant channel, ∆t is the duration for which we assume the density profile to remain approximately fixed and similar to what is considered for the matter profile, and W is the Gaussian energy resolution function which depends on the intrinsic properties of the detector.The experimentally reconstructed energy in the detector is given by E r and the actual neutrino energy is given by E t ≡ E ν .We fix the distance to the supernova at R = 10 kpc for our current work.The energy resolution function is defined as [101] W The sensitivity of the neutrino detectors for different values of C * and α is estimated using a simple χ 2 estimator where N obs is the number of neutrino events observed in the neutrino detector, N exp is the expected number of neutrino events in the detector.We fix exp to the case when C * = 0.1, and α = 1.65, corresponding to the Kolmogorov turbulence, unless otherwise specified.We consider ∆t = 1 s.

A. DUNE
The DUNE experiment [102] consists of a near and a far detector separated by 1300 km.The former is located at Fermilab, Illinois while the latter is located about 1.5 km below the ground at the Stanford Underground Research Facility (SURF), South Dakota.The far detector is a liquid argon time projection chamber (LArTPC) with 70 kt (fiducial mass of 40 kt) of liquid argon.We will focus on the far detector in this work.
The relevant interaction channel for supernova neutrinos is the charged current interaction of electron neutrinos with argon nuclei, ν e + 40 Ar → 40 K * + e − .We use N tar = 6 × 10 32 .The interaction cross-section σ νe associated with the ν e − Ar interaction cross section is taken from [103].The energy resolution for DUNE is taken to be σ E /MeV = 0.11 E r /MeV + 0.02E r /MeV.
The result for the total number of electron neutrino events N νe for DUNE is shown in the left panel of Fig. 6 as a density plot in the C * − α plane.The color map uses a rainbow template to denote the lowest number of events in violet and the highest in red.For a typical galactic scale supernova (R ∼ 10 kpc) we have ∼ 103 events in DUNE in the ν e channel.There are two interesting aspects that we note right away -(a) as α varies from 1.0 to 2.0, the number of events N ν does not change significantly, which implies that it would be difficult to put constraints on the value of α using the DUNE detector; (b) as C * is varied from 0.001 to 1.0 the number of events in the DUNE detector has a significant variation and goes from ∼ 1350 (red) for lower values of C * ∼ 0.001 to ∼ 900 (violet) for large values of C * ∼ 1.0.Thus, DUNE can be used to place constraints on the value of C * for a typical galactic supernova in the near future.We discuss the two points below using a χ 2 plot as discussed below.The significant variation of the number of events with C * is confirmed in Fig. 7a, which shows the variation of χ 2 versus C * for α = 1.65.The null hypothesis is considered to be C * = 0.1.We note that for low values of C * < 0.05, the χ 2 starts growing and it increases as the amplitude of turbulence decreases 3 .However, for larger values of C * > 0.05, constraints become weak.To get a better understanding of this behaviour, we compare the event spectra for different values of C * in Fig. 7c.The null hypothesis (C * = 0.1) is shown with a dotted dark blue line.The cases corresponding to C * = 0.005, 0.02, and 0.5 are shown as solid purple, orange, and red lines.On comparing the case of C * = 0.005 with the expected case we see a significant difference in the number of events in the energy bins above 17 MeV leading to the high value of χ 2 .As C * increases the number of events in the energy bins becomes similar to the expected case of C * = 0.1 explaining the trend of χ 2 as seen in Fig. 7a.
The variation over α does not have a specific trend for any value of C * and the χ 2 values obtained are extremely low and we do not show them here.This is evident from Fig. 8a where we show the event spectra for α = 1.05 (solid orange), 1.65 (dashed dark blue), and 2.0 (solid purple) for C * = 0.1.The expected scenario for this case is again chosen to be C * = 0.1 and α = 1.65 shown as the dashed dark blue line.We note that the spectra for all the three cases shown look very similar to the expected case in dotted dark blue, explaining the low values of χ 2 .This highlights the ineffectiveness of DUNE to put constraints on the value of α.

B. Hyper-Kamiokande
Hyper-Kamiokande (HK) [104] is an upcoming 260 kton (fiducial 187 kton) liquid scintillator detector, being built at the Kamioka Observatory in Japan.It is an upgrade to the Super-Kamiokande detector but with around 10 times the fiducial volume.It will consist of two tanks of ultrapure water acting as the medium for neutrino detection.The primary detection channel in HK is the inverse-beta decay (IBD), νe + p → e + + n.We use the IBD cross-section from [105].Similar to Super-Kamiokande, the energy resolution for HK is given by σ E /MeV = 0.6 E r /MeV [103].The number of targets for HK is assumed to be 1.2 × 10 35 .
The density plot for the number of νe events for HK in the C * − α plane is shown in the right panel of Fig. 6.The details of the colors used are discussed in Sec.IV A. We have ∼ 7 × 10 4 νe events for a typical galactic scale supernova (R = 10 kpc).This number is much higher than DUNE owing to the larger number of targets.This already hints at the possibility that HK would be a better detector than DUNE at constraining the parameter range for turbulence.However, similar to DUNE, the variation in N νe is insignificant across the different values of α for a given value of C * .
Using HK, relatively strong constraints can be placed on the amplitude of turbulence as seen by the variation in N νe with C * for a fixed value of α in the right panel of Fig. 6.This is also evident from Fig. 7b where we have very large values of χ 2 for C * ≲ 0.05, keeping α fixed at 1.65.To confirm our results, we show the event spectra for α = 1.65 and different values of C * in Fig. 7d.Once again, the expected case is assumed to be C * = 0.1 and α = 1.65 and is shown in dotted dark blue.We note that similar to the case of DUNE, for E ν ≤ 20 MeV, the spectra are similar among all the cases, however for E ν > 20 MeV, there is a significant difference in the spectra between C * = 0.1 and C * = 0.005 (solid purple line) and C * = 0.02 (solid orange line) leading to the comparatively high values of χ 2 .However, for C * = 0.5 (solid red line) the event spectra look similar to the expected case giving a low value of χ 2 as seen in Fig. 7b.
Finally in Fig. 8b, we show the event spectra for C * = 0.1 for different values of α.The spectra look similar to the expected case of C * = 0.1, indicating that HK, while sensitive, would not be able to put very strong constraints on α.

V. CONCLUSIONS AND IMPLICATIONS
Neutrino propagation through a turbulent medium leads to high degrees of non-adiabaticity.Hydrodynamical simulations of core-collapse supernova hint towards the presence of various instabilities in the region behind the shock, leading to the matter profile becoming turbulent.Such a turbulent medium can leave its signature on the neutrino spectra, thereby allowing future neutrino detectors to probe the strength of turbulence inside a core-collapse supernova.In this work, we performed a detailed study of the effects of turbulence on neutrino propagation and outlined the prospects for upcoming neutrino detectors like DUNE and HK to constrain the parameter space associated with turbulence.
We considered a usual supernova matter profile snapshot, characteristic of the cooling phase, and used the variant C formalism of the Randomization method to seed in the turbulence.This involves a stratified sampling of the modes to have a minimal effect associated with the physical scales.The two main parameters associated with turbulence are the amplitude (C * ) and the spectral index associated with the power spectrum (α).We numerically evolved the neutrino oscillation probabilities using the slab approximation, which is appropriate for a highly oscillatory profile as is the case in the presence of turbulence.It is important to note that for getting reasonable results, the discretization scale, which in this case is given by the width of the slabs, should be smaller than the length scales associated with the variation of the matter profile.To illustrate the difference between neutrino propagation in the presence and absence of turbulence, we contrasted various quantities like the electron neutrino survival probability, the effective mixing angle in matter, and the adiabaticity parameter.We showed that these quantities change significantly depending on whether the matter background has turbulence or not.
Depending on the location of the forward and reverse shocks, the neutrino spectra from a supernova can be very different due to the non-adiabaticity associated with the shock fronts.In particular, the H-resonance associated with ∆m 2 atm can break adiabaticity at the forward and reverse shocks.This can lead to a double-dip structure in the neutrino survival probability, as was discussed in detail in [77].We revisited this issue and studied if the double-dip feature still survives in the presence of turbulence between the shockfronts.We found that as the profile becomes more turbulent, flavor depolarization kicks in.In fact, in our study, if the strength of turbulence exceeds 30%, the doubledip structure is washed out, leading to complete flavor depolarization.This implies that the double-dip signature which can provide more information about the positions of the forward and reverse shocks may fail to do so in the presence of strong turbulence in between the shocks.
We considered the prospects of upcoming neutrino detectors like DUNE and HK to detect neutrinos from a galactic core-collapse supernova and constrain the parameters associated with turbulence in Sec.IV.Given the larger fiducial volume, HK is better in constraining the turbulence parameters than DUNE.Most importantly, we note that there is a significant variation of the number of events as the amplitude increases from sub-percent to O(1) for a fixed value of the spectral index.As a result, both DUNE and HK would be able to place reasonable constraints on the amplitude of the turbulence.However, the same cannot be said for the spectral index as the sensitivity is not too large.We conclude by discussing potential improvements in our results.With regard to the turbulence, we assumed a simplified model of one-dimensional to illustrate our results.A more realistic picture can be obtained using threedimensional computations where the turbulence samples k x , k y , and k z from given power spectral distributions.However, this would be computationally intensive.Furthermore, we also assumed that the matter density profile along with the turbulence is constant for the duration of neutrino propagation we consider, which is set to 1 second.This is a reasonable approximation since the density profiles from numerical simulations do not have significant changes over such durations in the cooling phase.
The results presented in this paper rely on an effective two-flavor framework.This can be further improved to include three-flavor effects.Furthermore, we also did not consider the effect of collective flavor oscillations on our results.Particularly, if the neutrinos are susceptible to fast flavor instabilities occurring near the core, then flavor depolarization might already have happened by the time the neutrinos reach the shock-front [41,62].In that case, the turbulence signature becomes degenerate and cannot be disentangled from the depolarization due to fast conversions.
Understanding the complex nature of turbulence along with its very existence still remains a challenging problem.This work highlights the possibilities of constraining turbulence in the era of large-scale neutrino detectors like DUNE and HK.The statistics from a future galactic supernova at the upcoming neutrino detectors should be good enough to gain invaluable insights into turbulence.In fact, for a very nearby supernova (within 1 kpc) [106,107] the complete spectrum can help place strong constraints on both the amplitude and the spectral index associated with turbulence, and allow us to investigate the double-dip in the neutrino spectra associated with shocks.A planned 100-kt neutrino detector THEIA [108] can also help place strong constraints owing to its large fiducial volume.We conclude that with the next galactic core-collapse supernova in the era of DUNE and HK, we can significantly improve our understanding of the evolution of shocks inside the SN, and the underlying hydrodynamic instabilities.This can help to answer the question about whether the turbulence inside a core-collapse supernova follows a Kolmogorov spectrum.

FIG. 1 .
FIG.1.The density profile snapshot at time t = 5 s post collapse[77].The various regions of interest in the shock wave profile are labelled.

FIG. 5 .
FIG. 5. (a) The electron neutrino survival probability as a function of the noise amplitude F for β = −5/6, to match with the results of [83] (see Fig. 1 there); (b) The electron neutrino survival probability as a function of β ranging from −1 to 0 for various values of F .Each data point in the numerics was repeated for 66 times and the average value was used.

FIG. 6 .
FIG.6.Density plot to show the total number of νe events for DUNE (left) and νe events for HK (right) on the C * − α plane.The source is assumed to be at a distance of R = 10 kpc.See text for details on the other parameters used.
FIG. 7. (a) Plot of χ 2 estimator for different values of C * for DUNE assuming α = 1.65 (corresponding to Kolmogorov turbulence); (b) same as (a) but for HK; (c) The expected number of νe events per energy bin for DUNE assuming α = 1.65 and C * = 0.005 (solid purple), 0.02 (solid orange), and 0.5 (solid red).The Nexp for the χ 2 estimator is considered to be the case when C * = 0.1 and α = 1.65, which is shown as a dotted dark blue line; (d) same as (c) but for HK where νe events are considered.

FIG. 8 .
FIG.8.The expected number of νe and νe events per energy bin corresponding to DUNE (left) and HK (right) for a fixed value of C * = 0.1 and different values of α = 1.05 (solid orange), 1.65 (dotted dark blue), and 2.0 (solid purple) respectively.The source is assumed to be at a distance of R = 10 kpc.