Flavour composition and entropy increase of cosmological neutrinos after decoherence

We investigate the evolution of the flavour composition of the cosmic neutrino background from neutrino decoupling until today. The decoherence of neutrino mass states is described by means of Lindblad operators. Decoherence goes along with the increase of neutrino family entropy, which we obtain as a function of initial spectral distortions, mixing angles and CP-violation phase. We also present the expected flavour composition of the cosmic neutrino background after decoherence is completed. Decoherence is proposed to happen after the two heaviest neutrino mass states become non-relativistic. We discuss how the associated increase of entropy could be observed (in principle). The physics of two- or three-flavour oscillation of cosmological neutrinos resembles in many aspects two- or three-level systems in atomic clocks, which were recently proposed by Weinberg for the study of decoherence phenomena.


Introduction
The indirect evidence for the existence of a cosmic neutrino background (CNB) is a major achievement of observational cosmology [1][2][3]. Its direct detection remains an experimental challenge [4,5].
Active neutrinos are in thermal equilibrium in the early Universe. Their momentum spectra follow the Fermi-Dirac distribution, which is preserved after their decoupling-even after they became non-relativistic. Non-instantaneous decoupling and neutrino heating by electron-positron annihilation introduce distortions to the perfect Fermi-Dirac spectrum for each neutrino flavour [6,7]. These spectral distortions are smoothed out by neutrino oscillations, but do not disappear completely. These distortions in the neutrino spectra imply a small increase of the neutrino density, giving N eff = 3.046 in terms of effective number of neutrinos. To arrive at this result, one assumes that there are no primordial lepton-flavour asymmetries.
Neutrino decoherence (i.e., the irreversible loss of quantum coherence) is an active research topic [8] with particular attention to its manifestation in neutrino-oscillation experiments [9][10][11] . In this work we study the cosmic evolution and decoherence of active neutrino flavours from their decoupling (when the Universe had a temperature of ∼1 MeV) to the present time. We ask the question what happens to the neutrino mass states when they become non-relativistic at late times. We propose in this work that the neutrino mass states undergo decoherence when the heaviest and second heaviest states become non-relativistic. The reasoning behind this hypothesis is that mass states in the non-relativistic regime have very different group velocities and therefore travel on different paths in space-time. The isotropic and homogeneous universe described by the ΛCDM model is homogeneous only in a statistical sense. There are fluctuations of the gravitational potential that vary in position and time. Coherent states traveling through this space-time may accumulate enough differences in their paths to reach an irreversible trajectory, causing the loss of information and the respective decoherence. Corresponding to the decoherence, a time averaging of the phase of these states traveling through different paths in space-time would also result in a breakup of time reversibility. As a consequence of this hypothetical decoherence, there is an associated increase of neutrino entropy. We also present the small corrections to the flavour distribution of the CNB today, which were seeded by primordial spectral distortions.
To track the evolution of the flavour content of individual neutrinos of given momentum we make use of the Wigner density matrix in a basis {|v a }: where the entries ρ ab (t, p) are defined at time t and momentum p.
To obtain the distribution function for an ensemble of neutrinos specified by the flavour state label α, we must trace the projected density matrix ρ p and multiply by the appropriate momentum spectrum and spin degrees of freedom: where f eq α (p, T) is the equilibrium momentum spectrum at temperature T in which we also account for the spin degrees of freedom. The density matrix contains the flavour dependent spectral distortions, which are also momentum dependent. Since we are only interested in the flavour composition of spectral distortions, we focus on the normalized density matrix spanned by the basis {|v a }. Moreover, we drop the sub-index p that indicates the momentum dependence for simplicity, but all the equations are momentum dependent. For this work, the momentum dependence is trivial as the comoving momentum is conserved after neutrino decoupling.
The time evolution of the density matrix ρ is given by the von Neumann equation. In an expanding Friedman-Lemaître model, after neutrino decoupling and after electron-positron annihilation (i.e., for T < m e /3 ≈ 0.2 MeV), it reads 1 : where H is the free Hamiltonian. The density matrix and the Hamiltonian are both functions of the cosmic scale factor a and the modulus of the comoving neutrino momentum q ≡ ap. Both ρ and H are Hermitian and trρ = 1. The differential operator: with H denoting the Hubble expansion rate and with a = 1 today. The cosmological redshift z = 1/a − 1.
A pure state is characterized by trρ 2 = 1, while for mixed states trρ 2 < 1. Thus, any physical system has trρ 2 ≤ 1. Since the von Neumann equation preserves trρ 2 , it cannot describe the process of decoherence. This follows as the von Neumann equation is a direct consequence of the unitary time evolution described by the Schrödinger equation.
The density matrix is connected to the (von Neumann) entropy via: 1 We use units in whichh = c = k B = 1 and M P = 1/ √ 8πG is the reduced Planck mass.
The von Neumann entropy vanishes for pure states and is maximal for maximally mixed states. For a two (three)-flavour state system, the maximally achievable family entropy (following Boltzmann) obviously is S = ln 2 (ln 3). Here we ignore the spin and momentum dependence of neutrinos. The von Neumann equation conserves this entropy and does not describe decoherence.
Here we assume that different mass states in the non-relativistic regime travel through different stochastic gravitational potentials, which leads effectively to a phase averaging. A time average along the neutrino path will consequently cause a definite decoherence. This loss of information cannot be recovered, it is not a de-phasing or kinetic decoherence that can be undone by a precise measurement. However, we will not model the different paths that the neutrino states travel and neither perform the time average, instead we will introduce a decoherence operator in the von Neumann equation encompassing the gravitational environmental fluctuations. The methodology of using an effective operator to introduce environmental sources of decoherence in an open quantum system is commonly adopted, including for neutrinos [12][13][14].
The modified von Neumann equation with a proper operator to describe decoherence is given by: which we refer to as Lindblad Equation [15]. The L a are so-called Lindblad operators, arising from tracking or averaging environment dynamics. The decoherence term is responsible for the fact that the quantum system can develop dissipation and irreversibility and lose quantum coherence. As we will show shortly, the entropy for the density matrix ρ L , as well as for a time-averaged density matrix over irreversible pathsρ, is not constant anymore. Many aspects of oscillating neutrinos in the early Universe have been discussed before. The focus of those works may be in the interplay with the primordial plasma [16,17], the matter effect [7,18], the role of a lepton asymmetry [6,19,20] or the back-reaction effects in primeval nucleosynthesis [21,22]. Important for this work are studies of distortions in the neutrino spectral distribution [21,23,24]. These works rely on an approach similar to ours, but are not identical, especially regarding the exact solution for three-flavour oscillations in the cosmological context, with and without including Lindblad operators to account for decoherence. In [25] a spontaneous suppression operator is adopted to account for decoherence, but the focus is not the evolution of entropy. Importantly, most existing studies are concerned with the interplay of neutrino oscillations and the plasma, therefore their focus lies on the epoch of neutrino decoupling. To our knowledge, the epoch after neutrino decoupling has not been explored in great detail, especially the effect of the transition from the relativistic to non-relativistic regime in a statistically homogeneous space, which we argue is the source of decoherence proposed in this work.
Often a prescription based on wave packets is adopted to account for the process of decoherence [26,27], e.g., for supernova neutrinos [28][29][30][31] or cosmological neutrinos [32]. In the latter work it is argued that cosmological neutrinos should not be treated as a classical, collisionless fluid. The author of [33] addresses the increase of neutrino entropy due to decoherence of mixed massive neutrinos, but without solving the Lindblad equation or using wave packets. As we show below, our findings go well beyond the studies presented in [32,33].
Lindblad operators have been used more recently to account for various decoherence processes [34][35][36][37][38], as well as in the context of neutrino oscillations in laboratory experiments, and both for two [12,35] and for three neutrino flavours [13,14,39]. However, the use of Lindblad operators for cosmological neutrinos in this work is novel, as well as the exact solution for three neutrino flavours that is valid throughout the relativistic and non-relativistic regimes.
Using Lindblad operators, one can mimic the damping of mixed states caused by a hypothetical decoherence phenomenon and predict the net change in entropy in a mathematical rigorous way. We also argue that cosmological neutrinos could be used to improve our understanding of decoherence in analogue to the recently proposed studies of decoherence in atomic clocks [37]. This work is structured as follows. In the next two sections we present exact solutions for twoand three-flavour neutrino oscillations in cosmology, including decoherence. We show that suitable Lindblad operators give rise to decoherence and that time averaging would lead to the same result for the asymptotic density matrix, we discuss the validity and generality of this result. In Section 4 we study the increase of family entropy associated with decoherence. We conclude with a discussion of our findings.
Whenever we use values for neutrino mixing angles and mass square differences, we use the best-fit values from [40] Where the values are for normal (inverted) hierarchy. The exception is the Dirac CP-violation phase which we set to zero for the sake of simplicity, although we explore its effect in the Appendix A. The best-fit value for the CP-violation phase is 0.8π(−0.03π), but the phase is not measured at statically significant level. The convention adopted on the mixing angles follows the standard of the Particle Data Group [41]. The neutrino oscillation parameters adopted here have been obtained from neutrino oscillation experiments only [40] and are in agreement with an independent recent compilation [42] that also uses cosmological information.

Neutrino Oscillations: Two-Flavour Case
We start by a discussion of the two-flavour case. Although some readers might think that this is a trivial exercise, we think that it is useful to clarify the concepts and the method for a discussion of the three-flavour case. We obtain an exact analytic treatment of neutrino vacuum oscillations from the relativistic to the non-relativistic regime in an expanding universe, their decoherence and time and/or momentum averaging for arbitrary initial conditions, neutrino parameters and cosmological model parameters.
The physical state of a system with two neutrino flavours is described by a two-dimensional Hilbert space (factored with the corresponding spaces for the other physical degrees of freedom-neutrino spin and momentum). The space of all Hermitian 2 × 2 matrices is spanned by the unit matrix I and the Pauli matrices σ i , where the Latin indices i run from 1 to 3. We write for any Hermitian matrix M = M 0 I + M i σ i , where we sum over repeated indices. We have trM = 2M 0 and trM 2 = 2(M 2 0 + M 2 i ). Please note that hermiticity implies that the components M i are real numbers. Expressing the density matrix in this matrix basis, the von Neumann Equation (3) becomes: with ijk denoting the totally antisymmetric symbol. The trace condition gives ρ 0 = 1/2. So far we did not specify a basis for the neutrino states. The free Hamiltonian is diagonal for the neutrino mass states and reads: with: There is a local minimum for the mixing angle θ 23 at 0.427 0.034 0.027 with a difference of ∆χ 2 = 0.02 when compared to the global minimum.
Without restriction of generality we assume that 0 ≤ m 1 < m 2 . For m 1 = m 2 we find that H 3 = 0 and ρ is a constant matrix.
Flavour-mixing is described by a two-dimensional rotation (U † U = I), written as: The Hamiltonian in flavour space reads: Thus, the solution in flavour basis is: where the initial condition for the mass states is given by the rotated conditions in flavour space:

Exact Solution
Neutrino production and detection involves neutrino interactions, i.e., flavour states. Thus, for both the initial conditions and the late time values of ρ we are interested in the flavour basis. Nevertheless, the Hamiltonian is diagonal in the mass basis, giving rise to simple time evolution of Equation (3). We thus first study the time evolution in the mass basis using Bloch vectors. Once the time evolution is known, we can specify the initial conditions in flavour basis, transform them to the mass basis, evolve in time and finally transform back to the flavour basis. As shown below, this can be done analytically without any assumption on neutrino masses, momenta or cosmological model.
In the mass basis the von Neumann Equation (7) is simply: Thus, ρ 0 and ρ 3 are constants. Furthermore, we define: With the new variable: we find the exact solution: where A q and φ q are to be fixed by the initial conditions. We find that trρ 2 = 1/2 + 2(ρ 2 3 + A 2 q ). As we saw already, the first necessary condition for neutrino oscillation (a non-trivial evolution of the density matrix) is m 2 > m 1 . A second necessary condition is θ = 0, as for θ = 0 the mass basis agrees with the flavour basis and only non-trivial values for ρ 0 and ρ 3 could be generated (under the assumption that neutrinos can only be generated in a pure flavour state). As was shown above, both ρ 0 and ρ 3 are preserved in the mass basis and thus for vanishing mixing angles no neutrino oscillations occur.
As we show below, there is also a third necessary condition for the oscillations of neutrino flavours to happen: at least one of the ρ fi = 0, otherwise the oscillation amplitude A vanishes.
To find explicit expressions we first study how the Pauli matrices are transformed from mass to flavour space. The transformation in the other direction is obtained by θ → −θ. This allows us to fix the constants A q and φ q . In the following we drop the explicit indication of the q-dependence; we find: and: Finally, at x > x ini , we may express the most general solution in flavour space as: It is interesting to check that indeed: is a preserved quantity. As trρ 2 ≤ 1 for any physical state, we find that the initial conditions must satisfy the constraint: Since the components ρ i are real, it follows that all individual components must come from the interval [−1/2, +1/2], for any initial conditions including arbitrary lepton-flavour asymmetries. Thus, we also see that trρ 2 ≥ 1/2.
For the special case of maximal mixing, given by θ = π/4, we find: with:

Initial Conditions
If at the initial time all neutrinos are in one of the two flavour states, we have: with δ = δ(q) ∈ [−1, 1] describing the asymmetry of the initial flavours. This means we assume that nature does not produce flavour-entangled (mixed) states at neutrino decoupling. Expanding in Pauli matrices: and replacing in (19) and (20), we have: and finally: In the following we choose x ini = 0, without restriction of generality.

Decoherence
To describe the decoherence of a system of two neutrino flavours in vacuum in an expanding Universe, we make use of a single Lindblad operator [35] and decompose it in Pauli matrices: where we have in principle four amplitudes l i . Any Lindblad operator L a must be Hermitian (L † = L) to guarantee a monotonic increase of the von Neumann entropy and it must commute with the Commutation with the diagonal Hamiltonian (in the mass basis in vacuum) requires that l 1 = l 2 = 0, resulting in the decoherence operator: Without restriction of generality we put l 0 = 0, as it does not show up in the Lindblad equation. Thus, the components of the master equation with the Lindblad operator L become: whose solution in mass basis can also be found analytically (using dt = −dx/2H 3 ): and are easily rotated to the flavour basis: Please note that l 3 can be an arbitrary real function of x q . The integral in the exponent is negative definite if l 3 = 0 and thus gives rise to a damping of all non-diagonal components in the mass basis. We interpret the function l 3 (x) as the influence of the environment, the expanding cosmos, that leads to the decoherence of neutrino states.

Averaging
The time averaging of (35) and (38) produces the same asymptotic result as the decoherence process by the Lindblad operator. The final effect is the suppression of terms with time dependence. For time averaging the fast oscillations terms take asymptotic values ( cos(x) = sin(x) → 0 and sin 2 (x) = cos 2 (x) → 1/2). Using Lindblad operators, one can see that the terms which have a time dependence in the mass basis, Equations (43) and (44), are exponentially suppressed once decoherence starts, extinguishing the off-diagonal contributions corresponding to mixed states. In the flavour basis, the off-diagonal terms are driven to a constant value. We use the effect of the Lindblad operator to describe decoherence, setting to zero the contributions from ρ 1 and ρ 2 in the mass basis and then rotate to the flavour basis:ρ For this averaged system we have trρ 2 = 1 2 [1 + δ 2 cos 2 (2θ)], which is independent of time or oscillation phase since we performed averaging and it depends on the mixing angle because it undergoes mixing. The result is numerically equivalent to a system that lost all coherence. For the exact, non-averaged solution we have trρ 2 = 1 2 (1 + δ 2 ), which is independent of mixing angles and time, as one expects for a system that undergoes unitary (deterministic) time evolution. The difference of the traces of the squared matrices, a measure for the difference of coherence of averaged and microscopic states, is then trρ 2 − trρ 2 = − 1 2 δ 2 sin 2 (2θ). For maximal mixing (θ = π/4), the time-averaged solution with arbitrary initial condition becomes:ρ Thus, the density matrix of maximally mixing neutrinos does not depend on the mixing of the initial flavour distortion. The probabilities to find a neutrino in the first or second flavour state are equal, and the amount of mixing is constant. For initial conditions in which the neutrinos are in pure flavour states (ρ f1 (x i ) = 0), the time-averaged density matrix is proportional to the unit matrix, as one would expect.

Relation of Lindblad Formalism and Averages
We observe that the expressions (50) and (53) are identical to the expressions (46) and (49) asymptotically (ρ L fi ρ fi ), i.e., if the oscillation phase x is large enough then decoherence has happened. The averaged density matrix agrees with the microscopic density matrix after decoherence. It might be tempting to conclude that decoherence and averaging over time (or perhaps momenta) would be equivalent. This has actually be proposed in previous literature, see e.g., [35]. However, a closer investigation reveals that this is not the case.
Let us apply our calculation to experiments which test solar, reactor or atmospheric neutrinos, to revisit the arguments given in [35]. The neutrinos of interest in oscillation experiments are relativistic, i.e., q ≈ E and thus ∆(a, q) ≈ ∆m 2 /2E (with energy E), and propagate in non-expanding space [(1/H)d ln a = dt, a = 1]. We find the oscillation phase (17) becomes: where the distance L is the distance traveled by neutrinos at the speed of light. The decoherence term becomes: For the analysis of neutrino-oscillation experiments, one usually measures the number of neutrinos that were emitted in one flavour and are detected in either the same or the other flavour. This corresponds to a maximal initial distortion (e.g., δ = 1 corresponds to ρ(x ini ) = |ν 1 ν 1 |). The well-known result, including a decoherence term, for the probability to measure the second flavour is now easily recovered, Let us compare this result with the probability obtained in a description in which wave packets instead of Lindblad operators (i.e., l 3 = 0) are considered to describe the process of decoherence [26,35]. The distribution of the oscillation phase x is often assumed to be a Gaussian, the phase averaged probability of flavour oscillations is then given by: where the phase width of the wave packet is σ = (x − x ) 2 . The integral gives: with x = (∆m 2 /2) L/E . Comparing the equation above with Equation (57), it has been argued that they have the same structure once we identify the decoherence term with the wave packet dispersion: This result is the basis of the argued equivalence between Lindblad decoherence and Gaussian averaging [35].
However, this argument is inconsistent already at the level of mathematical assumptions. To go from (58) to (59), σ was assumed to be a constant, i.e., independent of the phase x. However, the identification (60) suggests that either σ = σ(L), inconsistent with the assumption, or the integral on the l.h.s. should be constant. To make that integral a constant we can assume that l 2 3 is proportional to a Dirac-delta distribution on the time interval [0, L/c], which would correspond to a sudden, but incomplete decoherence, i.e., the asymptotic state remains to show (damped) neutrino oscillations. We conclude that the two mechanisms, decoherence and phase (or time or momentum) averaging are not physically equivalent.

Neutrino Oscillations: Three-Flavour Case
For the three flavour case the space of Hermitian matrices is spanned by the unit matrix and the Gell-Mann matrices λ i [43] where the index i runs from 1 to 8. In this case, Equation (3) becomes: where f ijk are the usual structure constants of the Lie algebra su(3). By unitarity, we mandatorily have ρ 0 = 1/3. In the mass basis the Hamiltonian is diagonal and is given by: with: Without restriction of generality we assume that 0 ≤ m 1 < m 2 < m 3 for normal hierarchy and 0 ≤ m 3 < m 1 < m 2 for inverted hierarchy. Once again, for equal or vanishing masses the Hamiltonian is proportional to the identity matrix.
The flavour mixing matrix is now a three-dimensional rotation (U † U = I) and can be written as: where c ij = cos θ ij , s ij = sin θ ij and δ CP is a Dirac charge-parity (CP) violating phase. In the main body of this work we assume vanishing CP-violation and include some results for a non-vanishing value in the Appendix A. Both Majorana phases are irrelevant for the aspects discussed in this work. The Hamiltonian in flavour space is obtained as in Equation (12).

Exact Solution
While the initial states and the states observable by means of a particle physics detector are given in flavour basis, the von Neumann equation and its solution is most suitable formulated in the mass basis. This approach simplifies the system of equations since the Hamiltonian is diagonal in the mass basis: The diagonal components of the density matrix (ρ 0 , ρ 3 and ρ 8 ) are constant. The remaining off-diagonal components give rise to oscillating solutions. The oscillation frequency is determined by: combinations of the asymmetric terms of the Hamiltonian (H 3 and H 8 ). There are six equations, forming three independent harmonic oscillators of two levels each, where their frequency is given by: for simplicity, we define three different oscillation phases to account for the three sectors of oscillation: where the only three independent combinations are ij = 21, 31, 32. We find the exact solutions: The amplitudes A, B and C and phases φ 12 , φ 45 and φ 67 are fixed by the initial conditions. For arbitrary initial conditions, we find trρ 2 = 1/3 + 2(ρ 2 3 + ρ 2 8 + A 2 + B 2 + C 2 ). During the radiation dominated epoch all three neutrinos are relativistic and the oscillation phase can be approximated by . When the neutrinos become non-relativistic, the oscillation phases start to evolve differently with redshift. In Figure 1 we show how the three different oscillation phases start to deviate from the phase evolution of relativistic neutrinos during the matter dominated epoch. Thus, we compare to the redshift dependence of the relativistic case, As can be clearly seen in the figure, at the moment when the most massive neutrino becomes non-relativistic, which happens at z ∼ 100 for the cases considered, the phases start to evolve quite differently, until they evolve again in parallel in the non-relativistic regime. It is worth noting that the heavier the neutrinos are, the earlier the transition to the non-relativistic regime happens. In Figure 1 we adopted the lightest mass state to be massless. Thus, one can infer that the transition happens at redshifts z 100.
The Gell-Mann matrices are transformed from mass to flavour space according to Uλ i U † . Instead of presenting the most general solution, we restrict our attention to initial states relevant to cosmology.

Initial Conditions
Just before neutrino decoupling muon and tau neutrinos interact via neutral currents only, while electron neutrinos also experience charge current interactions. Thus, the muon and tau neutrinos are expected to decouple slightly before the electron neutrinos [21,44]. Neutrino oscillations started in the early Universe slightly before their decoupling (defined as the moment when the interaction rate equals the Hubble expansion rate) from the primordial plasma. Non-instantaneous decoupling and quantum electro-dynamical corrections lead to spectral distortions. At this stage, the universe was filled with free electrons and positrons but not anymore with muons or taus, which had long annihilated or decayed to lighter leptons or photons. Thus, electron neutrinos end up with a different distortion of momentum and total number density than the other two active flavours. Later, during electron-positron annihilation, extra distortions are produced and again with different branching ratios for electron neutrinos and muon/tau neutrinos. x 32 x Q x 21 Figure 1. Evolution of the neutrino-oscillation phases (x 21 , x 31 and x 32 ) as a function of cosmological redshift z for the three Gell-Mann blocks for normal (left) and inverted (right) hierarchy with the lightest mass state assumed to be massless, shown for neutrinos at the peak of their thermal distribution (q = 3.15T ν ). The mass-squared differences are inferred from the global analysis of neutrino-oscillation data [40].
We assume that electron neutrinos are created with a distortion in the density matrix denoted by δ(q) and that muon and tau neutrinos can be described by the same distortion, −δ(q)/2. We further allow for a difference in the spectral distortion of muon and tau neutrinos described by β = β(q). Without primordial lepton-flavour asymmetry we would expect that β = 0. Without loss of generality we start to count oscillations at the moment when we set the initial conditions and thus have The assumption of pure initial flavour states gives then φ 21 = φ 31 = φ 32 = 0. The initial conditions are: with δ ∈ [−1, 2], β ∈ [−α, α], with α = √ 3 − 3δ 2 /4 from the condition trρ 2 ≤ 1. Setting the off-diagonal initial condition null is equivalent to assume that nature does not produce flavour-entangled states, which is an approximation since neutrinos decoupling is not instantaneous.
Using the above initial conditions in Equation (68), we are left with the following initial values for the density matrix in flavour basis (spanned by Gell-Mann matrices): where unshown entries are null. Rotating to the mass basis, we obtain an exact solution of the von Neumann equation for the Gell-Mann coefficients of the density matrix: [sin θ 12 sin(2θ 13 ) cos(2θ 23 ) + 2 cos θ 12 cos θ 13 sin(2θ 23 )] sin(x 32 ), where the coefficients ρ 0 , ρ 3 and ρ 8 are time-independent. The corresponding result including non-vanishing CP-violation phase is shown in the appendix. We do not present the general expressions in the flavour basis because they are lengthy and for the purpose of calculating solutions after decoherence, the solution in mass basis is all we need. Instead we restrict our presentation to the case of the normal neutrino hierarchy with vanishing Dirac CP-violation phase and for the best-fit values of the mixing angles from neutrino-oscillation data:

Decoherence
As with the two-flavour case, the Lindblad operator for the three-flavour case has contributions from the the same basis elements as the Hamiltonian ([H m , L] = 0), therefore its form in Gell-Mann matrices is: Consequently, the decoherence term reads: Apparently, there is no contribution of the Lindblad term to the evolution equations of ρ 0 , ρ 3 and ρ 8 , which thus remain constant in the mass basis. As above, we set l 0 = 0. The Lindblad Equation (68) can be written as: Its solution in mass basis acquires the following decaying modes: where the first terms on the right side are identical to the solutions without the Lindblad operator, as in Equations (76) and (82). The terms ρ 0 , ρ 3 and ρ 8 are constants and equal to their initial condition given by Equations (75), (79) and (83).

Averaging
The procedure to obtain time-averaged solutions is identical to the two-flavour case. We consider that the Lindblad operator acts on the solutions in mass basis, suppressing the time-dependent terms of the density matrix (i.e. ρ 1 , ρ 2 , ρ 4 , ρ 5 , ρ 6 , ρ 7 → 0). Then the averaged density matrix in flavour basis is obtained by rotating the remaining time-independent terms (i.e. ρ 0 , ρ 3 , ρ 8 ) to obtain expressions similar to (50) and (53).
Applying the Lindblad operator simplifies the solution in the mass basis, but the rotation to the most general flavour basis generates solutions that are again too lengthy to include here. However, we can show simple solutions using the best-fit values for the mixing angles in normal hierarchy and vanishing CP-violation:ρ This allows us to obtain the probability P αα to find a neutrino of the CNB in flavour state α: A graphical illustration of this result is provided in Figure 2. It shows that today's CNB does not necessarily have a 1:1:1 mix of the three active neutrino flavours. In fact, the expected mix of neutrino flavours depends on the initial values of the spectral distortions δ(q) and β(q). In the standard (minimal) scenario with vanishing lepton-flavour asymmetries we have β δ = O(10 −4 ) and deviations from flavour equality are tiny. Figure 2. Probabilities to find cosmological neutrinos at different flavour states before (dotted) and after (solid) decoherence as a function of the initial distortions δ = δ(q) (left, β fixed to zero) and β = β(q) (right, δ fixed to zero). The initial distortion for electron neutrinos is δ/3, for muon neutrinos −δ/6 + β/3, and for tau neutrinos −δ/6 − β/3. The mixing angles are given by the global fit to neutrino-oscillation data [40] for the normal hierarchy and we assume a vanishing CP-violation phase.

Discussion
The Lindblad operator introduces two degrees of freedom for the rate of decoherence (l 3 and l 8 ) but the formalism itself does not indicate when the operator should start to act. To better understand the evolution of the three-flavour system, we present the increase of the oscillation phase as a function or redshift in Figure 1. We see that the phases evolve in the same way as long as the neutrinos are effectively relativistic. This may give a first hint that the decoherence of cosmological neutrinos should start when the heaviest neutrino becomes non-relativistic, as already discussed for the case of two neutrino flavours.
We have learned from the discussion of the two-flavour case that decoherence via a (in general time and momentum dependent) Lindblad operator and averaging lead to the same asymptotic density matrix. The propagation speed of a wave packet is given by the respective group velocity. In the case of three neutrino masses we have three different group velocities, which are functions of redshift: In Figure 3 we plot the difference of these group velocities for pairs of neutrino mass states for q = 3.15T ν . We observe that the group velocities are identical as long as all neutrinos are relativistic. They differ when the heaviest neutrino becomes non-relativistic. Again, this suggests that coherent neutrino oscillations, as described by the von Neumann equation, take place during the relativistic neutrino propagation. Without decoherence that statement would hold until today.
One could argue that neutrino oscillations can be averaged shortly after neutrino decoupling and thus decoherence takes place in the early Universe; however, it seems that there is no justification for that as neutrinos do not scatter at T < 1 MeV and after the annihilation of positrons and electrons the number of potential scattering partners of neutrinos drops by another factor of 10 9 . The fact that the oscillating phase assumes high values does not necessarily mean that averaging is due automatically, since in principle one could still recover its exact value and obtain the corresponding microscopic state and trace it back to the initial state. A mechanism able to distinguish the mass states is still necessary. We suggest here that it is the transition from the relativistic to the non-relativistic evolution that induces decoherence of cosmological neutrinos and it is the difference in the inertial mass of the different neutrino states that "couple" in a different way to space-time. In a perfectly isotropic and homogeneous space-time, the different phase and group velocities would lead to a different world lines for the three mass states, but instead of evolving as independent "classical" states, they would evolve as entangled states. We think that situation must be different when we consider that the space-time is homogeneous and isotropic in a statistical sense only. We suggest that in such a situation the entanglement of the states will probably be lost due to their (gravitational) interaction with the environment.

Two-Flavour Case
Let us now have a closer look at the family entropy of the system. For density matrices with |ρ i | 1/2, we approximate the von Neumann entropy, such that the logarithm can be Taylor expanded: where we use the well-known properties of Pauli matrices, trσ i = 0 and trσ i σ j = 2δ ij . The allowed range of initial conditions gives ln 2 − 1 2 ≤ S ≤ ln 2 at the leading order. Presumably, higher order corrections change that into 0 ≤ S ≤ ln 2.
We can also obtain an exact expression for the family entropy. In the flavour basis: For small initial spectral distortions, |δ| 1, we find: Thus, the family entropy is constant as quantum coherence persists as long as the neutrino oscillations are determined by the von Neumann equation.
Decoherence due to whatever reason leads to the increase of entropy. If we calculate the family entropy from the asymptotic solutions to the Lindblad equation, we find: while for small distortions: Therefore, the increase in family entropy caused by decoherence is: It is maximal for maximal mixing, which also results in the maximum entropy state. For maximal mixing and an expected spectral distortion δ of order 10 −4 , we find an increase of family entropy of the order of 10 −8 .
Let us note that the time average of the entropy without decoherence does not change, as the calculation of entropy and averaging do not commute. Thus, the mathematical identity of the time-averaged density matrix and the asymptotic density matrix after decoherence does not correspond to a physical equivalence of decoherence and time averaging. Nevertheless, this mathematical identity can be useful in the calculation and is exploited in this work.

Three-Flavour Case
The entropy for three neutrino flavours and small spectral distortions can be Taylor expanded: where we made use of trλ i = 0 and trλ i λ j = 2δ ij . For the initial conditions specified in the previous section, we find the simple result: We can now calculate the increase of entropy due to decoherence. It is possible to obtain an analytical result, dependent on the mixing angles θ 12 , θ 13 , θ 23 and the initial flavour distortions δ and β. For vanishing distortion between muon and tau neutrinos (β → 0) there is no dependence on the angle θ 23 or in the CP-violation phase (complete result in the Appendix A). This second degree of freedom (β) is responsible for distinguishing the third-level state, without it the system becomes identical to a two-level system, when mixing between second and third state is irrelevant, and it is not possible to develop CP-violation.
To calculate the increase in flavour entropy of the neutrino due to decoherence, we make use of the fact that the von Neumann entropy does not depend on the basis in which the density matrix is given; tr[ρ a ln(ρ a )] = tr[ρ b ln(ρ b )], where (a, b) are different bases. We choose the most suitable basis to calculate the entropy difference of the initial (coherent) states and the final (decoherent) states: The initial and final entropy is calculated using the flavour and mass basis, respectively, for both the density matrix is diagonal in the suitable basis. Here we exploit the mathematical identity of the asymptotic solutions for the density matrix from the Lindblad equation and the time-averaged solutions of the von Neumann equation.
This is an excellent approximation up to δ ∼ 0.5 and, contrary to the exact solution, shows a simple dependence on the parameters.
For vanishing primordial lepton-flavour asymmetries, the special case of identical distortion for muon and tau neutrinos is theoretically well motivated. Asymmetries between these two flavours are not expected when the spectral distortions were produced. The combination of the expected initial distortion (δ = 4.45 × 10 −4 , β = 0) [21] with the measured mixing angles [40] gives rise to an increase of flavour entropy, ∆S = 3.43 × 10 −8 . In Figure 4 we present an interesting non-trivial case, when each flavour has a different distortion and we show the dependence on the only mildly constrained mixing angle θ 23 . We find that the latter affects the cosmological predictions only weakly.  δ and β (left). The values for the mixing angles are taken from the global fit to neutrino-oscillation data [40], assuming a normal neutrino hierarchy and a vanishing Dirac CP-violation phase. In the right panel the distortion δ is set to 10 −4 and the mildly constrained mixing angle θ 23 is allowed to vary (shaded area is the allowed 3σ region).

Conclusions
In this work we have studied cosmological aspects of the decoherence of mixed neutrino states. We have described decoherence phenomena via Lindblad operators in the von Neumann Equation (6).
The evolution of the family composition of the cosmic neutrino background, from the time when neutrinos are decoupled until today, has been investigated in Section 3. We obtain the expected flavour composition of the CNB as a function of arbitrary initial spectral distortions, mixing angles and Dirac CP-violation phase. The net effect on the current flavour composition is shown in Figure 2, where we see that neutrino oscillations and the effect of decoherence tend to equilibrate the initial flavour composition. For the measured mixing angles, the equilibration is not perfect and a small residual flavour imbalance is expected. The remnant (in general momentum dependent) spectral distortion of electron neutrinos is expected to be of the order of 10 −4 in the minimal scenario (no primeval lepton-flavour asymmetry).
The PTOLEMY experiment [5] is a proposal to detect cosmological neutrinos by looking for electron kinetic energies beyond the end point of the tritium β-decay spectrum [4]. According to our result, the fraction of electron neutrinos in the CNB carries information on the initial spectral distortion after e ± -annihilation and consequently about the state of the universe at that time. One could even speculate about futuristic detectors with such an exquisite sensitivity that even CNB intensity anisotropies [45], similar to the cosmic microwave temperature anisotropies, would be detected. The residual imbalance of neutrino flavours in the minimal scenario is one order of magnitude larger than the expected CNB anisotropies (apart from the dipole). The flavour imbalance would be increased for a lepton-flavour asymmetric universe.
We obtained exact solutions for the time-dependent Wigner density matrix, valid for any mass and momentum in a homogeneous and isotropic universe. The use of Lindblad operators results in a similar phenomenology as time-averaging for the asymptotic density matrices. We demonstrated that explicitly in Section 2 for a two-flavour example.
While it is interesting that Lindblad operators and averaging give rise to the same asymptotic results for the density matrix, the time of decoherence and the details of its mechanism are not provided by the formalism used. In the three-flavour case we are left with two unknown, non-trivial and real functions, l 3 (x) and l 8 (x). To obtain physical intuition, we also studied the behavior of the oscillation phase and the group velocities and their differences in Figures 1 and 3. We found that the group velocities start to differ significantly once the heaviest neutrino mass state becomes non-relativistic. This suggests that the transition to the non-relativistic regime might trigger decoherence in the mass basis via the stochastic inhomogeneities of space-time, which are experienced differently by the different mass states. Subsequently neutrinos would propagate in non-degenerate mass states, which means that they are in a frozen mix of flavour states.
The problem of identifying the decoherence time can also be approached from an experimental perspective by looking for observables related to the decoherence process. Recently, a similar idea has been proposed by Weinberg [37], where he suggested that decoherence of an atomic three-level system in the context of atomic clocks could provide enough information on the time scale of decoherence that one could measure or at least constrain the Lindblad coefficients. Likewise, decoherence of cosmological neutrinos could produce a traceable observable signal by the increase of CNB entropy that immediately follows the decoherence process.
Such an entropy increase of neutrinos due to decoherence of mixed states was pointed out for supernova neutrinos [29] and for cosmological neutrinos [33] previously. The latter work assumed a fixed initial distortion proportional to the cross-section of each neutrino at the time when the oscillations started. We went beyond and considered arbitrary spectral distortions and mixing angles. To the best of our knowledge, the entropy increase with its dependence on the initial spectral distortions, mixing angles and CP-violation phase in the cosmological context is obtained for the first time. This result is valid regardless the time of when the decoherence process happens. Therefore, if we could track the CNB entropy as a function of time, we would determine the moment of neutrino decoherence observationally.
It remains to figure out how the CNB entropy could actually be measured: Since the CNB is isotropic, any increase in entropy could only manifest itself macroscopically as an induced dissipative pressure (bulk viscosity [46,47]). Thus, an entropy increase due to decoherence would affect the cosmic neutrino equation of state. Any change of the equation of state gives rise to a contribution to the integrated Sachs-Wolfe (ISW) effect. It is thus interesting to ask if such an effect could be large enough to be observable. If decoherence happens when the neutrinos become non-relativistic, i.e., when neutrinos develop different group velocities as shown in Figure 3, then the decoherence contribution to the ISW effect would show up at z 100. However, we expect the effect to be tiny, since it is not only suppressed by δ 2 ∼ 10 −8 , but also by the ratio of neutrino density to matter density at z 100. Together this gives an effect of order 10 −10 in the temperature anisotropies. A lepton-flavour asymmetric universe might however give rise to a larger entropy increase (see Figure 4).
In this work we focused on the minimal scenario and proposed a small residual flavour imbalance of the CNB and a tiny increase of neutrino entropy at z ∼ 100. A more detailed investigation of lepton-flavour asymmetric models might reveal useful constraints on the primeval flavour composition of the Universe.
We observe that the phase δ CP introduces a difference of phase for each solution as long as the initial distortion β is non-vanishing, otherwise it becomes a global phase that can be absorbed in the initial phase.
For non-vanishing CP-violation phase, the difference between probabilities for neutrinos and anti-neutrinos in non-trivial. In the case of the transition from electron to muon neutrinos it is given by: P eµ − Pēμ = sin(2θ 12 ) sin(2θ 23 ) cos 2 (θ 13 ) sin(θ 13 ) sin(δ CP ) which is consistent with the literature [48]. In order to stress the role of coherence in this particular result, we use the Gell-Mann coefficients for the solutions contained in Equations (96) and (101). We note that differences between neutrinos and anti-neutrinos for any transition probability vanish in the limit of completed decoherence.
We obtain the limit after decoherence is completed of Equations (A1) and (A9) and then rotate to the flavour basis. For simplicity we replace the mixing angles for the best-fit in the global analysis [40] for normal neutrino hierarchy: which are consistent with the solution for vanishing CP-violation phase present in Equations (102) and (110).
The change in entropy is also consistent with the case of a vanishing CP-violation phase in Equation (124), and identical for vanishing initial difference between muon and tau (β = 0). In Figure A1 we show the probability to find a cosmological neutrino in a specified flavour state in the averaged limit as well as the expected change in the entropy, both as functions of the CP-violation phase. Figure A1. Change in entropy (left) and flavour imbalance (right) after decoherence as a function of CP-violation phase. Initial spectral distortions δ and β are both fixed to the value 10 −4 . The values for the mixing angles are given by the global fit of neutrino oscillation data [40] for normal hierarchy.