Decoherence in neutrino propagation through matter, and bounds from IceCube/DeepCore

We revisit neutrino oscillations in matter considering the open quantum system framework which allows to introduce possible decoherence effects generated by New Physics in a phenomenological manner. We assume that the decoherence parameters $\gamma_{ij}$ may depend on the neutrino energy, as $\gamma_{ij}=\gamma_{ij}^{0}(E/\text{GeV})^n$ $(n = 0,\pm1,\pm2) $. The case of non-uniform matter is studied in detail, both within the adiabatic approximation and in the more general non-adiabatic case. In particular, we develop a consistent formalism to study the non-adiabatic case dividing the matter profile into an arbitrary number of layers of constant densities. This formalism is then applied to explore the sensitivity of IceCube and DeepCore to this type of effects. Our study is the first atmospheric neutrino analysis where a consistent treatment of the matter effects in the three-neutrino case is performed in presence of decoherence. We show that matter effects are indeed extremely relevant in this context. We find that IceCube is able to considerably improve over current bounds in the solar sector ($\gamma_{21}$) and in the atmospheric sector ($\gamma_{31}$ and $\gamma_{32}$) for $n=0,1,2$ and, in particular, by several orders of magnitude (between 3 and 9) for the $n=1,2$ cases. For $n=0$ we find $\gamma_{32},\gamma_{31}<4.0\cdot10^{-24} (1.3\cdot10^{-24})$ GeV and $\gamma_{21}<1.3\cdot10^{-24} (4.1\cdot10^{-24})$ GeV, for normal (inverted) mass ordering.


Introduction
The accurate measurement of the mixing angle θ 13 by reactor neutrino experiments [1], with a small uncertainty comparable to that for θ 12 , has initiated a precision era for neu-CERN-TH-2018-041; IFT-UAM/CSIC-18-022; FERMILAB-PUB-18-067-T a pcoloma@fnal.gov b jacobo.lopez.pavon@cern.ch c ivanj.m@csic.es d nunokawa@puc-rio.br trino physics. In the standard three-family framework, the main remaining issues are the possible observation of leptonic CP violation, the determination of the ordering of neutrino masses and probing the Dirac or Majorana nature of neutrinos. Some hints currently exist in the latest data collected by NOvA and T2K which seem to point to maximal CP violation in the neutrino sector, but the statistical significance is still low [2,3]. Likewise, a global fit to neutrino oscillation data seems to show a mild preference for a normal mass ordering (see for instance [4,5]), which needs to be confirmed as more data become available.
At the same time, and in view of the precision of present and near future neutrino facilities, it is of key importance to verify if neutrinos have unexpected properties caused by New Physics (NP) beyond the standard three-family framework. In this work we study one of the possible windows to NP, the so-called quantum decoherence in neutrino oscillations, and update the existing bounds by analyzing IceCube and DeepCore data on atmospheric neutrinos. In particular, we are interested in a kind of decoherence effects in neutrino oscillations studied, for example, in [6][7][8][9][10][11][12][13] and, more recently, in [14][15][16][17][18][19]. These decoherence effects differ from the standard decoherence caused by the separation of wave packets (see e.g. [20,21]) and might arise, instead, from quantum gravity effects [22][23][24]. Throughout this work, for brevity, we will refer to such non-standard decoherence simply as decoherence.
The authors of ref. [7] derived some of the strongest available constraints on neutrino decoherence in neutrino oscillations up to date, using atmospheric neutrino data from the Super-Kamiokande (SK) experiment [25][26][27][28]. Moreover, they considered the general case in which the decoherence parameters could depend on the neutrino energy via a power law, γ = γ 0 (E/GeV) n , where n = 0, −1, 2. Nevertheless, these limits were obtained within a simplified two-family framework and without taking into account the matter effects in arXiv:1803.04438v2 [hep-ph] 21 Jul 2018 the neutrino propagation. Furthermore, only a reduced subset of SK data (taken, in fact, almost 20 years ago now) was analyzed [25][26][27][28].
In this work, we show that performing a three-flavour analysis which includes the matter effects is essential in order to correctly interpret such constraints. In particular, it is not obvious to which γ i j parameter the SK bounds derived in two families [7] would actually apply. We will show that it strongly depends on the neutrino mass ordering and on whether the sensitivity is dominated by the neutrino or antineutrino channels: for neutrinos the decoherence effects at high energies are mainly driven by γ 21 (γ 31 ) for normal (inverted) ordering, while in the antineutrino channel they are essentially controlled by γ 32 (γ 21 ) for normal (inverted) ordering. Concerning the solar sector, the authors of ref. [15] obtained strong constraints on γ 21 from an analysis of Kam-LAND data [29,30], for n = 0, ±1. 1 Finally, the authors of ref. [32] derived several bounds on the atmospheric decoherence parameters γ 32 and γ 31 from an analysis of MINOS data [33][34][35].
Non-standard decoherence has been invoked several times in the literature in order to decrease the tension in the parameter space among different sets of neutrino oscillation data. For example, in refs. [13,14] a solution to the LSND anomaly based on quantum decoherence, compatible with global neutrino oscillation data, was proposed. More recently, in [17] it was shown that the ∼ 2σ tension between T2K and NOvA on the measurement of the atmospheric mixing angle θ 23 could be alleviated through the inclusion of decoherence effects in the atmospheric neutrino sector, namely, γ 23 = (2.3 ± 1.1) · 10 −23 GeV. Such value of γ 23 would be close to the SK bound from ref. [7], γ < 3.5 · 10 −23 GeV (90% CL), but still allowed. This topic has recently brought the attention of a part of the community. In fact, several analyses of decoherence effects on present and future long-baseline neutrino oscillation experiments have been recently performed (albeit at the probability level only), see e.g. refs. [16,18,19]. In this work we will show that the reference value for γ 23 considered in [17] is indeed already excluded by IceCube data (we note however that, according to the latest results reported by NOvA, the significance of the tension has been reduced to less than 1σ [3]). Moreover, we find that Ice-Cube and DeepCore data are able to improve significantly over most of the constraints in past literature, both for solar and atmospheric decoherence parameters, in some cases by several orders of magnitude.
The paper is structured as follows. In sec. 2 we present the formalism and discuss the effects of decoherence on the oscillation probabilities. We first review the case of constant matter density profile, and then proceed to discuss the case of non-uniform matter. In particular we show that, within the adiabatic approximation, no significant bounds on the decoherence parameters can be extracted from solar neutrino data when the neutrino energy is assumed to be conserved. We then proceed to develop a formalism which permits a consistent treatment of the decoherence effects on neutrino propagation in non-uniform matter when the adiabaticity condition is not fulfilled, as is the case of atmospheric neutrino experiments. In sec. 3 we apply this formalism to the computation of the relevant oscillation probabilities in the atmospheric neutrino case, discussing the main features arising in presence of decoherence. Section 4 summarizes the main features of the IceCube and DeepCore experiments, the data sets considered in our analysis, and the details of our numerical simulations. Our results are then presented and discussed in sec. 5. Finally, we summarize and draw our conclusions in sec. 6. Appendix A and Appendix B discuss technical details regarding some of the approximations used in our numerical calculations.

Quantum decoherence: Density matrix formalism
The evolution of the density matrix ρ in the neutrino system can be described as where H is the Hamiltonian of the neutrino system and the second term D [ρ] parameterizes the decoherence effects. In vacuum, the diagonal elements of the Hamiltonian are given by h i = m 2 i /(2E), where m i (i = 1, 2, 3) are the masses of the three neutrinos and E is the neutrino energy. Here ρ is defined in the flavour basis, with matrix elements ρ αβ . Throughout this work, we will use Greek indices for flavor (α, β = e, µ, τ), and Latin indices for mass eigenstates (i, j = 1, 2, 3).
A notable simplification of eq. (1) can be achieved via the following set of assumptions. First, assuming complete positivity, the decoherence term D [ρ] can be written in the so-called Lindblad form [36,37] where D m is a general complex matrix. Second, avoiding unitarity violation, which is equivalent to imposing the condition d Tr[ρ]/dt = 0, requires D m to be Hermitian. Moreover, D m = D † m implies that the entropy S = Tr[ρ ln ρ] increases with time. Finally, a key assumption is the average energy conservation of the neutrino system, which is satisfied when In presence of matter effects, the Hamiltonian is diagonalized by the unitary mixing matrix 2Ũ (throughout this paper, in our notation the presence of a tilde denotes that a quantity is affected by matter effects). Therefore, after imposing the energy conservation condition given by eq. (3), we get This implies that the average energy is conserved along the whole neutrino propagation (through vacuum and matter). This assumption is indeed crucial for our analysis. It is expected to be fulfilled in vacuum and in very good approximation in matter. While we assume that the quantum decoherence itself does not cause the violation of energy conservation, due to the standard neutrino interaction with matter, for large neutrino energies the energy is not exactly conserved in presence of matter due to a small energy transfer to the background fermions. Therefore, in this case, eq. (4) does not hold exactly. This issue has been analyzed in detail in [16,19,31], where it was shown that in a more general framework in which energy conservation is not assumed, two types of effects in the neutrino oscillation probabilities can essentially be distinguished: pure decoherence effects which suppress the oscillating terms, and the so called relaxation effects which affect non-oscillating terms. In [10,31] it was shown that, for atmospheric neutrino oscillations, the relaxation effects which arise when the energy is not conserved are proportional to cos 2 2θ 23 . This suppresses relaxation effects with respect to pure decoherence effects by at least two orders of magnitude, since cos 2 2θ 23 is currently constrained by experimental data at the level cos 2 2θ 23 < 0.034 at 95% CL [4,5]. We will thus focus on the analysis of pure decoherence effects in atmospheric neutrino oscillations assuming that the neutrino energy is conserved, and therefore eq. (4) satisfied, along the whole propagation. From a model-independent point of view, the d j m are free parameters that could a priori depend on the matter effects. However, the most common assumption in the literature is to assume that the d j m are independent of the matter density even in presence of matter effects. 3 . In order to be consistent with most previous studies and to compare the bounds obtained in our analysis with the constraints derived in previous publications, we will also assume that this is the case. Notice that this assumption does not imply that the matter effects are not relevant when neutrino propagation is affected by decoherence: it just implies that the d j m are assumed to be constant during neutrino propagation through the Earth.

Neutrino propagation in uniform matter
Performing the following change of basis eq. (1) can be rewritten as If the matter profile is constant along the neutrino path, the system of equations becomes diagonal inρ i j where we have defined Therefore, the solution of eq. (1) for constant matter is simply given by with whereρ i j (0) is determined by the initial conditions of the system. For instance, if the source produces only neutrinos of flavor α, the initial conditions are given bỹ As a result, the oscillation probabilities in presence of decoherence (for a constant matter profile) read whereρ denotes the density operator. Finally, after some manipulation the above equation can be rewritten in the more familiar form with , where ∆m 2 i j ≡m 2 i −m 2 j are the effective mass squared differences of neutrinos in matter and we have used the approximation L ≈ t, L being the distance traveled by the neutrino as it propagates. Note that the power-law dependence on the neutrino energy given by eq. (14) breaks Lorentz invariance, except for the case with n = −1 which gives similar effects to the neutrino decay (see e.g. [38]). However, the effect encoded in γ i j only suppresses the oscillatory terms in the oscillation probability while a neutrino decay would also affect the non-oscillatory terms. Therefore, in the framework considered in this work the total sum of the probabilities adds up to 1, while this is not the case for neutrino decay.
From eqs. (13) and (14), one would expect to have a sizable effect in neutrino oscillations for γ i j L ∼ 1. This condition gives an estimate of the values of γ i j for which an effect may be experimentally observable: Nevertheless, we would like to remark that fulfilling this condition is not enough to have sensitivity to decoherence effects, as we will discuss in the next subsection. Even though in our simulations we will numerically compute the exact oscillation probabilities, in order to understand qualitatively the impact of decoherence on the oscillation pattern it is useful to derive approximate analytical expressions. In this work, we will be focusing on the study of atmospheric neutrino oscillations, for which the oscillation channel P µ µ is most relevant. Recently, in [39,40] approximated but very accurate analytical expressions for the standard oscillation probabilities in presence of constant matter density were derived. For the ν µ → ν µ oscillation channel including decoherence effects, using the same parametrization as in ref. [40], we find: where A i j ≡ A i j (θ 23 ,θ 12 ,θ 13 , δ ) = 2|U µi (θ 23 ,θ 12 ,θ 13 , δ )| 2 |U µ j (θ 23 ,θ 12 ,θ 13 , δ )| 2 , (17) and the effective mass splittings and mixing angles in matter can be expressed as [40]: cos 2θ 13 = cos 2θ 13 − a/∆ m 2 ee (cos 2θ 13 − a/∆ m 2 ee ) 2 + sin 2 2θ 13 , cos 2θ 12 = cos 2θ 12 − a /∆ m 2
2.2 Neutrino propagation in non-uniform matter: adiabatic regime Equation (13) applies for constant density profiles (which is a very good approximation in the case of long-baseline neutrino oscillation experiments such as T2K or NOvA), but if the matter density is not constant the analysis becomes more complicated. Nevertheless, when the adiabaticity condition dŨ/dt 1 is fulfilled, as in the solar neutrino case, the solution of the evolution equations given by eqs. (9) and (10) is still a good approximation. In such a case, the oscillation probability is given by where ν e f f i denotes the effective mass eigenstates at time t. In the case of solar neutrinos, the initial flux of ν e is produced in the solar core and the initial conditions are given by: whereŨ 0 denotes the effective mixing matrix at the production point. On the other hand, since the evolution is adiabatic, when the neutrinos come out from the Sun we have |ν e f f i = |ν i and thus Finally, for solar neutrinos observed at the Earth we obtain, after averaging over the oscillating phase: which coincides with the standard result. In other words, the decoherence effects encoded in γ i j can not be bounded by solar neutrino oscillation experiments. This is due to the standard loss of coherence in the propagation from the Sun to the Earth, which strongly suppresses the oscillating terms.
In a similar fashion high-energy astrophysical neutrinos at IceCube are not sensitive to decoherence either, since these neutrinos are produced in distant astrophysical sources and thus the oscillations will have averaged out by the time they reach the detector.

Neutrino propagation in non-uniform matter: layers of constant density
In the atmospheric neutrino case, the matter profile cannot be considered constant since the neutrinos propagate through the Earth crust, mantle and core, which have different densities. The adiabaticity condition is not fulfilled either. In this case, eq. (6) should be solved including the non-adiabatic terms, which give non-diagonal contributions. Even though this can be done numerically, we will show that dividing the matter profile into layers of constant density considerably simplifies the analysis and reduces the computational complexity of the problem. In particular, this is crucial in the case of atmospheric neutrino oscillation experiments, for which numerical studies are already computationally demanding even in the standard three-family scenario. Dividing the matter profile into layers of different constant densities has proved to be a very good approximation in the standard three-family scenario and, therefore, we expect the same level of accuracy in presence of decoherence. Since the matter is constant in each layer, the evolution equations can be solved for each layer M as in sec. 2.1: where ∆t M ≡ t M −t M,0 , and t M,0 and t M denote the initial and final times for the propagation along layer M, respectively. Now the problem of computing the probability is just reduced to performing properly the matching among the evolution on the different layers. Let us first consider the simplest case of two layers A and B. The oscillation probability when the neutrino exits the second layer (at time t B ) is given by The key point here is that the matching should be done between the solutions of eq. (1) at the frontier between the two layers and in the flavor basis, as After imposing the matching condition, the elements of the density matrix in the second layer at t B,0 can be written in the matter basis as: where we have considered that the initial flux is made of ν α as initial condition for the first layer. After substituting this result into eq. (24) we finally obtain It can be easily checked that, in the limit γ i j → 0, the standard oscillation probability is recovered. In the three-layer case, following the same procedure we find The procedure can be easily generalized to an arbitrary number of layers. Indeed, under the approximation L ≈ t, and defining the probabilities can be written in a more compact way as for two layers for three layers for N layers 3 Atmospheric oscillation probabilities with decoherence Atmospheric neutrino oscillations take place in a regime where matter effects are significant and can even dominate the oscillations. The relevance of matter effect increases with neutrino energy and is very different for neutrinos and antineutrinos, as the sign of the matter potential changes between the two cases. Matter effects also depend strongly on the neutrino mass ordering. In order to understand better the numerical results shown in this paper, it is useful to derive approximate expressions for the oscillations in the ν µ → ν µ andν µ →ν µ channels in the presence of strong matter effects.
From the results obtained in refs. [39,40], for neutrino energies E > ∼ 15 GeV matter effects drive the effective mixing angles in matterθ 12 andθ 13 to either 0 or π/2, depending on the channel (neutrino/antineutrino) and the mass ordering. It is easy to show that, in this regime, the oscillation probabilities in eq. (16) can be approximated as: 23 1 − e −γ 21 L cos∆ 21 (33) for neutrinos, and for antieneutrinos, assuming a normal ordering (NO). For inverted ordering (IO) we get instead for neutrinos, 23 1 − e −γ 21 L cos∆ 21 (36) for antineutrinos. In the limit γ i j = 0, eqs. (33)-(36) reassemble the standard neutrino oscillation probabilities derived under the one-dominant mass-scale approximation [41]. From eqs. (33)- (36) it is easy to see that the approximated oscillation probabilities for an inverted mass ordering can be obtained from the corresponding ones for normal mass ordering, just performing the following transformation: Moreover, note that since the three decoherence parameters and the three mass splittings are related (see eqs. (8) and (14)), these two transformations automatically imply that Equations (33)-(36) illustrate why a proper consideration of the matter effects in the context of three families is of key importance in order to correctly interpret the bounds extracted within a simplified two-flavour approximation (as done in e.g. ref. [7]). According to our analytical results, which will be confirmed numerically below, the constraints obtained from SK in a two-family approximation cannot be simply applied to γ 31 or γ 32 , contrary to the naive expectation. In fact, the interpretation of such limits depends strongly on the ordering of neutrino masses and on whether the sensitivity is dominated by the neutrino or antineutrino channels: for neutrinos the decoherence effects at high energies would be mainly driven by γ 21 (γ 31 ) for normal (inverted) ordering. Conversely, in the antineutrino channel decoherence effects are essentially controlled by γ 32 (γ 21 ) for normal (inverted) ordering. Therefore, we conclude that in order to avoid any misinterpretation of the bounds from atmospheric neutrinos, a three-family approach including matter effects should be considered. Figure 1 shows the numerically obtained ν µ → ν µ (top panels) and ν µ → ν µ (bottom panels) oscillation probabilities for NO (left panels) and IO (right panels), with and without decoherence, as a function of the neutrino energy for a three-layer model (details on the accuracy of our three-layer model and the specific parameters used in our simulations can be found in Appendix A). For the sake of simplicity, in this section we focus on the case n = 0, where the γ i j do not depend on the neutrino energy (the results for different values of n show a similar qualitative behavior). The standard oscillation parameters have been fixed to the best-fit values given in [4,5]. Figure 1 clearly shows how the decoherence tends to damp the oscillatory behavior, in qualitative agreement with eq. (16) and the corresponding approximate expressions given by eqs. (33)- (36). Note that eq. (16) has been obtained under several approximations and, in particular, considering only one layer with constant matter density. However, we should stress that in our simulations the computation of the probability has been done numerically, considering a three-layer matter profile (see Appendix A for details).  The ν µ → ν µ (top panels) and ν µ → ν µ (bottom panels) oscillation probabilities with (n = 0) and without decoherence effects as a function of the neutrino energy. The probabilities have been computed for normal (left panels) and inverted (right panels) neutrino mass ordering, using a three-layer model for the Earth matter density profile, and correspond to the case in which the neutrinos cross the center of the Earth. In this figure, in the cases where decoherence effects are included we have set the parameters listed in the legend to the same constant value γ = 2.3 · 10 −23 GeV.
In Appendix B, we will show that the bounds derived in these limits correspond to the most conservative bounds that can be extracted in the general case. As a reference value for the decoherence parameters, in this section we have considered γ = 2.3 · 10 −23 GeV, for each of the three limiting cases listed above.
The results in fig. 1 show that, for neutrinos with a NO (top left panel), the impact of decoherence is essentially controlled by γ 21 , in good agreement with eq. (33): no signifi-cant effects are seen in the atmospheric limit (A), while a similar impact is obtained in the solar limits I (B) and II (C). Conversely, for IO (top right panel) the effects are dominated by γ 31 instead: no effect is observed for the solar limit II (C), while in scenarios (A) and (B) the effect is very similar. This can be qualitatively understood from the approximate probability derived in eq.  ing equivalence for the results obtained in the three limiting cases listed above: This is also confirmed at the numerical level, as it can be easily seen by comparing the different lines shown in the left (NO) and right (IO) panels in fig. 1 for the three limiting cases.
It is also remarkable that, for both normal and inverted mass orderings, even when the standard oscillations turn off (at very high energies), there is still a large effect on the probability due to decoherence effects, that could potentially be tested with neutrino telescopes like IceCube. In particular, for E > ∼ 200 GeV one can approximate cos∆ i j ≈ 1, ∀i, j.
Therefore, in the standard case (with γ i j = 0) the last three terms in eq. (16) approximately vanish, leading to P µ µ ≈ 1. However, in presence of decoherence those terms will not vanish completely, as e −γ i j L cos∆ i j = 1. This leads to a depletion of P µ µ , which is no longer equal to 1 in this case. The size of the effect will of course depend on the baseline of the experiment. Since at high energies the oscillation probability does no longer depend on the neutrino energy, at oscillation experiments with a fixed baseline the effect may be hindered by the presence of any systematic errors affecting the normalization of the signal event rates. However, at atmospheric experiments this effect can be disentangled from a simple normalization error by comparing the event rates at different nadir angles.
The dependence of the neutrino probabilities with the zenith angle θ z is illustrated in fig. 2, assuming a normal mass ordering and fixing the standard oscillation parameters to the best-fit values given in [4,5]. The results are shown as a neutrino oscillogram (see for instance [42]), which represents the oscillation probability in the P µ µ channel as a function of neutrino energy and zenith angle θ z (which can be related to the distance traveled by the neutrino). Figure 2 shows the oscillation probability P µ µ in the three limiting cases described above, comparing it to the results in the standard scenario (γ i j = 0). As expected, the effects depend on the direction of the incoming neutrino and they are more relevant in the region −1 cos θ z −0.4, this is, for very long baselines. This was to be expected, since the decoherence effects are driven by e −γ i j L . In addition, the dependence of the oscillation probability with the zenith angle at very high energies (E > 100 GeV) is clearly visible in the bottom panels of fig. 2. As we will show in sec. 5, this will lead to an impressive sensitivity for the IceCube setup. Finally, note that the results for inverted ordering show similar features to those in fig. 2, once the mapping in eq. (40) is applied, and are therefore not shown here.

IceCube/DeepCore simulation details and data set
The IceCube neutrino telescope, located at the South Pole, is composed of 5160 DOMs (Digital Optical Module) deployed between 1450 m and 2450 m below the ice surface along 86 vertical strings [43]. In the inner core of the detector, a subset of these DOMs were deployed deeper than 1750 m and closer to each other than in the rest of IceCube. This subset of strings is called DeepCore. Due to the shorter distance between its DOMs, the neutrino energy threshold in DeepCore (∼ 5 GeV) is lower than in IceCube (∼ 100 GeV). This allows DeepCore to observe neutrino events in the energy region where atmospheric oscillations take place, see fig. 1, whereas IceCube only observes high-energy atmospheric neutrino events.
As outlined in sec. 2, for high energy astrophysical neutrinos the effect of non-standard decoherence in the probability would be completely erased by the time they reach the detector. Therefore, in this work we will focus on the observation of atmospheric neutrino events at both IceCube and DeepCore, in the energy range ∼ 10 GeV to ∼ 1 PeV. In particular, we have used the three-year DeepCore data on atmospheric neutrinos with energies between ∼ 10 GeV and ∼ 1 TeV, taken between May 2011 and April 2014 [44], and the one-year IceCube data taken between 2011-2012 [45][46][47], corresponding to neutrinos with energies between 200 GeV and 1 PeV.
At IceCube and DeepCore, events are divided according to their topology into "tracks" and "cascades" [48]. Tracks are produced by the Cherenkov radiation of muons propagating in the ice. In atmospheric neutrino experiments, muons are typically produced by two main mechanisms: (1) via charged-current (CC) interactions of ν µ with nuclei in the detector, and (2) as decay products of mesons (typically pions and kaons) originated when cosmic rays hit the atmosphere. Conversely, cascades are created in CC interactions of ν e or ν τ 4 : in this case, the rapid energy loss of electrons as they move through the ice is the origin of an electromagnetic shower. At IceCube/DeepCore, cascades are also observed as the product of hadronic showers generated in neutral-current (NC) interactions for neutrinos of all flavors. Our analysis considers only track-like events observed at both IceCube and DeepCore although, as we will see, some small contamination from cascade events can be expected (especially at low energies).

IceCube simulation details
For IceCube, the observed event rates are provided in a grid of 10 × 21 bins [46], using 10 bins for the reconstructed energy (logarithmically spaced, ranging from 400 GeV to 20 TeV), and 21 bins for the reconstructed neutrino direction (linearly spaced, between cos θ rec z = −1.02 and cos θ rec z = 0.24). The muon energy is reconstructed with an energy resolution σ log 10 (E µ /GeV) ∼ 0.5 [45], while the zenith angle resolution is in the range σ cos θ z ∈ [0.005, 0.015], depending on the scattering muon angle.
The number of events in each bin is computed as: where E, θ z denote the true values of energy and zenith angle, while E rec , θ rec z refer to their reconstructed quantities. Here, φ atm µ,± is the atmospheric flux for muon neutrinos (+) and anti-neutrinos (-), P ± µ µ (E, θ z ) is the neutrino/antineutrino oscillation probability given by eq. (28), and A eff i,±,µ (E, θ z ) is the effective area encoding the detector response in neutrino energy and direction (which relates true and reconstructed variables), the interaction cross section and a normalization constant, and has been integrated over the whole data taking period. In our IceCube simulations, we have used the unpropagated atmospheric flux (HondaGaisser) provided by the collaboration [45,50], and for the effective area we have used the nominal detector response from refs. [45,50]. On the other hand, the exponential factor takes into account the absorption of the neutrino flux by the Earth, which increases 4 Technically, a CC ν τ event could be distinguished from a ν e CC event, e.g., by the observation of two separates cascades connected by a track from the τ propagation [49]. However, for atmospheric neutrino energies the distance between the cascades cannot be resolved by the DOMs at IceCube/DeepCore, leaving in the detector a signal similar to a single cascade.  [46] are represented by the black dots, and the error bars indicate the statistical uncertainties for each bin.
with the neutrino energy. Here, X(θ z ) is the column density along the neutrino path and σ ± (E) is the total inclusive cross section for ν µ orν µ . Note that in eq. (41) no contamination from cascade events is considered since the mis-identification rate is expected to be negligible at these energies [51]. Similarly, the number of atmospheric muons that pass the selection cuts can also be neglected, given the extremely good angular resolution at these energies [45]. Figure 3 shows the expected number of events for Ice-Cube from our numerical simulations including decoherence, for γ 21 = γ 31 = 2.3 · 10 −23 GeV (solid blue lines) and γ 21 = γ 31 = 10 −22 GeV (dot-dashed green lines), as a function of cos θ rec z , for events in different reconstructed energy ranges. For simplicity, we have considered the n = 0 case (that is, γ i j independent of the neutrino energy). The expected result without decoherence is also shown for comparison (dashed red lines), while the observed data [46] are shown by the black dots.
For the analysis of the IceCube data we have performed a Poissonian log-likelihood analysis doing a simultaneous fit on the following parameters: ∆ m 2 32 , θ 23 and γ i j . The rest of the oscillation parameters have been kept fixed to their current best-fit values from ref. [4,5]. The most relevant systematic errors used in the fit are summarized in tab. 1, and have been taken from refs. [45,47,50]. For each systematic uncertainty a pull term is added to the χ 2 following the values listed in the table, except in the cases indicated as "Free" (when the corresponding nuisance parameter is allowed to float freely in the fit).

DeepCore simulation details
In the case of DeepCore, the observed event rates [44] are provided in a grid of 8 × 8 bins, using 8 bins for the reconstructed neutrino energy and 8 bins for the reconstructed neutrino direction. The energy resolution σ E/GeV is in the range of 30%-20% while the zenith angle resolution improves with the energy, from σ θ z = 12 • at E ν = 10 GeV to σ θ z = 5 • at E ν = 40 GeV [44].
In each bin, the number of events is computed as   [44] are represented by the black dots, and the error bars indicate the statistical uncertainties for each bin.
Unlike for IceCube, at DeepCore muon tracks can be produced from ν µ → ν µ and ν e → ν µ events 5 . Moreover, the track-like event distributions at DeepCore will also receive partial contributions from cascades which are mis-identified as tracks: hence the sum over β = e, µ, τ in eq. (42). Therefore, here φ atm α,± stands for the atmospheric flux for neutrinos/antineutrinos of flavor α (where we have used the fluxes from ref. [52]), and P ± αβ refers to the neutrino/antineutrino oscillation probability in the channel ν α → ν β for neutrinos (+) (orν α →ν β , for antineutrinos (-)). The rejection efficiencies for the contamination are included in the detector response function A eff i,±,β , which now depends on the flavor β of the interacting neutrino. Finally, an estimate of the atmospheric muons that overcome the selection criteria (taken from refs. [44,50]) is also added for each bin in reconstructed variables, N i,µ . Figure 4 shows the expected number of events for Deep-Core obtained from our numerical simulations including decoherence, for γ 21 = γ 31 = 2.3 · 10 −23 GeV (solid blue lines) and γ 21 = γ 31 = 10 −22 GeV (dot-dashed green lines), as a function of cos θ rec z , for events in different reconstructed energy ranges. For simplicity, we have considered the n = 0 case (that is, γ i j independent of the neutrino energy). The expected result without decoherence is also shown for comparison (dashed red lines), while the observed data [44] are shown by the black dots. 5 The flux from ν τ can be considered negligible at these energies. Table 2: Systematic errors used in our analysis of DeepCore data, taken from refs. [44,48].

Source of uncertainty Value
Flux -normalization Free Flux -energy dependence as (E/E 0 ) η ∆ η = 0.05 Flux -(ν e +ν e )/(ν µ +ν µ ) ratio 20% Background -normalization Free DOM efficiency 10% Optical properties of the ice 1% In this work a Gaussian maximum likelihood is used to analyze the DeepCore data, performing a simultaneous fit on the following parameters: ∆ m 2 32 , θ 23 and γ i j . The rest of the oscillation parameters have been kept fixed to their current best-fit values from refs. [4,5]. The systematics used in the fit are those associated with the flux, the detector response and the atmospheric muons given in ref. [44] and are summarized in tab. 2. For each systematic uncertainty a pull term is added to the χ 2 following the values listed in the table, except in the cases indicated as "Free" (when the corresponding nuisance parameter is allowed to float freely in the fit). We have checked that our analysis reproduces the confidence regions in the ∆ m 2 32 − θ 23 plane obtained by the DeepCore collaboration in ref. [44] to a very good level of accuracy.
Finally, it should be noted that our fit does not include the latest atmospheric data recently published by the Deep-Core collaboration [53]. The new analysis uses a different data set (from April 2012 to May 2015) and a new implementation of systematic errors, which lead to smaller confidence regions in the ∆ m 2 32 − θ 23 plane. However, the detector response parameters and systematic errors used in the latest release have not been published yet. In view of the better results obtained for the standard three-family oscillation scenario, a similar improvement is to be expected if the analysis performed in this work were to be repeated using the latest DeepCore data.

Results
Following the procedure described in sec. 4 we have obtained the χ 2 for every point in the parameter space. Marginalizing over the relevant mixing and mass parameters, namely, ∆ m 2 32 and θ 23 , the sensitivity of the data to γ i j parameters is determined by evaluating the ∆ χ 2 , with ∆ χ 2 ≡ χ 2 − χ 2 min , where χ 2 min is the value at the global minimum. In this section we will only show the results obtained for NO, since we have checked that extremely similar results are obtained for IO after applying the mapping given in eq. (40). Nevertheless, in sec. 6 we will also provide the 95% confidence level (CL) bounds obtained in our numerical analysis for the IO case. The bounds obtained are in very good agreement with the mapping given in eq. (40). Figure 5 shows the obtained ∆ χ 2 as a function of γ 0 for the three limiting cases defined in sec. 3: (A) atmospheric limit, γ 0 = γ 0 32 = γ 0 31 (red curve); (B) solar limit I, γ 0 = γ 0 21 = γ 0 31 (green curve); and (C) solar limit II, γ 0 = γ 0 21 = γ 0 32 (blue curve). In all cases, the solid (dashed) lines correspond to the results obtained from our analysis of the IceCube (DeepCore) data, and each panel shows the results obtained assuming a different energy dependence for the decoherence parameters, see eq. (14): n = 0 (top panel), n = 1 (middle panel) and n = 2 (bottom panel). The shaded regions are disfavored by previous analysis of SK [7] (90% CL) and KamLAND [15] data (95% CL). As explained in sec. 3, the KamLAND constraints derived in [15] apply to γ 0 12 (solar limits) while it is not clear to which γ i j the bounds from SK obtained in [7] would apply, since this depends on the true neutrino mass ordering (which is yet unknown).
Note that the size of the atmosphere has been neglected in our calculations (see Appendix A for details). This is a good approximation for small values of the decoherence parameters, but it starts to fail if the decoherence effects are large enough to affect neutrinos with cos θ z > 0. Therefore, in the case of IceCube we have shown our results only in the region where this approximation holds. In the case of Deep-Core, due to the smaller energies considered, our approxi-mation has no impact on the results even for large values of the decoherence parameters. Therefore, the approximation has only an impact on the IceCube results in a region of the decoherence parameter space which is already ruled out either by DeepCore or other experiments. Figure 5 shows that for both DeepCore and IceCube the best sensitivity is achieved for the solar limits (B) and (C) while the weakest limit is obtained in the atmospheric limit (A). In particular, the strongest bound is obtained for (C). This is in agreement with the behaviour of the oscillation probability in presence of strong matter effects, discussed in sec. 3. On one hand, as shown in sec. 3, for NO the decoherence effects are mainly driven by γ 21 in the neutrino channel and γ 32 in the antineutrino channel. On the other hand, the number of antineutrino events is going to be suppressed with respect to the neutrino case, due to the smaller cross section and flux. Hence, the best sensitivity is expected for case (C), where γ 0 = γ 0 21 = γ 0 32 , since both neutrinos and antineutrinos are sensitive to decoherence effects. Conversely, in case (B), where γ 0 = γ 0 21 = γ 0 31 , only neutrinos are sensitive to decoherence effects, and therefore some sensitivity is lost with respect to the results for case (C). Finally, in case (A), with γ 0 = γ 0 32 = γ 0 31 , the bounds come mainly from the impact of decoherence on the antineutrino event rates and, since these are much smaller than in the neutrino case, the obtained bounds are much weaker when compared to the results obtained in case (B). Figure 5 shows a flat asymptotic feature of the Deep-Core's ∆ χ 2 for large values of γ 0 , where the sensitivity becomes independent of γ 0 . Conversely, for IceCube there is a decrease in sensitivity for values of γ above a certain range: for example, for n = 0 the best sensitivity is achieved for γ 0 ∼ O(10 −22 ) GeV while it decreases for higher values. This behaviour can be understood as follows. For the neutrino energies observed at IceCube (above 100 GeV) the oscillation phases do not develop and the probabilities do not depend on the energy (cos∆ i j ≈ 1 in eq. (16)). Therefore, at IceCube the sensitivity to the decoherence effects comes from the observation of a non-standard behaviour of the number of events with the zenith angle alone. Naively, eq. (15) gives the values of L and γ that yield a large effect. Considering n = 0, for example, where there is a one-to-one relation between the two, we get that for γ 0 ∼ 10 −22 GeV the effect will be maximal for distances of the order L ∼ O(10 3 ) km. This is the typical distance traveled by atmospheric neutrinos crossing the Earth and therefore the sensitivity of IceCube is maximized in this range. Conversely, for larger (smaller) values of γ 0 , only neutrinos coming from the most horizontal (vertical) directions are affected, leading to a reduced impact on the χ 2 .
From the comparison between the different panels in fig. 5 we can see that the limits change considerably with the value of n, which parametrizes the energy dependence Values of the ∆ χ 2 as a function of the decoherence parameter for the Atmospheric limit (red), Solar limit I (green) and Solar limit II (blue) defined in sec. 3. The results obtained from our analysis of IceCube (DeepCore) data are denoted by the solid (dashed) lines. The three panels have been obtained for NO, assuming a different dependence on the neutrino energy: n = 0 (top panel), n = 1 (middle panel) and n = 2 (bottom panel). The shaded regions are disfavored by previous analysis of SK [7] and KamLAND [15] data, see text for details. The horizontal black line indicates the value of the ∆ χ 2 corresponding to 95% CL for 1 degree of freedom. of the decoherence parameters (see eq. (14)). In particular, we observe in fig. 5 that the sensitivity improves with n and that, as it is increased, the results for IceCube improve much faster (compared to DeepCore) due to the much higher neutrino energies considered in this case. The behaviour of the sensitivities with the value of n is better appreciated in fig. 6, where we show the bounds obtained at 95% CL (for 1 degree of freedom) as a (discrete) function of the power-law index n, for n = −2, −1, 0, 1 and 2. The DeepCore bounds are rep-resented by solid circles while the IceCube constraints are given by the solid triangles. The results seem to follow the linear relation where E 0 2.5 TeV (30 GeV) for IceCube (DeepCore). This can be understood as follows. Decoherence effects enter the oscillation probabilities only through the factor γL = γ 0 (E/GeV) n L, for any value of n. Naively, we expect that : 95% CL bounds on the decoherence parameters γ 0 , for NO, as a (discrete) function of the power-law index n for the Atmospheric limit (red), Solar limit I (green) and Solar limit II (blue). The solid circles (triangles) correspond to the DeepCore (IceCube) analysis.
the sensitivity limit is obtained for γL ∼ O(1) (although the precise value will eventually depend on the neutrino mass ordering, on the particular γ i j which drives the sensitivity, and on the data set considered). Taking the logarithm of γ 0 (E/GeV) n L = constant, we reproduce eq. (43). At first approximation, the value of E 0 in eq. (43) can be estimated as the average energy of the IceCube and DeepCore event distributions, E , as where dN/dE is the event number distribution. This leads to E 4 TeV (40 GeV) for IceCube (DeepCore), which are in the right ballpark although somewhat different from the values of E 0 giving the best fit to the data shown in fig. 6. Nevertheless, we find these to be in reasonable agreement, given our naive estimation of E 0 as the mean energy for each experiment.

Summary and Conclusions
In this work, we have derived strong limits on non-standard neutrino decoherence parameters in both the solar and atmospheric sectors from the analysis of IceCube and DeepCore atmospheric neutrino data. Our analysis includes matter effects in a consistent manner within a three-family oscillation framework, unlike most past literature on this topic. In sec. 2 we have developed a general formalism, dividing the matter profile into layers of constant density, which permits to study decoherence effects in neutrino oscillations affected by matter effects in a non-adiabatic regime. Our analysis shows that the matter effects are extremely relevant for atmospheric neutrino oscillations and their importance in or-der to correctly interpret the two-family limits obtained previously in the literature, as outlined in sec. 3.
We have found that the sensitivity to decoherence effects depends strongly on the neutrino mass ordering and on whether the sensitivity is dominated by the neutrino or antineutrino event rates. For neutrinos, the decoherence effects at high energies are mainly driven by γ 21 (γ 31 ) for normal (inverted) ordering, while in the antineutrino case they are essentially controlled by γ 32 (γ 21 ) for normal (inverted) ordering. This means that, considering a three-family framework including matter effects is essential when decoherence effects in atmospheric neutrino oscillations are studied. Our results are summarized in tab. 3, together with the most relevant bounds present in the literature. Table 3 provides the 95% CL bounds extracted from our analysis of DeepCore and IceCube atmospheric neutrino data, for both normal and inverted ordering, and for the three limiting cases considered in this work: (A) atmospheric limit (γ 21 = 0), (B) solar limit I (γ 32 = 0) and (C) solar limit II (γ 31 = 0). In Appendix B we show that the bounds derived in these limits correspond to the most conservative results that can be extracted in the general case.
In this work, we considered a general dependence of the decoherence parameters with the energy, as γ i j = γ 0 i j (E/GeV) n with n = ±2, 0, ±2. Our results improve over previous bounds for most of the cases studied, with the exception of the n = −1 case. For n = −1, KamLAND gives the dominant bound on γ 21 while MINOS gives the strongest constraints on γ 31 and γ 32 6 (indeed, both KamLAND and MINOS are also expected to give the strongest bound for n = −2, although to the best of our knowledge no analysis has been performed for this case yet). We have found that DeepCore considerably improves the present bounds in the solar sector (γ 21 ) for n = 0, 1, 2 and gives a constraint in the atmospheric sector comparable to the SK limit, although a factor 2 weaker, in the n = 0 case. Our results show that, for n = 0 (which is the case most commonly considered in the literature), Ice-Cube improves the bound on γ 31 and γ 32 in (more than) one order of magnitude with respect to the SK constraint, obtained in a simplified two-family approximation, and by more than one order (almost two orders) of magnitude for NO (IO) with respect to the KamLAND constraint on γ 21 . In particular, we find that the reference value for γ 23 considered in ref. [17] to explain the small tension previously reported between NOvA and SK data is indeed already excluded by IceCube data. Regarding the cases with n = 1, 2, we find that the sensitivity of IceCube is particularly strong.   Table 3: DeepCore/IceCube bounds on γ 0 i j in GeV (γ i j = γ 0 i j (E/GeV) n ), at the 95% CL (1 degree of freedom), for both normal and inverted ordering as indicated. Previous constraints are also provided for comparison, and the dominant limit in each case is highlighted in bold face (notice that we considered the most conservative bound from the two solar limits). in 4 (5) orders of magnitude with respect to the SK limit for NO (IO). MSCA-RISE-2015.This manuscript has been authored by Fermi Research Alliance, LLC under Contract No. DE-AC02-07CH11359 with the U.S. Department of Energy, Office of Science, Office of High Energy Physics. The publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes.

Appendix A: Computation of oscillation probabilities in three layers
The simulation of atmospheric neutrino experiments is computationally demanding in the standard three-family scenario, and even more if decoherence effects are included in the analysis. Therefore, due to the cost of implementing a large number of layers for the PREM profile density, in this work we consider a simplified three-layer model for the Earth matter density profile assuming a core and Earth radii of 3321 km and 6371 km, respectively. The values of the matter densities of the inner layer (core) and the outer layer (mantle) are taken to be around ρ = 12 g/cm 3 and 4.6 g/cm 3 , respectively. However, their values are slightly adjusted depending on the distance traveled by the neutrinos to match as close as possible the profile of the PREM model [54]. Note that, in our simulations, we have not considered the atmosphere as an additional layer. This is a good approximation for neutrinos going upwards in the detector (cos θ z < 0), but is not a valid approximation in the region cos θ z > 0. This has only an impact on the IceCube results for extremely large values of the decoherence parameters, which are already ruled out either by other experiments or by DeepCore. In fig. 7 we compare the results obtained for the oscillation probability for our modified three-layer approximation (left panel) against the exact numerical results using the full Preliminary Reference Earth Model (PREM) profile [54] (right panel), which divides the Earth into eleven layers where the matter density in each layer is given by a polynomial function of the distance traveled. In this figure, the results are shown for the standard three-family scenario with no decoherence, in order to illustrate the accuracy of our three-layer approximation. The results are shown as a neutrino oscillogram, which represents the oscillation probability in the P µ µ channel in terms of energy and the zenith angle θ z of the incoming neutrino. In this figure, a normal mass ordering was assumed, together with the following input values for the oscillation parameters [4,5]: ∆ m 2 21 = 7.4 · 10 −5 eV 2 , ∆ m 2 31 = 2.515 · 10 −3 eV 2 , θ 12 = 33.62 • , θ 13 = 8.54 • , sin 2 θ 23 = 0.51, and δ = 234 • .
As can be seen from the comparison between the two panels, some small differences take place but only in a restricted range of values of energy and zenith angle. Therefore, we conclude that the agreement between the probabilities obtained using the exact PREM model (right) and our approximate three-layer model (left) is sufficiently good for the purposes of this work. We have also checked that, using our simplified three-layer model applied to the standard case without decoherence, we are able to reproduce up to a very good approximation the DeepCore oscillation fit for the atmospheric parameters θ 23 and ∆ m 2 32 [44].

Appendix B: Five-dimensional analysis
The γ i j are not completely independent parameters, see eq. (8).
In order to simplify the analysis, in this work we have focused on three different representative cases: (A) Atmospheric limit, γ 21 = 0 (γ 32 = γ 31 ); (B) Solar limit I, γ 32 = 0 (γ 21 = γ 31 ); and (C) Solar limit II, γ 31 = 0 (γ 21 = γ 32 ). Considering these one-γ i j -dominated cases is expected to be a very good approximation in view of eqs. (33)- (36). Nevertheless, in this appendix we will show that the results obtained in these simplified scenarios also apply to the more general case in which the three γ i j are different from zero. Let us assume that just one D m matrix contributes to the decoherence term of the evolution equations given by eq. (2). In such a case, one of the γ i j parameters is a function of the other two γ i j . Without loss of generality, if we choose γ 21 and γ 31 as our free parameters, γ 32 is then given by In order to understand how general are the results presented in sec. 5, we have performed a five-dimensional analysis varying γ 21 , γ 31 , θ 23 and ∆ m 2 32 in the fit, and imposing the constraint given by the equation above. In fig. 8 we show the ∆ χ 2 obtained from the five-dimensional DeepCore analysis as a function of γ 21 (dashed green curve) and γ 31 (dashed red curve), marginalizing over the rest of the free parameters, for the n = 0 case (the same conclusions apply to the other cases studied in this work). For the sake of comparison, the ∆ χ 2 associated to the atmospheric (solid red  Fig. 8: ∆ χ 2 obtained from the five-dimensional DeepCore analysis as a function of γ 21 (dashed green curve) and γ 31 (dashed red curve), marginalizing over the rest of the free parameters, for n = 0 and NO. The ∆ χ 2 for the Atmospheric (solid red curve), Solar I (solid green curve) and Solar II (solid blue curve) limits is also shown. curve), solar I (solid green curve) and solar II (solid blue curve) limits is also included in the same figure. NO was assumed but the results can be easily extrapolated to the IO case using the mapping given in eq. (40). Figure 8 shows that the five-dimensional ∆ χ 2 distribution projected into γ 31 coincides with the Atmospheric limit one, while when it is projected into γ 21 resembles the most conservative of the two solar limits. This is due to the marginalization over the parameters which are not shown. For instance, in the case of γ 21 the marginalization selects, between the two solar limits, the most conservative result. We conclude therefore that our analysis distinguishing the three limits (A), (B) and (C), provides the most conservative bounds that can be applied to the general case in which the three γ i j are different from zero.