Testing quantum decoherence at DUNE

We address some theoretical issues of the quantum decoherence phenomenon within the neutrino oscillation framework and carry out various tests under DUNE simulated experimental environment. On the theoretical side, we provide a general expression for an invariant decoherence matrix under a quantum basis rotation. On the simulated experimental side, considering a rotation invariant and a non-invariant decoherence matrix, we study the impact on the fitting of the standard oscillation parameters, the sensitivity in the mass hierarchy and CP violation, combining the neutrino and antineutrino mode and all available neutrino oscillation probabilities channels. Furthermore, a sensitivity for the decoherence parameter of order $10^{-24}\text{GeV}$ at 3$\sigma$, is obtained for our best case.

In this paper we will present a detailed demonstration of the behaviour of the decoherence matrix under rotations and obtain the general shape of a (rotation) invariant decoherence matrix. Furthermore, we will probe the neutrino oscillation probabilities for an invariant and non-invariant decoherence matrix under simulated experimental conditions in the context of DUNE, which will have unprecedented sensitivity to the identification of the mass hierarchy and the measurement of the CP phase δ [70]. Considering the two aforementioned kinds of decoherence matrix, we will perfom various tests such as: the sensitivity to the decoherence parameter, the effects of decoherence in the measurement of oscillation parameters and in the sensitivity to the mass hierarchy and δ. We include the study of the degeneracy between Γ and δ, in agreement with our second motivation point.

A. Density matrix formalism
The description of the neutrino system weakly interacting with the environment using the open quantum system approach is given by: This is the well known Lindblad Master Equation, wherê ρ(t) is the neutrino's density matrix and H is the Hamiltonian of the system. The term D encloses the dissipative/decoherence effects and is written as: where ρ is an eight-dimensional column vector consisting of the ρ k , the components of the aforementioned SU(3) expansion, while the 8 × 8 matrices M H and M D encode the Hamiltonian and decoherence components of the same expansion, respectively. The dissipative/decoherence matrix M D , that contains all the decoherence parameters, has to be a symmetric, positivesemidefinite matrix and its entries should satisfy a set of inequalities (see [71] for a full list). Thus, we can get the oscillation probabilities P (ν α → ν β ) ≡ P ναν β that can be obtained via inner products where ρ να (t) is the time evolved state from an initial neutrino flavour ν α and ρ ν β is the final neutrino flavour ν β to be detected. Considering that the neutrinos are ultrarelativistic, we have t = L where L is the baseline.

B. Neutrino Hamiltonian, rotation and CP phase
The neutrino Hamiltonian in the VMB for a neutrino of energy E ν is given by where G F , n e are the Fermi coupling constant and electron number density, respectively. Here U is the PMNS matrix, parametrized as with U δ = diag(1, 1, e −iδ ). Considering that the matrix U 23 U † δ conmutes with diag(A CC , 0, 0) and the matrix U δ conmutes with U 12 and diag(0, ∆m 2 21 , ∆m 2 31 ) separately, the neutrino Hamiltonian can be re-written as follows where H V (0) is given by The Hamiltonian H V (δ) can be diagonalized with the unitary matrix U T (δ, φ 1 , φ 2 ) that rotates the mass eigenstates in the MMB into the mass eigenstates in the VMB.
Using Eq. (6) and the property that the diagonal matrix H M (δ) commutes with the unitary phase operator U † φ (φ 1 , φ 2 ) = diag(1, e −iφ1 , e −iφ2 ), being φ 1 and φ 2 arbitrary phases, the latter equation can be rewritten as follows: Now that H M (δ) is shown as a unitary transformation of H V (0), it follows that the eigenvalues of the former do not depend on δ. We thus write: Comparing Eqs. (9) and (10) we get This result implies that once we find a matrix U T (0, 0, 0) that diagonalizes H V (0), a general U T (δ, φ 1 , φ 2 ) that diagonalizes H V (δ) can be constructed.

C. Invariant decoherence matrix analysis
We can apply the results of the preceding section to the problem of finding a decoherence matrix invariant under U T . As mentioned in Ref. [67], the relationship between VMB and MMB decoherence matrices goes as follows: where the orthogonal matrix P performs the rotation from the VMB to the MMB and is related to how the SU (3) generators transform under U T . The matrix P (δ, φ 1 , φ 2 ) is defined as: where the expressions for P 0 = P (0, 0, 0), R(δ) and F (φ 1 , φ 2 ) are given in the appendix. Eq. (13) is the analog of Eq. (11) after performing the SU (3) expansion.
Although we are emphasizing the dependence of M M D on δ and the arbitrary phases φ 1 and φ 2 in Eq. (12), we must notice that this matrix, in general, may depend on the neutrino energy and the matter potential A.
Now we want to prove that, given an invariant ma- This situation is very convenient to deal with since the neutrino matter oscillation probability formulae in the presence of quantum decoherence are easily found from their corresponding expressions in vacuum. To prove this, we rewrite Eq. (12) explicitly as and set the invariance condition For the case δ = 0, the invariance condition for Imposing the following relation to the right hand side of the preceding equation together with Eq.(15), we get: In Ref. [67], it is found that for δ = 0, π the invariant matrix M V D is given by (17), we find that the corresponding invariantM V D for any δ is described by the following block diagonal matrix: From the latter equation we can calculateM V D for δ = π/2, 3π/2 obtaining: D. Adding Hamiltonian terms in the weak coupling limit At this point it is interesting to make a comparison between the approach given in [68] for defining invariant matrices, and the one presented here and in [67]. With the purpose of demanding energy conservation, the following conditions are imposed in [68]: [H V AC ,V j ] = 0 and [Ĥ M ,V M j ] = 0, being thatV j andV M j are the dissipator operators written in the VMB (for the pure vacuum case) and in the MMB, respectively. To achieve the aforementioned conditions, it is required thatV j changes when the matter potential termÂ is added to the vacuum Hamiltonian H V AC ; all of this in the VMB. The new dissipator operator is defined by T , depending on the matter potential and the neutrino energy, in a way that it recovers its vaccum form when it is rotated to the MMB, this is: This situation is incongruent with the underlying regime used for describing the neutrinos as an open quantum system, which is the Born-Markov approximation (BM). The Born approximation is applied when the environment interacts weakly with the neutrino system [72,73], while the Markovianity assumption implies the use of the Lindblad form (Eq. (2)) [74,75]. In this evolution, the system is modified by the environment and parameterized by the dissipator operators (i.e. the decoherence matrix), but not viceversa. Thereby, the environment is considered unaffected by the system and the addition of the matter potential term to the vacuum neutrino Hamiltonian should not change (or in a negligible way) the dissipator operatorV j .

E. Selected decoherence matrix models
Our analysis will be limited to diagonal decoherence matrices due to their relative simplicity. While analytical expressions are available for these matrices in the context of vacuum oscillations [53], oscillations in matter can be found through perturbative expansions. For a generic matrix we can make an expansion in the small quantities α = ∆m 2 21 /∆m 2 31 , θ 13 andΓ i = Γ i L, assuming the latter is <0.1. The lowest order contribution to the transition probability P νµνe is found to be where P νµνe is the standard oscillation transition probability. We also notice that only three decoherence parameters of such a matrix dominate the contributions to the probability. Furthermore, these terms do not depend on the neutrino energy or matter potential, leading to a probability shift across all energies -a property that will help us qualitatively explain our results. Using the current values of the standard oscillation(SO) parameters given in Table I, this shift turns out to be positive.
For such a matrix we do not have a transition probability formula akin to (25), however, we may still use (23) to get a qualitative understanding of our results referred to this model.

III. SIMULATION AND RESULTS
For DUNE simulations, we will use the optimized beam as described in the CDR [70] and make use of the files from Ref. [76]. We use an exposure of 150 kt MW year exposure each for the neutrino and antineutrino channels. Both the neutrino mode (Forward Horn Current -FHC) and antineutrino mode (Reverse Horn Current -RHC) are used for 3.5 years on each. We are also using both the appearance (ν µ → ν e ) and disappearance (ν µ → ν µ ) channels on neutrino and antineutrino modes. All NC and CC background rates discussed in detail in Ref. [70] Section 3.6.1 are included.  Throughout this section, we use the SO parameters from the global fit based on data available in January 2018 [77,78] assuming normal ordering (NO) and some modifications. The ranges provided are not symmetric about the fit values, but we will take them as symmetric by taking the average of the one-sided deviations. Due to the relatively unconstrained value of δ, we consider different values for this parameter with special emphasis on δ = −π/2, which is rather similar to the current best fit value. As such, no priors are assigned to δ. We summarize this information in Table I. For all subsequent statistical analyses we will use the GLoBES package [79,80], while oscillation probabilities are calculated using NuSQuIDS [81].
This study relies on the χ 2 test statistic to compare data generated by a set of "true" parameters ξ true against a hypothesis that assumes a set of "test" values ξ. Unless otherwise stated, the parameters ξ true will have the values in Table I. The values of δ and Γ, will be provided on a case-by-case basis. In the context of DUNE and GLoBES, χ 2 is defined as where the N i is the expected number of events in the i'th energy bin. When including priors, GLoBES modifies χ 2 by adding an extra term where the summation in j is over all parameters for which σ j = 0. In this and all subsequent definitions for χ 2 , it is implicit that Γ and γ are used interchangeably depending on the assumed model.

A. Determination of standard oscillation parameters
Our first test consists on performing SO fits on data that includes decoherence. The data is generated assuming the true parameter values in Table I, for δ true = −π/2 and Γ true and no priors. The test hypothesis is Γ = 0 and we obtain χ 2 min by marginalizing over the SO parameters, attaining its minimum for a set of parameters which we label as "fit". We project the χ 2 function on the θ 13 − δ plane where all unmentioned test parameters are fixed to their fit values. We plot our results in Fig. 1. As we increase the value of Γ true , the SO fit shifts towards higher values of θ 13 for both models. We can explain this behavior through Eq. (23) where: as long as we increase Γ the ν µ → ν e transition probability grows in an energy-independent way, which, if it is fitted under the SO assumption, can be misconstrued as a larger mixing angle θ 13 . Both plots have similar shapes, where model B has a more notable shift to higher θ 13 ; a feature that is explained by looking at Eq. (23), where the parameters Γ 3 , Γ 8 are non-zero and increase the value of the shift.
On the other hand, the general shape of the contour plots remains the same. We also made the same analysis using the FHC and RHC separately, where the former provided stronger constraints compared to the latter, because of the higher statistics of the former. Similarly, we found that the appearance channels restrict the parameter space better than the disappearance channels.
We point out that short baseline neutrino oscillation experiments can measure the value of θ 13 quite precisely and will not be confused with the SO+decoherence scenario, since ΓL is much smaller than in DUNE. The precision on θ 13 is also expected to improve by the time DUNE becomes operational. An interesting possibility arises: if DUNE's measurement of θ 13 is incompatible with those obtained from reactor experiments, then decoherence would provide an explanation for this discrepancy, particularly if DUNE's measured θ 13 is larger than the accepted value.

B. Effects of decoherence on constraining δ
To test the ability of DUNE to constrain decoherence parameters, we generate data assuming true values for δ and Γ and make a χ 2 plot on the δ, Γ plane. Our χ 2 in this case assumes δ and Γ as free parameters and all remaining parameters are fixed to their true fit values (there is no marginalization here). Our findings are displayed in Fig.  2.
If we assume Γ true = 10 −24 GeV, we find that the data is still compatible with Γ/10 24 GeV = 2.0, 4.3 and 7.0 (1.7, 2.9 and 4.3) for model A (B) at 1, 3 and 5σ, respectively. We note that model B's limits are more stringent than model A, a feature explained mostly by the increased contribution from the first order correction to the transition probability (see Eq. (23)). For this assumed value of Γ true , both models are able to exclude the standard oscillation scenario at the 1σ level.
Even in the presence of decoherence, δ is well constrained. In [67] we pointed out that the CP assymmetries for δ = ±π/2 were similar for low neutrino energies and it was possible to confuse these two values of δ in the presence of decoherence. We also mentioned that after taking the whole of DUNE's energy range into account there was no similarity remaining. The aforementioned point is consistent with our current results, where we found that our plots have no secondary contours appearing in the vicinity of δ = π/2.

C. Constraining the decoherence parameter
For this study the simulated data assumes Γ true = 0 and a fit is performed for a given test parameter Γ. Our test statistic is where we marginalize over all standard oscillation parameters. We show the results in Fig. 3. The sensitivity to Γ/(10 −24 GeV) for the model A is in the range 1.5-2.1, 5.4-6.7 and 8.2-10.3 at the 1σ, 3σ and 5σ levels. On the other hand, the mode of the invariant matrix, model B, has a sensitivity almost invariant in δ and it is stronger than one predicted for model A, as we have already seen it in Fig. 2. This is a result of the increased contribution from the first order decoherence term for model B, which does not depend on the CP phase; the effect of other δ-dependent terms in a perturbative expansion are suppressed, such as the Γθ 13 term present in Eq. (25). For either model, DUNE's sensitivity to Γ is superior to the previously reported KamLAND's 95% C.L. sensitivity of 6.8×10 −22 GeV [56] and the 90% C.L. sensitivity of 1.2 × 10 −23 GeV [68]. Our current results are also of the same order of magnitude as those reported in a recent IceCube study [69]. We do remark that the decoherence models used in these previous studies are slightly different than ours and act as benchmarks.

D. Mass ordering sensitivity
Mass ordering sensitivity is obtained by comparing generated data from an NO assumption with decoherence against an IO hypothesis without decoherence. When generating data, a value of Γ true is assumed and the test statistic is marginalized in all oscillation parameters while keeping Γ = 0 fixed, meaning that The sensitivities are presented in Fig. 4 for different Γ true . We see that, regardless of δ true and our assumed decoherence model, the sensitivity is well above 5σ for our chosen values of the decoherence parameter and in fact improves it. To explain this feature, we point out that the IO fit yielded similar values for the SO parameters for all δ true , with the fit value δ fit in the vicinity of δ = −π/2. This occured for all the Γ true that we analyzed. When we studied the IO fit, we noticed that the probability in the ν µ → ν e channel is lower, and well separated, than that corresponding to the NO data in the energy range 2 GeV-4 GeV. Being that, as long as Γ true increases, the NO transition probabilities do as well. In return, the resulting IO fit transition probabilities separate more and more from the corresponding NO assumptions from the data. This separation leads to the observed increase in sensitivity. The impact of Γ on the sensitivity is strongest as we approach δ true = π/2.

E. CP violation sensitivity
To obtain DUNE's sensitivity to CP violation, we adopt a definition similar to the one presented in [70] χ 2 CP = min[χ 2 (δ = 0, Γ = 0, δ true , Γ true ), and we marginalize over all unmentioned test parameters. We checked that, after marginalization, the oscillation parameters remain approximately the same, with the exception of θ 13 that settles at higher values to account for the increased event rates in the appearance channels due to decoherence. The IO alternative was never preferred during the minimization procedure of χ 2 . This observation is easily explained by Fig. 4 because DUNE can distinguish between the hierarchies with a large significance.
The data generated using NO has an additional fake CP violation contribution coming from the matter potential present in the term proportional to Γθ 13 , which is mistaken in an SO fit as CP violation. For this reason, even when δ true = 0, π, decoherence may cause us to reject the null hypothesis (no CP violation). Due to this decoherence induced CP violation, decoherence lets us reject the null hypothesis at a higher confidence level throughout all values of δ.

IV. SUMMARY AND CONCLUSIONS
We have tackled diverse aspects of the decoherence phenomenon in the context of neutrino oscillations. From the theoretical side, we derived a general expression for invariant decoherence matrices that are diagonal (under rotations from VMB to MMB). The advantage of this type of matrix is its simple implementation into the matter oscillation probabilities formula as clarified in [67]. On the other hand, we have probed an invariant and noninvariant decoherence matrix in the context of DUNE. We have achieved a sensitivity for the decoherence parameter of O(10 −24 GeV) at 3σ and for δ = −π/2, in the best case scenario (the invariant matrix model). A very interesting result is the presence of an interplay between Γ and δ, predicted at the theoretical level in [67]. Additionally we found that, if decoherence plus standard oscillation is embodied within the data, the pure standard oscillation fit tends to select higher values of sin 2 θ 13 compared to the pure oscillation case. Finally, we have also seen that the decoherence phenomenon mildly affects the DUNE sensitivity on the mass hierarchy, while the CP sensitivity is clearly improved but with the disadvantage that even a true lack of CP violation (δ true = 0, π) can be confused with CP violation.