Revisiting quantum decoherence for neutrino oscillations in constant matter density

We re-examine the matter neutrino oscillation probabilities considering the decoherence phenomenon as a sub-leading effect. In this paper we point out the relevance of having the correct interpretation of the decoherence matrix in the different quantum bases, within the framework of neutrino oscillation probabilities in matter. Based on this treatment we develop an analytical formula for matter neutrino oscillation probabilities for three generations, with a range of application up to the decoherence parameter $\Gamma \sim 10^{-23}$ GeV. We observe that, due to decoherence, the amplitudes of the neutrino/antineutrino oscillation probabilities increase in an energy independent way. We also find that decoherence can reduce the absolute value of the CP asymmetry, relative to its value at the pure oscillation case. As a side effect we have found a degeneracy between the decoherence parameter $\Gamma$ and the CP violation phase $\delta$.

Even though all the evidence indicates that the neutrino oscillation relies on the interference between different neutrino mass eigenstates, the presence of an as yet undetected sub-leading mechanism is still possible. Within the context of new physics, there are various candidates that can coexist with oscillations induced by mass, that we will call from now on standard oscillations. Among them we have neutrino decay  , nonstandard interactions [33][34][35][36][37][38][39], decoherence in oscillations [41-53, 55, 56], and other new physics effects [57][58][59][60][61].
In particular, the general consequences of considering a quantum system in interaction with its environment are irreversibility and decoherence. The decoherence phenomenon tends to destroy the interference pattern, through the introduction of damping terms of the type of e −ΓL (where Γ is a decoherence parameter and L is the neutrino source-detector distance or baseline). It is also possible for this phenomenon to modify the oscillation frequencies through the appeareance of new coherence terms. It has been pointed out that the source of decoherence could be originated by strings and branes [62,63], as well as quantum gravity effects [64]. There have been several papers that have included dissipative effects into the neutrino system, treating this as an open quantum system, developing the oscillation probabilities formulae in two and three generations, for vacuum and matter [41][42][43][44][45][46].
The main goal of this paper is revisiting the treatment for obtaining a three neutrino oscillation formula in matter, when dissipative effects are included. In this framework, the neutrino Hamiltonian in matter can be written in the vacuum mass eigenstates basis (VMB) or in the matter mass eigenstates basis (MMB). The latter is the basis that diagonalizes the neutrino matter Hamiltonian, for a constant matter density. When solving the system in the MMB we have to check if the decoherence matrix we propose in this basis can be generated from a rotation of the corresponding one in the VMB, where the decoherence matrix is in fact defined. This very relevant detail has been overlooked in some papers by assuming that the decoherence matrix can be written in MMB as diagonal and energy independent. We will show that, in general, the latter assumptions are fulfilled only in a few cases.
Considering this rotation appropiately, we provide a three generation analytical formula valid for a decoherence parameter Γ ≤ 10 −23 GeV, that corresponds to an upper limit for ΓL of O(10 −1 ), with a baseline L of O(10 3 ) km, a source-detector distance compatible with long baseline scenarios. We study the behaviour of these probabilities, using the DUNE baseline and energy range [65], and also explore how a CP violation measurement would be affected due to the presence of decoherence.

A. General neutrino Hamiltonian
We work with N neutrino generations going through matter with constant density. Sterile neutrinos get an additional contribution from the neutral current matter potential that cannot be ignored, since the latter cannot be factorized out as a global phase as we do with the purely active neutrinos scheme. The neutrino Hamiltonian in the VMB for a neutrino of energy E can be arXiv:1711.03680v2 [hep-ph] 27 May 2018 written as: Here, ∆M 2 = Diag(0, ∆m 2 21 , . . . , ∆m 2 N 1 ) is the mass term and ∆m 2 ij = m 2 i − m 2 j is the squared mass difference.
The matter potential term is 2EG F n e and A N C = √ 2EG F n n are the charged current (CC) and neutral current (NC) potentials. The parameters G F , n e , and n n are the Fermi coupling constant, the electron and neutron number densities respectively. The mixing matrix U rotates the VMB mass eigenstates into the flavour eigenstates.
The Hamiltonian H V given above is non-diagonal, but is diagonal in the MMB and can be written as: where∆ ij becomes the effective value of ∆ ij = ∆m 2 ij /2E in matter.
We define U M as the matrix that rotates the MMB mass eigenstates into the flavour eigenstates and can be expressed in terms of effective mixing anglesθ ij and effective phasesδ. We also introduce the matrix U T , which relates operatorsÔ V in the VMB to operatorsÔ M in the MMB viaÔ where we have defined U T = U † U M , the matrix that rotates the mass eigenstates in the MMB into the mass eigenstates in the VMB.

B. Density matrix formalism
We will consider the neutrino system coupled with the environment and treated as an open quantum system. Thus, its evolution is described by the Lindblad Master Equation whereρ(t) is the neutrino's density matrix and H is the Hamiltonian of the system. The dissipative term D is written as where {V j } is a set of dissipative operators with j = 1, 2, . . . , N 2 − 1 for N neutrino generations. The presence of the operatorsV j causes the evolution ofρ to be non-unitary. The dissipative term must satisfy the requirements of complete positivity and a von Neumann entropy that increases with time. The second condition is achieved by requiringV j to be Hermitian [40]. The usual approach is to rewrite Eq. (4) by expanding all terms in the basis for Hermitian matrices, which consists of the identity operator I and the SU (N ) generators t i , with i = 1, 2, . . . , N 2 − 1. For this aim, we decompose all operatorsÔ, such asρ, H andV j , aŝ where the generators t i satisfy [t i , t j ] = i k f ijk t k and f ijk are the structure constants of SU (N ). Then, we geṫ being the elements of the matrix M H given by where ρ k and h i are the components ofρ and H, respectively, written on the basis of Equation (6). The matrix M D , which contains all the decoherence parameters, will be called the decoherence matrix. In general, it satisfies the following From the second property, it follows that the diagonal elements provide an upper bound to the off-diagonal ones. The inequalities mentioned in the third property will only be explicitly written for a few select cases to be discussed later. We point out that the entries of the decoherence matrix are sometimes assumed to be of the form (M D ) ij = γ ij E n , where typical values of n are -1,0, and 1. In this paper, we restrict ourselves to analyze the case where n = 0. The component ρ 0 = 1/N is constant in time and is only relevant when evaluating oscillation probabilities. The evolution of the ρ k can then be written in a compact formρ where ρ is an eight-dimensional column vector consisting of the ρ k . The solution to the differential equation is For a neutrino in a flavour eigenstate ν α , we can use Eq. (10) to find the evolved vector ρ να (t). Finally, the probability P (ν α → ν β ) ≡ P ναν β for a neutrino to be detected in the flavour state ν β is calculated via inner products We emphasize that for neutrinos, which are ultrarelativistic, we have t = L where L is the baseline.

C. Decoherence matrix relation between different quantum bases
When we take into account matter effects in the open neutrino system, it is expected that the decoherence matrix in the VMB (M V D ), where the decoherence parameters are defined, has a non-trivial relationship with the one written in the MMB (M M D ). The issue of moving to the MMB has been pointed out previously [46]. However, the description and the consequences of this relationship has been overlooked and not properly treated. We will get the connection between both decoherence matrices and explore the validity of different forms for the decoherence matrix in the MMB in the context of the neutrino oscillation system. We will also show that some decoherence matrices presented in literature are not allowed.
Starting from the Lindblad master equation and Eq.  (7), where the coefficients are now replaced with their corresponding ones in the MMB. We also point out that the properties of M D in section II B are not affected by unitary transformations.
The matrix M H given in the MMB is the simplest to deal with since the Hamiltonian is diagonal. For the decoherence term, we know that The labels V, M stand for VMB and MMB respectively. After a unitary transformation, the Gell-Mann matrices will be a new superposition of the generators t k where P kj will be an O(N 2 − 1) matrix. Using this substitution, we find that Doing a similar treatment to ρ M j and writing it in terms of ρ V j , we can prove that the decoherence matrix in the MMB is obtained by performing an orthogonal transformation with the matrix P From Eq. (15) a powerful property arises when the elements of the M M D are independent of the matter potential A CC , which implies that they are also independent of energy. In that case, any value of A CC , shall satisfy Eq. (15), in particular A CC = 0. As we can take P T A CC =0 = P | A CC =0 = I we have that: For instance, the trivial case that fulfills the condition above is when M V D is proportional to the identity I. When M M D = M V D the matter+decoherence oscillation probabilities can be directly obtained by replacing the standard oscillation angles and masses in the vac-uum+decoherence probability formula given in [43] with their effective values in matter. From now on, this substitution will be referred to as effective matter parameter substitution (EMPS). We point out that if Eq. (16) is not fulfilled, M M D is A CC and energy dependent. As we said before, the matrices M M D that satisfy Eq. (16), for the case of the potential A CC with constant matter density, would also fulfill this equation regardless the value of A CC . This observation implies that the aforementioned M M D matrices are the same even if A CC (x) depends on the position x or the non-constant matter density. In particular, we can conclude that these matrices satisfying the condition M V D = M M D could be used within the framework of the adiabatic or nonadiabatic case [42,56].
We must emphasize that all the results and discussions that will be presented in this paper have been developed within the scheme of a potential A CC independent of the position x.

III. TWO GENERATION MIXING
We first analyze the simple case of two generation mixing to illustrate the effects of rotating the decoherence matrix from the VMB to the MMB. The Hamiltonian in the VMB H V and the mixing matrix U are given by with ∆ = ∆m 2 /2E, and A = V CC /∆ = A CC /∆m 2 . In the MMB, we have the effective value of ∆, given by∆ and an effective mixing angleθ The matrix U M (U T ) is obtained directly from U by performing the substitution θ →θ(φ). The matrix P that performs the rotation from the VMB to the MMB in the SU (N ) basis is Using the notation in [42], we can rewrite P using the following correspondences where Taking a generic decoherence matrix with a, α and γ non-negative, we have the following inequalities that must be satisfied [42] 2R The corresponding matrix in the MMB is Given the form of P , we see that α is naturally not affected, while the effective decoherence parameters arẽ From this equation, we see that a diagonal decoherence matrix in the VMB, which is b = c = β = 0, can have off-diagonal entries in the MMB. Also, since the matrix P depends onθ, the decoherence matrix in the MMB is inherently dependent on E and A. It is interesting to note that some models provide a power law energy dependence to the decoherence parameters, which will not hold in the MMB. We continue by analyzing some features of 2ν oscillations in the presence of decoherence. We will also present a few numerical examples, which will use the standard oscillation parameters summarized in Table I. Decoherence parameters are assigned on a case-by-case basis.

A. Conditions for effective matter parameter substitution (EMPS)
From Eq. (25), we see that the only decoherence matrix that satisfies Eq. (16) for the two neutrino generations case is: Baseline L 1300km constrained by the condition α ≤ 2γ.
In order to exemplify what happens when the aforementioned replacement is misused, we take the following decoherence matrix, given in [46]: This matrix, written in the MMB is: Then, the identity given in Eq. (16), that allows us to make the EMPS, breaks down. This is ilustrated through three curves presented in Figure 1: one is the vacuum+decoherence oscillation probability; the other is the exact matter + decoherence oscillation probability, obtained numerically; the last one is doing the EMPS, this is the substitution ∆ →∆ and θ →θ into the vacuum formula P νµνe = sin 2 2θ (1 − e −at cos ∆t) /2. Comparing these curves, we see how the EMPS fails to describe the decoherence phenomenon, producing a fake peak in the probability around the energy region wherẽ θ = π/4 that allows maximal neutrino mixing. At lower energies we observe that the exact matter+decoherence formula and the direct susbtitution formula coincide very well. This is explained because at this part of the energy spectrum we have thatã → a,γ =b =β =c = 0, making the direct substitution valid for these energies.

B. Relevance of off-diagonal decoherence parameters
As we have pointed out in section II B, the decoherence matrix parameters are bound by the diagonal elements. Based on this, it is tempting to say that off-diagonal parameters can be neglected. In light of Eq. (25), we see that in the MMB the diagonal entries of the decoherence matrix also receive contributions from c.
We revisit this point to view the impact of off-diagonal parameters on oscillation probabilities, by considering the matrix proposed in [42] with b 2 + β 2 ≤ γ 2 /4, from which the following transition probability is obtained where we defined ν = ∆ sin 2θβ = β cos 2θ + b sin 2θ The interesting feature is that b, β affects the oscillation frequency. A simple extension where c = 0 is with the following constrains For this scenario, there is no simple form for the matrix exponential (Eq. (10) alternative is to assume that the decoherence parameters are small γ ∆, γt 1. Expanding to first order in the decoherence parameters P νeνµ = sin 2 2θ sin 2 Ωt 2 + γ + cµν where µ, ν are defined in (21) and Ω = (µ 2 + ν 2 ) 1/2 . Note that the parameters b, β do not contribute to the oscillation probability at first order. This feature is independent of the particular form of M V D . The vacuum limit is recovered by setting µ → −∆, ν → 0, Ω → ∆ and θ → θ We appreciate that the off-diagonal parameter c has a contribution comparable to the diagonal element γ even in the vacuum. In Figure 2, we present the probability when b = β = 0 and γ = 10 −23 GeV, which is appropriate for the DUNE baseline. The values chosen for c are zero and γ/2 and we display the results after propagation in vacuum and in matter. In both cases, it is clear that the contribution from c cannot be neglected. The discrepancy between c = 0, c = γ/2 starts from approximately 4.0 GeV and increases with energy.

C. Conservation of energy
In literature, the possibility of choosing a decoherence matrix such that energy is conserved has also been suggested. Conservation of energy implies that Tr(ρH) = 0 which requiresρ M 3 = 0 if we work in the MMB. This conditions reads For this to be valid at all t requiresc =β =γ = 0. We also need these three parameters to be E-independent, meaning that a = b = c = β = γ = 0 and α arbitrary (see Eq. (25)). The decoherence matrix must also satisfy the inequalities 2R ≥ 0 and 2S ≥ 0, given in Eq. (23), which is only achieved if α = 0 (the remaining inequalities are automatically satisfied). We therefore conclude that there is no decoherence phenomenon that allows for energy conservation in matter at all energies.
We must note that all Γ-independent terms (ρ (0) , ρ (θ) ,etc.) in the solution of (39) do not evolve in time and will contribute to the standard neutrino oscillation probability in matter. Although we do not write these terms explicitly in the probability formula, since the standard oscillation contribution is calculated numerically, these Γ-independent terms are still required to calculate terms of higher order that involve Γ. For example, we need to find ρ (α) to calculate ρ (Γα) .
The initial conditionρ(0) = ρ(0) = ρ ν β depends on U M only and is Γ-independent. Any term proportional to Γ (ρ (Γ) , ρ (Γα) , etc.) will therefore vanish at t = 0. After calculatingρ, oscillation probabilities are found using Given our method, the probability is also a power series in α, θ 13 and the elements in the decoherence matrix.

B. Decoherence matrix
Starting from a diagonal matrix M V D , with all its elements different, and demanding M V D = M M D in Eqs. (B2) and (B4) given in Appendix B we find that for δ = 0, π: and for δ = π/2, 3π/2: subject to the condition Γ 1 /3 ≤ Γ 2 ≤ 5Γ 1 /3. We have verified numerically that the diagonal cases above are the only ones that remain unchanged after being rotated into the MMB. Besides these cases we have the trivial one when M V D is proportional to the I. As we have already mentioned in Sec II C, only for these diagonal M M D matrices we can apply the EMPS to obtain the matter+decoherence oscillation probabilities.
We notice that the latter procedure has been misleadingly applied in recent papers. For instance, the M M D given in [46]: used in [71] do not depend on A or E and must follow the relation M V D = M M D . The matrix changes form when rotated back to the VMB and the resulting M V D is A-dependent. The latter does not make sense since decoherence, which is actually defined in the VMB, is an effect independent of the matter potential. Therefore the matrices M M D given in Eqs. (43) and (44) can not be derived from a realistic decoherence scenario. On the other hand, if we consider the matrices in Eqs. (43) and (44) as M V D , decoherence matrices defined in the VMB, they certainly meet all the physical conditions required to be a decoherence matrix in vacuum. However, they are not suitable for the application of the EMPS.
For simplicity, we will assume for our calculations the decoherence matrix in the VMB as diagonal. We assume one of the simpler forms which is subject to the condition Γ 2 ≤ 4Γ 4 .

C. Matter+decoherence oscillation probabilities
We set t = L and, in the context of DUNE, L = 1300 km is going to be the baseline of the experiment. As a first scenario, we will assume there is no decoherence and define ∆ = ∆ 31 L/2. DefiningΓ i = Γ i L, we obtain the neutrino oscillation probabilities, where the appearance probability is given by (46) and the survival probability is The probabilities P νµνe are, respectively, the appearance and survival probabilities in the absence of decoherence (Γ i = 0). These can be calculated numerically by any standard neutrino oscillation package. The antineutrino channels are found via the replacement δ → −δ and A → −A. We have also compared this result to numerical simulations, with an error lower than 5% in the energy range 0.3 ≤ E/GeV≤ 10.0 when Γ ≤ 10 −23 GeV for all channels, assuming the central values of the standard oscillation parameters. We therefore use these  formulas to describe new features arising from decoherence.
In this case, we see that Γ 4 does not contribute to the survival probability at the lowest order. Note that both the appearance and survival probabilities exhibit A and E-independent terms proportional to Γ 2 or Γ 2 2 . We will call these effects a shift. A second term proportional to Γ 2 2 is also present, which depends on sin(A∆), and becomes a matter-dependent term that vanishes for vacuum oscillations.
One of the more interesting features is CP violation induced by decoherence. Looking at Eq. (46), the CP violating terms are proportional to Γ 2 α, Γ 2 θ 13 , in addition to the standard CP violation arising from P (0) νµνe . In the limit A → 0, decoherence induced CP violating terms vanish at second order and the effect becomes subdominant.

D. Results
To present our results, we calculate the 3ν neutrino oscillation probabilities, including decoherence effects, using the nuSQuIDS package [69]. We compare the case of no decoherence with the decoherence bi-parameter model given in Eq. (45). We set Γ 2 = 10 −23 GeV and Γ 4 to the following values: 2.5 × 10 −24 , 10 −23 and 10 −22 GeV. These values of Γ 4 are chosen in order to satisfy Γ 2 ≤ 4Γ 4 , where Γ 4 = 2.5 × 10 −24 GeV is the limit value. Standard neutrino oscillation parameters are fixed to the central values in [65] (see Table II) with the exception of δ. Normal hierarchy (NH) is assumed throughout.
The appearance probabilities with the aforementioned set of values are shown in Figure 3, for both neutrinos and antineutrinos. As we can see the appearance probability has no significant change in shape after introducing decoherence. We do observe that the probabilities are higher in the presence of decoherence at all energies, whether it is the neutrino or antineutrino channel. The term responsible for this shift is the energy independent term proportional to Γ 2 which dominates the decoherence contribution appearing in our (46). It is important to recall that this formula is valid up to values Γ ∼ 10 −23 GeV. We also note that for values of Γ 4 between 2.5 × 10 −24 and 10 −23 GeV, this parameter has a negligible influence in the probabilities. Our analytical formula (46) reflects this since it has no decoherence contribution involving Γ 4 . On the other hand, when Γ 4 is raised to 10 −22 GeV, past the validity of our formula, the damping effect due to decoherence beigns to take control in the oscillation probabilities, diminishing its amplitudes. Additionally, using our numerical approach, we have tested the results given in [46]. We do not reproduce the different peaks exhibited in the plots given in this reference, and its probabilities do not have the energy independent shift predicted by our numerical and analytical result (see Figure 4 in [46]). Similarly to the two generation case, in [46] the effects of rotating the decoherence matrix in the three generation framework were not been properly treated. The latter is the root of the appearence of these strange peaks.
After assessing the decoherence effects on the probabilities, we look into its effects on CP violation. For this study we use the CP asymmetry defined as In Figure 4, we show four contour plots of the CP asymmetry A, in the plane δ and Γ, corresponding to four neutrino energies: 0.5, 2.7, 4.0 and 6.0 GeV. The latter values are representatives of the energy range of DUNE's neutrino beam spectrum. We set the decoherence parameters as follows: Γ = Γ 2 = 4Γ 4 . In each plot of this Figure, we have highlighted the iso-contours A produced by evaluating all the combinations between δ = {π/2, 3π/2} and Γ = {0, 1, 2} × 10 −23 GeV. We do not display the iso-contour that contains (3π/2, 0), given that this is practically a dot at this pair.
First, we observe a common feature of the asymmetry A for these energies: as we increase the Γ, the asymmetry decreases with respect to its values at Γ = 0 (the pure standard oscillation case), for any δ. For example, looking at 2.7 GeV, we see that for Γ = 0 and δ = π/2(3π/2), the asymmetry has values of 0.11(0.64). Meanwhile, for Γ = 2 × 10 −23 and δ = π/2(3π/2), the value of the asymmetry is 0.08(0.42), turning out in a decrement of the asymmetry by of 28%(35%) . Therefore, the tendency of this observable to vanish could be a signature of the decoherence phenomena. The degradation of the asymmetry can be explained by the term proportional toΓ 2 θ 13 in our approximate formula that was verified numerically.
Second, the iso-contour curves, which contain points with the same CP asymmetry A, reveal the existence of a degeneracy in the values of (δ, Γ). Taking, for instance, the 2.7 GeV curves, we appreciate that the iso-contour produced by the pair (3π/2, 2 × 10 −23 GeV) gives the same value of A as pairs like (1.12π, 0) and (1.84π, 0), which refer to the pure standard oscillation case. The value of A at these pairs are lower than the one at (3π/2,0). Likewise, for (π/2, 0) we note that the same value of A is obtained for (0.6π, 2 × 10 −23 GeV) and (0.36π, 2 × 10 −23 GeV). For this last case, the value of A is higher than the one at (π/2, 2 × 10 −23 GeV). Now we explore the possibility of lifting of these degeneracies for the following (δ, Γ) pairs: (π/2, 10 −23 GeV) and (3π/2, 10 −23 GeV), shown in Figure 5. In these Figures we displayed various iso-contours corresponding to different energies, 0.5 GeV, 2.7 GeV, 4.0 GeV and 6 GeV. We observe that for both (δ, Γ) pairs these four curves intersect only at the true point, thus solving the degeneracy. The clarity of the solution to the degeneracy problem in (δ, Γ) would be affected if we take consider a detailed experimental context. For the latter case we expect allowed regions around the true point and possibly around others (at some confidence level), being difficult to predict at what confidence level we may identify the true solution. An implication of this is the possibility of imposing upper limits in Γ.
In the latest global analysis of neutrino data, using the standard oscillation hypothesis, a hint for δ = 1.3π has been found [72]. If decoherence is present as a subleading effect, the latter result has to be taken with care, due to the degeneracy introduced by decoherence in the measurement of δ. In fact, it would be valuable to assess this distortion in the context of a simulation that convolutes the neutrino probabilities with the cross-section, fluxes, efficiencies, resolutions, etc, in facilities such as DUNE [73].

V. SUMMARY AND CONCLUSIONS
We have re-examined neutrino oscillation probabilities in matter in the presence of subleading decoherence effects. The effect of rotating from the VMB to the MMB was mentioned, pointing out that such a rotation inevitably changes the form of the decoherence matrix in the new basis. The inability of substituting oscillation parameters in the vacuum oscillation formula with their effective values in matter for arbitrary decoherence matrices has been heavily emphasized, providing strict conditions which must be satisfied for this method to be viable. In the context of three generation mixing, we have presented a perturbative approach to the decoherence problem, valid for Γ < 10 −23 GeV for the DUNE baseline, which explains the prominent features of the oscillation probabilities. More importantly, we show that a term proportional to Γθ 13 connects decoherence with CP violation. This term causes a (δ, Γ) degeneracy associated with a reduction of the CP asymmetry, when the decoherence parameter increases. We have shown that it is possible to lift these degeneracies, at the level of neutrino oscillation probabilities, combining a set of energies from the DUNE neutrino energy beam. However, to give a more precise answer, a realistic analysis is needed. Without a doubt, a future measurement of the CP asymmetry is going to be a useful tool for either bring to light the decoherence phenomena or constrain it. and X are expressions that we are not currently interested in, since they belong to the main diagonal. If M Γ is diagonal, it follows that Γ 1 = Γ 3 , Γ 4 = Γ 6 and Γ 5 = Γ 7 . Imposing this condition on M Γθ , we find that For this matrix to be diagonal, M Γθ must vanish, this turns out into three possibilities: 1. All Γ i are equal and M M D (or M V D ) are proportional to the identity.
plus the extra conditions given for a diagonal M Γ , relevant for the second and third case. The explicit expressions for M M D (or M V D ) are given in Eqs. (41) and (42) for δ = 0, π and δ = π/2, 3π/2, respectively.