Neutrino oscillation bounds on quantum decoherence

We consider quantum-decoherence effects in neutrino oscillation data. Working in the open quantum system framework we adopt a phenomenological approach that allows to parameterize the energy dependence of the decoherence effects. We consider several phenomenological models. We analyze data from the reactor experiments RENO, Daya Bay and KamLAND and from the accelerator experiments NOvA, MINOS/MINOS+ and T2K. We obtain updated constraints on the decoherence parameters quantifying the strength of damping effects, which can be as low as $\Gamma_{ij} \lesssim 8 \times 10^{-27}$ GeV at 90% confidence level in some cases. We also present sensitivities for the future facilities DUNE and JUNO.


I. INTRODUCTION
The outstanding experimental progress of the past years has allowed to establish a clear picture of neutrino oscillations [1,2] and to measure most neutrino oscillation parameters with unprecedented accuracy [3][4][5].While few open questions remain to be answered, such as the determination of possible leptonic CP violation, the neutrino mass ordering and the octant of the atmospheric angle, the next generation of neutrino experiments will likely allow to close up the neutrino physics picture and to determine all oscillation parameters with good precision.
The established oscillation paradigm, and the consequent interpretation of current and forthcoming neutrino oscillation data, relies on the three-neutrino mass and mixing scheme.According to this scenario, neutrinos propagate as a superposition of mass and flavor eigenstates.During propagation, each component mass eigenstate evolves with a different frequency thus resulting in the phenomenon of flavor conversion.Since neutrinos interact only very weakly with other matter particles and since they are stable, this quantum superposition effect is maintained over macroscopic distances.In general, the study of neutrino oscillations considers the neutrino to be isolated from its environment, and the oscillation effects to be coherent.However, neutrino eigenstates may loose their quantum superposition, for example if the neutrino system interacts with a stochastic environment, thus resulting in a loss of coherence and in the damping of neutrino oscillation probabilities -a phenomenon known as neutrino decoherence.
We will work in the framework of open quantum systems [24], commonly used in quantum decoherence studies [14,.We will assume that the oscillation phenomenon can be changed by the interaction of the neutrino subsystem with the environment.The fluctuating nature of space-time in a quantum theory of gravity is a commonly cited potential source of a stochastic background that might produce neutrino decoherence effects [48][49][50][51].In this work, however, we will adopt a phenomenological approach and remain agnostic about its origin.By analyzing several neutrino oscillation experiments, we will set stringent limits on possible decoherence effects.
Experimental searches for neutrino quantum decoherence effects have been performed with a vast range of facilities, and assuming different neutrino sources.In general, the loss of coherence does not prevent neutrino flavor conversion, but modifies the oscillation phenomenon through distance-and energy-dependent damping terms, as the neutrinos propagate away from the source.This suggests that long-baseline experiments are better suited for probing neutrino quantum decoherence.Regarding neutrino sources, astrophysical neutrinos are therefore excellent candidates.Solar neutrinos have been used to place the strongest limits on energy-independent quantum decoherence, as shown in [31,52], in combination with data from the Kamioka Liquid Scintillator Antineutrino Detector (KamLAND).It should be noted, though, that depending on the model of decoherence, all sensitivity to measure decoherence parameters might be lost if neutrino oscillations become averaged, as in the case of solar neutrinos (or the diffuse high energy extragalactic neutrino flux).We will comment on this important point later on.Earlier works made use of atmospheric neutrinos using Super-Kamiokande data to obtain bounds on the quantum decoherence parameters [27].Currently, the most stringent limits on decoherence parameters with positive energy dependence (Γ ∝ E n , n > 0) have been obtained using atmospheric neutrinos observed at the IceCube Neutrino Observatory [39].Note that if the origin of the decoherence effects is related to quantum gravity, one would naturally expect n > 0. Long-baseline accelerator and reactor experiments, like the Main Injector Neutrino Oscillation Search (MINOS) and KamLAND have also been proven capable of setting bounds on decoherence parameters [35,53].Further analyses include long-baseline accelerator data from MINOS and Tokai to Kamioka (T2K) [43], as well as from NuMI Off-Axis ν e Appearance (NOvA) and T2K [36].Sensitivity studies have been presented for the Deep Underground Neutrino Experiment (DUNE) in [40,54] and for the Jiangmen Underground Neutrino Observatory (JUNO) [55].
Building upon these earlier works, we work in the open quantum system framework and analyze data from multiple experiments in order to set up-do-date stringent constraints on the decoherence parameters.Keeping a phenomenological approach, we will parameterize the decoherence effects and their energy dependence in a model-independent way, so to allow for both n ≥ 0 and n < 0. Since there are many possible sources of neutrino decoherence, it is desirable to search for these effects in as many experiments and energy regimes as possible.In this work we compare data from a range of reactor and accelerator experiments to maximize the sensitivity across a broad range of parameter space.As we will show, depending on the energy dependence and on the distance travelled by neutrinos from the source to the detector, different facilities will be more sensitive to the damping effects than others.We analyze data from reactor experiments such as the Reactor Experiment for Neutrino Oscillation (RENO) [56] and Daya Bay [57], which allow for high-resolution, high-statistics measurements of the antineutrino flux and consequently of possible damping features.We also consider KamLAND [58,59], which observes neutrinos from a large number of nuclear reactor sites with longer baselines, thus providing complementary information.We further analyze data from accelerator experiments, which are sensitive to higherenergy neutrinos.We include the MINOS [60], NOvA [61,62] and T2K [63] experiments.The analysis of reactor and accelerator data yields sensitivity to decoherence effects across the MeV to GeV regimes.Finally, we provide sensitivity analyses for two future facilities: JUNO [64,65], a reactor experiment sensitive to the solar oscillation parameters and the atmospheric mass splitting, and especially designed for measuring the neutrino mass ordering, and for the DUNE long-baseline neutrino experiment [66][67][68][69], also expected to determine the neutrino mass ordering, among other goals in its extended physics program.
This paper is organized as follows.We introduce in Section II the open quantum system formalism and the phenomenological scenarios that we analyze.In Section III we present technical details of the experiments and of our statistical analyses.We discuss the results of our analyses in Section IV, presenting updated bounds on the decoherence parameters and sensitivities for near-future experiments.We finally draw our conclusions in Section V.

II. THEORETICAL FRAMEWORK
In this section we present the theoretical framework that we use to obtain the neutrino oscillation probabilities in the presence of decoherence effects.We adopt the open quantum system approach, treating propagating neutrinos as a subsystem which interacts very weakly with a larger, unknown environment [24,70].The interaction of the neutrino system with the quantum environment is the source of damping effects in the oscillations.
The evolution in time of neutrinos in the quantum (sub-)system can be described through the Lindblad equation [71,72] where ρ ν (t) is a Hermitian density matrix describing the neutrino states in the mass basis and H is the Hamiltonian of the neutrino subsystem, which encodes both the standard oscillation effects resulting from the unequal neutrino masses, and any conventional matter effects.According to the density matrix formalism, the density matrix of mixed quantum states is ρ ν = j P j |ψ j ⟩⟨ψ j |, for a system of j states of probability P j .The dissipative term D[ρ ν (t)] encodes the decoherence effects, such that when D[ρ ν (t)] → 0 the standard neutrino oscillations without damping effects are recovered.It is also a Hermitian matrix which can be expressed as [26,27,45,73] where N refers to the dimension of the SU (N ) Hilbert space defining the neutrino subsystem (N = 3 when assuming three neutrino flavors) and O j are dissipative operators described as N × N complex matrices which characterize the coupling of the neutrino subsystem to its environment (j max = 8 for three neutrino flavors), for example in quantum gravity space-time [28,[74][75][76].The operators O j act only on neutrino states in such a way that j O † j O j = I and the trace of ρ ν (t) is unchanged.Let us notice that, while being Hermitian, the matrix D[ρ ν (t)] can be nonunitary due to the interactions of the neutrino subsystem with the environment including possible neutrino loss.The general expression of D[ρ ν (t)] in Eq. (2) would in principle allow for a model-independent analysis of neutrino decoherence effects in oscillation data.However, given the large number of free parameters that it contains, such an analysis is realistically not viable from a practical point of view.To simplify the theoretical framework, focusing on a N -neutrino system, we can expand the operators in Eq (2) in terms of the Gell-Mann matrices from the SU (N ) group [26,41,42,45,73] In the previous expression, c k are the coefficients of the expansion, and the index k runs from 0 to N 2 − 1.In the specific case of three neutrino flavors: k = 0 . . .8, λ 0 is the identity matrix and λ k are the Gell-Mann matrices satisfying [λ a , λ b ] = i c f abc λ c , f abc being the structure constants of SU (3).We can then express the decoherence term as D[ρ ν (t)] = (D kℓ ρ ℓ ν )λ k , where ρ ℓ ν are the coefficients of the neutrino density matrix, expanded in the SU (N ) basis (see Eq. ( 3)), and D kℓ ρ ℓ ν = c k are the elements of a N 2 × N 2 matrix representing the free parameters of the system.Focusing again on a three-neutrino scenario, the matrix D can be parameterized by 45 parameters as follows where all entries are real scalars and the Γ i are positive in order to satisfy the relation Tr(ρ ν (t)) = 1.
As already mentioned, while a general model-independent analysis might be contemplated, it proves convenient to reduce the number of free parameters by imposing general physical conditions on D [26,27,45,54,73,77].These include: • Unitarity of the system.Probability conservation implies D k0 = D 0ℓ = 0 [42], given that f ab0 = 0.
• Complete positivity of the time-evolution ρ ν place conditions on the diagonal elements, thus making them not completely independent.
• Energy conservation of the neutrino subsystem is satisfied through the commutation relation [H, O j ] = 0.This bound includes the decoherence effect in the evolution.
Although in principle energy nonconservation in the neutrino subsystem could be envisaged [26,77], while of course always being satisfied for the global system, in the following we work under the requirement of energy conservation on the neutrino subsystem.These assumptions allow for a remarkable simplification of the dissipator in Eq. ( 4), which eventually assumes the diagonal form [34,35,54] in the case of three neutrino flavors.Comparing with Eq. ( 4), we see that imposing the abovementioned conditions in practice leads to In principle, also the so-called relaxation parameters Γ 3 and Γ 8 could be present in Eq. ( 5).However, since the terrestrial experiments that we consider here are not expected to have competitive sensitivity with bounds from solar neutrinos [52], we set Γ 3 = Γ 8 = 0 throughout the paper.As anticipated, the interaction of the neutrino subsystem with the environment, through the matrix Eq. ( 5) will induce damping effects of the form e −Γ ij L in the oscillation probabilities, L being the experiment baseline and Γ ij quantifying the strength of the damping effects, resulting in a coherence length L coh = 1/Γ ij .While remaining agnostic of the (quantum-gravity) origin of such dissipating terms, we follow a phenomenological approach and introduce a general form for the energy-dependence of the dissipator matrix controlling the decoherence effects.We assume a power-law energy-dependence following previous studies in the literature [27,31,39,41,43,45,79] where E 0 is a pivot energy scale which we set to E 0 = 1 GeV, as previously done in the literature, and n is a power-law index to be tested experimentally.Given the lack of a fundamental theory for quantum gravity, the dependence of Γ ij on the neutrino energy is currently unknown.We will make the hypothesis that it can assume the following integer values: Note that this power-law dependence on the neutrino energy breaks Lorentz invariance, except for the case with n = −1 which imitates the oscillation energy dependence.Our choice for n thus accommodates: a case where decoherence might be induced by lightcone fluctuations [46] (n = −2); a Lorentz-invariant case (n = −1) [27]; the energy-independent case (n = 0); n = +1 which mimics a cross section-like energy-dependence, or as suggested in Ref. [80]; and the case (n = +2) that can arise in quantum-gravity models, in which Γ GeV is expected [81][82][83][84][85].For the sake of completeness, let us also recall that the loss of coherence due to wave-packet separation would correspond to the case n = −4.While we do not cover this scenario here, we refer the reader to Refs.[22,23] for recent analyses.We also note that the n = −1 case is phenomenologically similar to invisible neutrino decay [45].Finally, it has recently been suggested that scenarios with an extreme energy dependence (n ≤ −10) could explain the Gallium anomaly [86] without the need of light sterile neutrinos and hence without entering in conflict with bounds from other short-baseline experiments [87].The choice of D defines the flavor structure of the decoherence effects, and ultimately the flavor ratio of a neutrino ensemble once coherence is completely lost.In the analyses of this paper we consider several scenarios, reflected in activation of different Γ ij simultaneously.The models that we consider are summarized in Table I, and can be grouped into three categories; all three Γ ij activated (A), two Γ ij activated (B, C, D), and only a single Γ ij activated (E, F, G).In each case we test each of the above mentioned energy dependencies.

A. Phenomenological models
Activating differing Γ ij results in differing matrix elements of the [H, ρ ν (t)] term in Eq. (1) being damped.In particular, specific active pairs of Γ ij damp differing oscillations frequencies resulting from the solar (∆m 2  21 ) and atmospheric (∆m 2 3j ) mass splittings.Model D for instance results in the (faster) atmospheric sector oscillations being damped whilst the (slower) solar sector oscillations are not.Such a scenario could therefore represent new physics that depends on the mass splittings (or equivalently oscillation wavelengths), such as lightcone fluctuations [46] or a mass-dependent coupling.Model A represents the case where all oscillation channels are damped together, and can result from perturbations to the neutrino wavefunction phase [45].In general, combining experimental data from various oscillation channels/experiments will be required to optimally test all models.

B. Neutrino oscillation probability and decoherence
The neutrino oscillation probability in presence of the dissipator under consideration here, Eq. ( 5), has been derived many times in the literature, see for example Ref. [54].The result reads where Γ ij (E) is given in Eq. ( 6) and Ũ and ∆ m2 kj are the mixing matrix and mass splittings in matter.In our analyses we calculate these quantities using the expressions provided in Ref. [88], which currently constitute the most precise approximation for neutrino oscillations in constant matter [89].In the following, matter effects will be relevant in the analyses of KamLAND, NOvA and DUNE.Note that the expression of the neutrino oscillation probability is equivalent to the standard expression when Γ ij (E) → 0 (i.e. when no decoherence effects are present).It should be noted that Eq. ( 7) is an approximation, since when including matter effects non-diagonal elements can appear in D after rotation to the matter basis [38].However, this approximation is valid for the analyses performed in this paper, because of several reasons: Firstly, matter effects are never very strong as they would be, for example, in an analysis of atmospheric neutrinos.Secondly, most of the experiments considered here are mainly using disappearance channels, which are much less affected by matter effects than appearance channels.In experiments that use both disappearance and appearance channels, the statistics of the appearance channels are much smaller than those of the disappearance channels.Additionally, the non-diagonal terms would only effect the neutrinooscillation probability at large energies in the tails of the spectra which do not have many events.
We show in Fig. 1 the ν µ (left panel) and νe (right panel) survival probabilities using baselines and energies relevant for the experiments considered in this paper.As can be seen, different choices of n can induce different forms of spectral distortions.We fix Γ ij (E) to values in the ballpark of current bounds, which we will infer in Sec.IV.Experiments using low-energy neutrinos will be more sensitive to scenarios with negative n, while high-energy neutrinos will allow us to set stronger bounds when n > 0, as it is particularly clear from the left plot.In Fig. 2 we keep n = 0 fixed and show the impact of the different models on the oscillation probabilities.From the right panel (representing reactor experiments) we see that the models can be separated into two groups: One group, which has Γ 21 activated and one which has not.If Γ 21 is active we obtain, in addition to possible spectral distortions, also an overall shift of the oscillation probability.Note that a simple shift might not result in an observable effect if one only takes ratios of events at different detectors.Therefore we expect the reactor experiments to have best sensitivity to models which have Γ 21 and at least one other Γ activated, while also using n < 0. Also in the case of long-baseline accelerator experiments (left panel) we find that there exist two groups of models.One group which has Γ 32 active, and one without.As can be seen from the left panel of Fig. 2, the deviation from the standard three-neutrino oscillation probability is particularly large for the former models (A-C-D-G).We can therefore expect accelerator experiments to have best sensitivity to models with (at least) Γ 32 ̸ = 0 and positive n.Of course, this expectation based on oscillation plots is somewhat simplified, since in the real analyses correlations with standard neutrino oscillation parameters and systematic uncertainties of the experiments must be taken into account.Nevertheless, as will be shown later, the strengths of the bounds for the models under consideration will follow this simplified expectation in many cases.

III. EXPERIMENTS INCLUDED IN THE ANALYSIS
In this section we discuss the experimental data that have been considered in this paper and give some details of the analyses performed.We have used data from several reactor and long-baseline accelerator experiments, discussed in Sections III A and III B, respectively.We also detail in Section III C the analysis procedure for the expected sensitivities at the next-generation experiments JUNO and DUNE.It should be noted that the decoherence effects discussed in this paper do not affect solar neutrino measurements.We will therefore always impose a prior on the solar mixing angle, fixing sin 2 θ 12 = 0.318 ± 0.016 [3].
A. Reactor experiments: KamLAND, Daya Bay and RENO KamLAND was a reactor experiment which was placed at the site of the former Kamiokande experiment and was used to measure electron antineutrinos from more than 50 reactors, located at distances ranging from ∼ 100 km to ∼ 1000 km.Due to the long-baselines, matter effects are not negligible at KamLAND, and must be included in the calculation of the neutrino oscillation probability.For the present analysis we include the data presented in Refs.[58,59].The χ 2 function used for KamLAND reads where the index i runs over the energy bins.In the previous expression, N dat,i denote the observed event numbers per energy bin, while N exp,i (⃗ p, ⃗ α) indicate the expected event numbers for a given set of oscillation parameters ⃗ p and systematic uncertainties ⃗ α.The second term is a penalty term on the total number of events, while the third term contains penalty factors for all of the systematic uncertainties α k , with expectation value µ k and standard deviation σ k .The systematic uncertainties that are included in our analysis account for reactor uncertainties (normalizations related to the different reactors and an uncorrelated shape error) and detector uncertainties (detection efficiency and energy scale).The last term contains the above mentioned penalty for sin 2 θ 12 .We use GLoBES [90,91] to compute the event numbers and to perform the statistical analysis.The reactor fluxes are parameterized as in Ref. [92] and the inverse beta-decay cross section is taken from Ref. [93].We have found that the measurement of ∆m 2  21 obtained from the analysis of KamLAND data is robust under the decoherence effects discussed here.We will therefore include another penalty when marginalizing over ∆m 2  21 in all experiments discussed in the following subsections, corresponding to ∆m 2  21 = (7.53± 0.22) × 10 −5 eV 2 .
Apart from KamLAND, we also include data from the reactor experiments Daya Bay and RENO.We analyze the data sets corresponding to 2900 days of running time at RENO [56] and 1958 days of running time at Daya Bay [57].In the statistical analyses, we include uncertainties related to the thermal power for each core, to the detection efficiencies, uncertainties on the fission fractions, a shape uncertainty for each energy bin, and an uncertainty on the energy scale.The χ 2 function for RENO is Here, R F/N i = F i /N i , is the ratio of events at the far (F i ) and near (N i ) detectors in the ith energy bin.In particular, R dat,i are the background-subtracted observed event ratios, while R exp,i (⃗ p, ⃗ α) are the expected event ratios for a given set of oscillation parameters ⃗ p.The uncertainty for each bin is given by σ RENO i .The second term contains penalty factors for all the systematic uncertainties α k , again with expectation value µ k and standard deviation σ k .The last term χ 2 solar contains the solar penalty for sin 2 θ 12 and the KamLAND penalty for ∆m 2  21 .Similarly, for Daya Bay, we define Here, we take the ratios between the far and the first near detector and between the two near detectors.We use the same reactor flux model and inverse beta-decay cross section as in the analysis of KamLAND.

B. Accelerator experiments: MINOS/MINOS+, T2K and NOvA
MINOS was an accelerator-based neutrino oscillation experiment studying muon neutrinos produced at Fermilab at the NuMI beam facility and detected at two detectors located at 1.04 km and 735 km away from the source.First, during the MINOS phase, the neutrino beam peaked at an energy of 3 GeV.Later, during the MINOS+ phase, the energy of the beam was increased, peaking at 7 GeV.Here we consider data corresponding to an exposure of 10.56 × 10 20 proton-ontarget (POT) in MINOS and 5.80 × 10 20 POT in MINOS+, collected in the same detectors [60].We adopted the analysis procedure followed by the experimental collaboration for the search of active-sterile neutrino oscillations in Ref. [60], by modifying the public MINOS/MINOS+ code to account for decoherence effects instead of active-sterile oscillations.The statistical analysis is performed by using the following χ 2 definition: where N dat,i and N exp,i are the observed and the predicted event rates, in the energy bin i.The covariance matrix V cov.accounts for the contributions from several sources of systematic uncertainties [60], and the first sum runs over the two detectors.The last term, χ 2 solar , is the same as in Eqs. ( 9) and (10).
We further consider data collected by T2K [63] and NOvA [62].The T2K collaboration observed events induced by neutrinos and antineutrinos, corresponding to an exposure at Super-Kamiokande of 1.97×10 21 POT in neutrino mode and 1.63×10 21 POT in antineutrino mode.NOvA has instead reached 13.6×10 20 POT in neutrino mode [61] and 12.5×10 20 POT in antineutrino mode.For the energy reconstruction in both experiments we assume Gaussian smearing adding bin-to-bin efficiencies, which are adjusted to reproduce the best-fit spectra reported by the experimental collaborations.Our statistical analysis includes several sources of systematic uncertainties, related to the signal and background predictions.The χ 2 function for the analysis of T2K and NOvA data (and DUNE, as we will see later) is given by where again N dat,i and N exp,i are the data and prediction in bin i and the first sum is taken over the different oscillation channels: ν µ → ν µ , ν µ → ν µ , ν µ → ν e and ν µ → ν e .The last two terms are equivalent to all other cases.

C. Future experiments: JUNO and DUNE
In this subsection we discuss the details for the sensitivity analyses at JUNO and DUNE.When performing these sensitivity analyses we generate a mock data set for both experiments using the neutrino oscillation parameters from Ref. [3] and without decoherence effects.In the statistical analysis we vary the decoherence parameters and marginalize over the standard parameters, equivalently to the other experiments.
JUNO [64] is a next-generation reactor experiment, which aims to precisely measure the solar neutrino oscillation parameters, the atmospheric mass splitting, and the neutrino mass ordering.JUNO will consist of eight reactors [65,94] located at around 52 km distance from the main detector.Our simulation of the experiment follows the descriptions in Refs.[64,65,94,95] assuming 8 years of running time.Our statistical analyses are performed with where the definition of the quantities is equivalent to the other cases, barring that this time N dat,i consists of the mock data set.
DUNE [66][67][68][69] is going to be the successor experiment of NOvA.It will also consist of two detectors, which are going to be exposed to a megawatt-scale neutrino beam produced at Fermilab, composed of (nearly) only muon neutrinos or antineutrinos.The near detector will be placed approximately 600 meters away from the source of the beam, while the far detector, divided into four modules, each using 10 kton of argon as detection material, will be installed 1285 kilometres away, deep underground at the Sanford Underground Research Facility in South Dakota.To simulate the neutrino signal at DUNE we use the latest configuration file for GLoBES provided by the DUNE collaboration [96], which assumes 6.5 years of running time in both neutrino and antineutrino modes.Our analysis includes disappearance and appearance channels, simulating both signals and backgrounds.The simulated backgrounds include contamination of antineutrinos (neutrinos) in the neutrino (antineutrino) mode, and also misinterpretation of flavors.In addition, we also consider the DUNE-High Energy (HE) flux configuration, which has been proposed to optimize the study of tau neutrinos at DUNE [97].In this option the beam will peak at higher energies in comparison to the nominal configuration.When considering the sensitivity of DUNE we will compare the sensitivity using the standard DUNE beam with the high energy beam [97].The χ 2 function for DUNE is the same as for T2K and NOvA, see Eq. (12).

IV. RESULTS
In this section we discuss the results from the analyses performed.Given the form of the neutrino oscillation probability provided in Section II B we expect a long-baseline experiment to be especially suited to bound the decoherence effects discussed in this paper.Note, however, that if the baseline is too long, oscillations can become fully averaged and any sensitivity to decoherence is lost.This explains why solar neutrino experiments have no sensitivity to the decoherence parameters of interest in our analyses1 .Due to our choice of E 0 in Eq. ( 6) experiments with low (high) energy neutrinos will be sensitive to models with negative (positive) energy dependence (encoded in the parameter n).
The results of our analyses for all the models under consideration are presented in Figs.3-9.In the left panels we show the 90% C.L. upper limits that can be obtained with current data from the different experiments described in Sec.III, for each n.In the right panels we show the sensitivities for JUNO and DUNE (using the standard and the high-energy flux configuration), always in comparison with the current strongest bound, obtained from the different experiments at each n and represented by the grey-shaded region.
In Fig. 3 we consider Model A (Γ ≡ Γ 21 = Γ 31 = Γ 32 ).As can be seen, for negative n the strongest bound is obtained from the analysis of KamLAND data.For n = 0 all accelerators produce similar results, while for positive n the bound obtained from MINOS/MINOS+ is the strongest by far.The bounds obtained from the analysis of Daya Bay and RENO data are very similar and always about one order of magnitude weaker than the bound from KamLAND.Note that the slope of the corresponding limits, as a function of n, is different for each experiment.This is due to the pivot point E 0 = 1 GeV chosen in Eq. (6).Notice that the choice of a different reference value would change the bounds obtained in this work by a factor . Reactors use energies much smaller and hence they are mostly sensitive to negative n.The energies relevant for T2K are close to E 0 and therefore the curve is nearly horizontal.NOvA and MINOS use larger energies and hence they obtain stronger bounds for positive n.This feature is particularly striking in the case of MINOS/MINOS+.In this case the overall bound is dominated by the analysis of KamLAND and MINOS data.In the right panel of Fig. 3 we compare the current best bound for each n with the sensitivities for JUNO and DUNE/DUNE-HE.We find that JUNO will be able to improve the bound for negative n, while DUNE will improve the bound for n ≥ 0. Note that the sensitivity for DUNE using the high-energy beam can improve or worsen when comparing with   6)).
the standard beam depending on n.Since the standard energy beam peaks at smaller energies, stronger bounds can be obtained for negative energy dependencies than with the high-energy beam.
The bounds and sensitivities for Models B and C (Figs. 4 and 5) are very similar to the ones obtained for Model A, for n < 0. As anticipated in Section II B we see that as long as Γ 21 and either Γ 31 or Γ 32 are active, the effect on the oscillation probability at reactors is similar to the case when all Γ are switched on, and hence the results are nearly the same.For positive n the bound is dominated by the analysis of data from accelerators.In this case we find that the sensitivity is better in Model C than in Model B, since Model C uses Γ 32 and Model B Γ 31 , which has a smaller effect on the oscillation probability as we saw in Section II B. The situation changes for model D  (Γ 31 = Γ 32 and Γ 21 = 0), shown in Fig. 6.Since the fast oscillations due to ∆m 2  31 and ∆m 2 32 are averaged at KamLAND, KamLAND has no sensitivity at all to bound this model.For the other experiments the bounds on this model are again very similar to the ones obtained for Model A. Note that since JUNO will be able to resolve the fast oscillations, the JUNO sensitivity does not disappear as shown in the right panel of the figure .For the remaining models only one parameter is active, Γ 21 (see Fig. 7), Γ 31 (see Fig. 8), and Γ 32 (see Fig. 9).In the first case we find that all experiments except NOvA are capable of bounding Γ 21 .It should be noted that for all experiments under consideration, except KamLAND and JUNO, the ∆m 2  21 -driven oscillation is not developed yet, and the effect of activating Γ 21 can be cancelled by  6)).modifying sin 2 θ 12 accordingly.These experiments are capable of bounding Γ 21 because of the prior from the analysis of solar data [3] imposed in our analyses.If we left sin 2 θ 12 free in the analysis the bounds in this model would become weaker or would eventually disappear completely as in the case of NOvA.In the case of KamLAND and JUNO the main oscillation channel is directly affected and therefore we could always set a bound.The other two cases (where we switch on only Γ 31 or Γ 32 ) are very similar.The differences in the bounds come from the fact that these models are either effecting the oscillation term which goes with ∆m 2  31 or the one with ∆m 2 32 , which have slightly different amplitudes, as already discussed.In the case of the reactor neutrinos only in a few cases a bound can be obtained, because the data can be fit well using only one oscillating term with a modified amplitude.However, experiments of next generation will be capable of probing all models for every energy dependence.
We summarize all bounds that we have obtained for each model in Table III.As can be seen, the most constraining limit is often obtained in KamLAND and MINOS, although other experiments also contribute in a few cases.TABLE III: Summary of results: each column shows the most constraining upper limit on Γ ij , in GeV, for each model (A to G) and value of n.We also clarify, within parenthesis, which experiment sets the bound (KL = KamLAND, R = RENO, M = MINOS/MINOS+, N = NOvA).good agreement on the bounds with the authors of Ref. [35] and in addition we provide the bounds for n = ±2.Our analyses provide the strongest bounds to date for negative values of n.It should be further noted that very strong bounds were obtained from a combined analyses of solar and KamLAND data in Ref. [31].However, these limits are model-dependent and in particular do not apply to a scenario in which only decoherence effects are included (as it is always the case in this paper), as was already argued in Refs.[35,79].We find that the bounds for negative n will be improved by about one or two orders of magnitude, depending on the model under consideration, at JUNO.Regarding n ≥ 0 the bounds obtained from MINOS/MINOS+, T2K, and NOvA, can be improved using DUNE and DUNE-HE.The capabilities of DUNE where already discussed in Ref. [54] focusing only on n = 0.For completeness, we note that very strong bounds on quantum decoherence can also obtained from atmospheric neutrinos [39].
As a final note, in this paper we have performed analyses of each experiment individually.We only imposed a prior on sin 2 θ 12 from solar neutrino measurements, and after finding that Kam-LAND's ∆m 2  21 measurement is robust under decoherence we added another prior on this parameter in the analysis of the remaining experiments.The measurements of the standard parameters at the other experiments might not be robust anymore.We have found, for example, that the measurement of sin 2 θ 13 at Daya Bay and RENO is affected, and therefore this parameter must be left to vary freely in the analyses of the remaining experiments.A truly combined analysis might improve a bit the bounds that we have obtained here.The measurements of the standard parameters might be more precise from a combination of experiments and might therefore break the degeneracies with the Γ ij leading to slightly stronger bounds.Since the improvement from such a complicated analysis is not expected to be very large, we decided to focus on the experiments individually.

Γ 2 ΓΓ 2 Γ 2 FIG. 1 :FIG. 2 :
FIG. 1: Neutrino oscillation probabilities for the baselines, energies and flavors relevant for several of the experiments considered in this paper, in the case of Model A and for various n.The left panel is characteristic of long-baseline accelerator neutrino experiments, whilst the right panel represents reactor neutrino experiments.

TABLE I :
The phenomenological models considered in this paper.