Accounting for the Heisenberg and Pauli principles in the kinetic approach to neutrino oscillations

While oscillations of solar neutrinos are usually studied using the single-particle quantum-mechanical approach, flavor conversions of supernovae neutrinos are typically analyzed using the kinetic equation for the matrix of densities due to the necessity of including also the scattering processes. Using the Wigner formulation of quantum mechanics we show the equivalence of the quantum-mechanical and kinetic approaches in the limit of collisionless neutrino propagation (in a background medium). We also argue that solutions of the kinetic equation account for the Heisenberg uncertainty principle and the related effect of wave packet separation (for single neutrinos), as well as the Pauli exclusion principle, if the initial conditions are consistent with these fundamental quantum principles. Such initial conditions can be constructed e.g. by identifying the matrix of densities with the (reduced) single-particle Wigner function computed using initial conditions for the neutrino wave function. These constraints on the initial conditions may have an impact on the phenomenology of supernovae neutrinos.


Introduction
Neutrinos play an important role in the stellar dynamics and provide us with an invaluable tool to study stellar evolution. The detection of neutrinos produced by the Sun, pioneered by Ray Davis and his collaborators [1] and later verified by the Super-Kamiokande [2] and SNO [3] experiments, confirmed the Standard Solar Model [4] and led to the discovery of flavor neutrino oscillations. The observation of neutrinos from the supernova SN1987A by the Kamiokande II [5], IMB [6], and BNO [7] collaborations provided an insight into the process of the core collapse and the subsequent formation of a neutron star, as well as put constrains on properties of Beyond the Standard Model physics. The long-awaited [8] detection of new bursts is expected to provide further valuable information on the mechanism of supernovae explosion.
Neutrino propagation in the solar interior is almost collisionless and can be described by the Schrödinger equation for the neutrino wave function. Although scattering processes changing the neutrino momentum are practically irrelevant, coherent forward scattering on the ambient matter strongly affects flavor conversion of solar neutrinos through the MSW effect [9,10]. Furthermore, the effect of decoherence by wave packet separation, related to the neutrino momentum uncertainty and thus to the Heisenberg uncertainty principle, that is naturally accounted for in the quantummechanical approach, damps the oscillations. As a result, even a hypothetical experiment detecting one neutrino at a time would not observe a seasonal variation of the neutrino flavor composition.
On the other hand, in some phases of supernovae evolution neutrino collisions with particles of the ambient medium can play a dominant role. Due to the necessity of including also the scattering processes, flavor conversions of supernovae neutrinos are usually analyzed using the kinetic equation for the matrix of densities [11]. The oscillation term of the kinetic equation naturally incorporates coherent forward scattering and the related MSW effect. The Pauli-blocking factors in its collision term ensure that the exclusion principle is respected in the scattering processes. However, because the matrix of densities is a function of coordinate and momentum, the kinetic equation is often perceived as classical, i.e. describing an ensemble of particles with definite coordinates and momenta. For this reason, its consistency with the uncertainty principle and the ability to account for the effect of wave packet separation (at the level of individual neutrinos) are questioned.
The goals of the present work are twofold. First, we show that for collisionless neutrino propagation (in a background medium) the quantum-mechanical and kinetic approaches to neutrino oscillations produce equivalent results. Second, based on this observation, we argue that solutions of the kinetic equation account for the Heisenberg principle and the related effect of wave packet separation, as well as the Pauli principle, if the initial conditions are consistent with these fundamental quantum principles. To this end, we make use of the Wigner formulation of quantum mechanics and show that the form of the evolution equation for the (reduced) single-particle Wigner function matches the form of the kinetic equation for the matrix of densities. Constructing initial conditions for the matrix of densities using the neutrino wave function we also show that their form matches the form of the initial conditions for the (reduced) single-particle Wigner function. Therefore, (in the collisionless limit) the quantum-mechanical and kinetic approaches to neutrino oscillations produce equivalent results. Because the quantum-mechanical approach accounts for the uncertainty and exclusion principles, the same applies also to the kinetic approach provided that the kinetic equation is supplemented by appropriate initial conditions. This line of reasoning is elaborated on in sections 2 to 4. In section 2 we 'translate' the by now standard quantum-mechanical approach to solar neutrinos into the language of the single-particle Wigner function. Furthermore, we derive an evolution equation for the Wigner function and show that it accounts for the uncertainty principle if supplemented by initial conditions constructed from the neutrino wave function. In section 3 we generalize this analysis to N -particle neutrino systems. In particular, we derive an evolution equation for the reduced single-particle Wigner function and show that it accounts for the exclusion principle if supplemented by initial conditions constructed from the N -particle wave function. In section 4 we derive the kinetic equation in the limit of collisionless neutrino propagation (but including coherent forward scattering) and show that its form matches the form of the evolution equation for the (reduced) single-particle Wigner function. We also construct initial conditions for the matrix of densities using the neutrino wave function and show that they match those for the Wigner function. The obtained results are summarized in section 5. Finally, appendices A to D contain supplementary technical material.
2 Single-particle quantum-mechanical approach Thermonuclear reactions that power the Sun produce a large flux of electron neutrinos. Due to the flavor conversion that occurs during neutrino propagation in the Sun, in vacuum, and in the Earth neutrino flavor composition measured by terrestrial experiments differs from that at the production point. As the neutrino propagation in these three environments is (almost) collisionless, its evolution is well described by the Schrdinger equation. The Wigner function constructed from solutions of the latter can be represented as a convolution of a shape and a phase factor. The shape factor describes position of the neutrino wave packet and satisfies the Liouville equation. The phase factor describes conversion of the neutrino flavor and satisfies a Schrdinger-like equation. A rather accurate estimate of the phase factor can be obtained taking into account that for the experimentally known neutrino parameters: i) mixing angle θ 13 and the medium corrections to θ 13 are small; ii) neutrino propagation in the Sun is close to adiabatic; iii) neutrino propagation in the Earth is close to adiabatic within the layers, but the adiabaticity is maximally violated at the boundaries between the layers. Due to the effect of decoherence by wave packet separation solar neutrinos arrive at the Earth as an incoherent superposition of the mass eigenstates. For this reason their flavor composition does not exhibit a seasonal variation. Crossing the Earth surface, the neutrino mass eigenstates split into the matter propagation eigenstates and start to oscillate. This leads to a small difference between the flavor composition of the night-time and day-time neutrinos. Below we briefly review the single-particle quantum mechanical approach to solar neutrinos paying particular attention to the interplay of the shape and phase factors and the averaging prescription, used to take into account details of the neutrino energy spectrum and of the spatial distribution of the production points.
Solar neutrino fluxes. The Sun is a main-sequence star. Roughly 99% of its energy is produced in the so-called pp-chain, while the remaining 1% is produced in the CNO cycle [12]. Both chains transform protons into helium simultaneously creating electron neutrinos, see e.g. figure 10 of reference [13]. The neutrino fluxes produced by the individual branches of the pp-chain and CNO cycle (the latter recently observed by the Borexino collaboration [14]) are predicted by the standard   solar model [4]. Their radial distribution and energy spectra, calculated in reference [15], are shown in figure 1.
In addition to the energy spectrum characterizing the entire neutrino ensemble, individual neutrinos are characterized by the size of their wave packets. The latter is determined by the production process and differs for different flux components. For the pp and 7 Be neutrinos σ x ∼ 6 · 10 −8 cm [16]. On the other hand, for the 8 B neutrinos the estimated size of the wave packet is a factor of two larger, σ x ∼ (1 − 2) · 10 −7 cm [16]. This corresponds to the energy scale σ p ∼ 100 eV.
Formal solution of the Schrödinger equation. Because propagation of solar neutrinos is almost collisionless, its is well described by the Schrdinger equation. In the neutrino mass basis it reads x) and the kinetic energy operator given by . Because the width of the neutrino wave packet in coordinate space is usually much smaller than the length scale for the variation of the matter potential in the interior of the Sun or the Earth, one can approximate neutrino propagation by propagation in a spatially homogeneous but timedependent potential V ij (t). The value of this effective matter potential at a time t is determined by the potential at the position the neutrino wave packet reaches by that time. This approximation is widely used in the literature. Let us emphasize that it is used in the present section to analyze solar neutrinos, but is not relied upon in sections 3 and 4. For a homogeneous but time-dependent potential the momentum-representation counterpart of the Schrodinger equation (2.1) reads where H ij (t, p) = δ ij ω i (p) + V ij (t). As can be verified by substitution, its (formal) solution reads where t P is the neutrino production time and the time evolution operator is given by The Sun produces only electron-flavor neutrinos. The neutrino flavor states are linear combinations of the mass eigenstates, ψ α = U αi ψ i , where U is the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix [17,18]. Expressions for wave packets of the neutrino mass eigenstates have been derived in reference [19] by analyzing the neutrino production in the QFT framework and assuming that the external particles are described by Gaussian wave packets. In particular, it has been found that the (effective) momentum uncertainties of the produced neutrino state are different in different directions. However, because the resulting flavor composition of the day-time and night-time neutrinos does not depend on the exact shape of the neutrino wave packets, here we follow reference [20] and neglect the difference of the momentum uncertainties in different directions as well as the difference in the shape of ψ i (t P , p) for different mass eigenstates. In this approximation the initial wave function is given in the flavor basis by where σ p is the effective momentum uncertainty of the produced neutrino and p w -its characteristic momentum.
Single-particle Wigner function. The collision integral, that describes neutrino interactions at a particular point x of the detector, contains an integration over the neutrino momentum p. To treat x and p on equal footing it is convenient to describe neutrino propagation and detection in terms of the Wigner function [21]. A single-particle Wigner function, written in terms of the momentumrepresentation wave function, reads Substituting the formal solution equation (2.3) into equation (2.6) we arrive at To leading order in the gradients H(t, p ± ∆/2) ≈ H(t, p) ± ∂ p H(t, p)∆/2 where, given the form of the Hamiltonian, ∂ p H(t, p) = ∂ p ω(p) ≡ v(p). Because solar neutrinos are ultrarelativistic, v(p) can be approximated by the unit vector in the direction of the neutrino propagation, which we denote by v p in the following, the corrections being of the order of m 2 /p 2 . In this approximation the velocity matrix is proportional to the identity matrix and hence commutes with the Hamiltonian. From equation (2.4) it then follows that U(t, t P , p ± 1 2 ∆) ≈ U(t, t P , p)e ∓ i 2 ∆vp(t−t P ) . The resulting approximate Wigner function factorizes into a product a phase and a shape factors [22], with the shape factor given by Whereas the phase factor U(t, t P , p) (the terms phase factor and evolution operator are used here interchangeably) depends on the matter profile along the neutrino trajectory, the shape factor is determined (in the ultrarelativistic limit, see reference [22] for an analysis of subleading corrections) only by the initial conditions. For the initial conditions specified in equation (2.5) the flavor-basis counterpart of the shape factor, g αβ = U αi g ij U † jβ , reads [22] g αβ (t, where σ x is the coordinate uncertainty defined by σ x σ p = 1 2 . At t = t P the shape factor is centered in the vicinity of the production point x P . Note that d 3 x d 3 p/(2π) 3 g αβ (t, x, p) = δ αe δ eβ , i.e. the shape factor is normalized to unity.
An experiment detecting neutrinos via the charged-current interactions [23,24] is sensitive to the individual components of the flavor-basis Wigner function, αα = U αi ij U † jα . For the initial conditions specified in equation (2.10) we find where U αe = U αi U ij U † je is the flavor-basis counterpart of the phase factor. On the other hand, an experiment detecting neutrinos via the neutral-current interactions [3] does not distinguish between the neutrino flavors and is sensitive to α αα . Unitarity of the PMNS matrix and of the phase factor implies that α αα (t, x, p) = g ee (t, x, p). The ratio P α (t, p) ≡ αα (t, x, p)/ α αα (t, x, p) = |U αe (t, t P , p)| 2 (2.12) gives the probability of detecting at a time t the momentum mode p of the neutrino wave packet in a flavor α. Unitarity of the PMNS matrix and of the phase factor implies that α P α (t, p) = 1.
Because energies of solar neutrinos are below the muon and tau production thresholds, in the remainder of this section we will focus on calculating P e (t, p) and ee (t, x, p). To emphasize that neutrino propagation in a spatially inhomogeneous potential is approximated here by propagation in a spatially homogeneous but time-dependent potential we will keep the time argument of P e instead of expressing it through the distance traveled by the neutrino wave packet.
Propagation in the Sun. The neutral-current neutrino interactions with matter are flavor-diagonal and only induce an overall phase in U, which does not affect |U ee (t, t P , p)|. The charged-current neutrino interactions with electrons produce a potential V e = √ 2G F n e , where n e is the electron number density [13]. The latter is predicted by the Standard Solar Model and is well approximated by n e /N A = 245 exp(−10.54 R/R ) [4], with N A being the Avogadro constant. The resulting Hamiltonian is given in the mass basis by H ij (t, p) = δ ij ω i (p) + U † ie V e (t)U ej . In the standard parametrization the PMNS matrix reads [25] U (θ 23 , θ 13 , θ 12 where O ij is the orthogonal rotation matrix in the ij-plane which depends on the mixing angle θ ij , and Γ δ = diag(1, 1, e iδ ), where δ is the Dirac-type CP-violating phase. Following the simple yet accurate approximation scheme developed in reference [26] we neglect the (1, 3), (3, 1), (2,3) and (3, 2) elements of the Hamiltonian, see appendix A for more details. This yields This Hamiltonian is diagonalized by an orthogonal transformation in the 12-plane with the rotation angle ϑ 12 determined by For ∆ω 21 (p) = V e (t)c 2 13 ·cos 2θ 12 mixing of the neutrino mass eigenstates is maximal, ϑ 12 → π/4, i.e. we encounter the celebrated MSW resonance [9,10]. The energies of the propagation eigenstates are given by Note that ∆E reaches its minimal yet nonzero value at the MSW resonance. This effect is known as avoided level crossing [27,28]. Some further details can be found in appendix A. Neutrino propagation in the Sun is very close to adiabatic [16,[29][30][31][32][33][34][35][36]. In the propagation basis the evolution operator is given in the adiabatic limit by exp[−iφ(t S , t P , p)], where and t S is the time when the neutrino reaches the solar surface. Projecting the propagation eigenstates on the mass eigenstates at t P and t S we arrive at the evolution operator in the mass basis, A formal derivation of equation (2.18) is presented in appendix B. At the solar surface the matter density is equal to zero and the Hamiltonian is diagonal in the mass basis, therefore ϑ 12 (t S , p) = 0. For low-energy neutrinos ∆ω 21 (p) V e (t)c 2 13 ·cos 2θ 12 and hence ϑ 12 (t P , p) is close to zero as well. That is, the matter effects are small and the oscillations proceed as in vacuum. On the contrary, for the high-energy neutrinos produced close to the center ∆ω 21 (p) V e (t)c 2 13 · cos 2θ 12 and therefore ϑ 12 (t P , p) ≈ π/2 − θ 12 , see appendix A for more details. For the electron-neutrino initial conditions, see equations (2.5) and (2.10), in the limit θ 13 → 0 the resulting Wigner function reads in the mass basis ij (t S , x, p) = δ i2 δ 2j g ee (t S , x, p). Thus, for θ 13 → 0 this fraction of solar neutrinos is adiabatically converted into the second mass eigenstate by the matter effects. For a non-zero θ 13 we obtain 22 ∝ c 2 13 and 33 ∝ s 2 13 , while 23 and 32 are proportional to c 13 s 13 . Hence, to leading order in θ 13 this fraction of solar neutrinos leaves the Sun as a superposition of the second and third mass eigenstates.
Day-time neutrinos. Using the composition property we can represent the evolution operator for a neutrino reaching the Earth surface during day-time in the form where t D is the moment when the neutrino reaches the Earth. Because in vacuum the Hamiltonian is diagonal in the mass-basis, Because two subsequent rotations, by ϑ 12 (t P , p) and by θ 12 , in the 12-plane combine to a rotation by θ P 12 ≡ θ 12 (t P , p) = θ 12 + ϑ 12 (t P , p) we obtain in the flavor basis Its absolute square contains non-oscillating contributions, as well as oscillating terms with phases e −i∆φ ij (t D ,t P ,p) . The collision integral that describes neutrino interactions in the detector contains an integration over the neutrino momentum, models the detector's momentum resolution function. In particular, the momentum resolution function filters out those neutrinos whose momentum lies below the detection threshold. Momentum resolution of realistic detectors is far below σ p . Hence, the resolution function is practically constant in the range σ p of momenta around the characteristic momentum p w and the integral can be approximated by f (p, p w ) d 3 k/(2π) 3 ee (t, x, k). The remaining integration over k yields the neutrino density matrix ρ ee (t, x). Taking into account that ϑ 12 (t, p) weakly depends on the ∼ σ p deviations of p from p w we obtain for the contribution of the non-oscillating terms, see equation where now θ P 12 = θ 12 + ϑ 12 (t P , p w ), and To estimate the contribution of the oscillating terms we expand ∆φ ij (t D , t P , p) around the char- In this approximation and for Gaussian initial conditions the momentum integration can be done analytically and yields terms proportional to Because the distance from the Sun to the Earth is more than 200 times larger than the solar radius, . For low energy neutrinos the matter effects are small and the second term is also well approximated by ∆φ ij (t S , t P , p) = ∆ω ij (p)(t S − t P ). On the other hand, the matter effects can affect ∆φ 12 (t S , t P , p) for high-energy neutrinos, see appendix C for more details. However, the first term alone is sufficient to ensure that |∂ p ∆φ ij (t D , t P , p w )| σ x . Hence, the right-hand side of equation (2.25) is exponentially suppressed. Thus, in a (realistic) detector unable to resolve momenta that differ by less than σ p the oscillating terms average to zero. This is a manifestation of the effect of kinematic decoherence by wave packet separation [37]. These arguments imply that it is sufficient to consider only the contribution in equation (2.23). The latter can be written in the form is the survival probability of the electron neutrino. Expressed in terms of the (effective) mixing angles it reads Note that the right-hand side of equation ( Night-time neutrinos. Crossing the Earth surface, the mass eigenstates split into the matter propagation eigenstates and start to oscillate [35]. Using the composition property we can represent the respective evolution operator as where t N is the moment when the neutrino reaches the detector after crossing the Earth. Recalling equation (2.20), we obtain in the flavor basis Its absolute square contains terms in which the propagation phases φ(t D , t P , p) cancel out, as well as terms proportional to e −i∆φ ij (t D ,t P ,p) . Because the neutrino trajectory in the Earth is at least four orders of magnitude shorter than the distance between the Sun and the Earth, φ(t N , t D , p) is negligibly small compared to φ(t D , t P , p). Hence, it can be ignored when analyzing the contribution of the e −i∆φ ij (t D ,t P ,p) terms. As has been argued in the preceding paragraph, the latter are exponentially suppressed due to the effect of wave packet separation and can be neglected. This , with the survival probability P e (t N ) given by The red curve depicts the night-day asymmetry (averaged over the neutrino production chains and regions) as function of the nadir angle η for the neutrino energy E = 10 MeV calculated using the PREM model [38]. The dashed blue curve denotes the night-day asymmetry computed in the two layer (mantle-core) approximation. The blue-shaded areas correspond to the trajectories crossing only the mantle, whereas the green-shaded areas correspond to the trajectories crossing also the core. Each shaded region marks one of the nine layers of the PREM model.
In the limit U(t N , t D , p w ) → 1, where 1 is the identity matrix, P e (t N ) → P e (t D ). As the Earth matter effects are small, instead of computing P e (t N ) it is convenient to analyze P e (t N ) − P e (t D ) instead. The evolution operator U(t N , t D , p w ) is determined by the matter density profile along the neutrino trajectory. The Preliminary Reference Earth Model (PREM) [38] distinguishes nine layers with the matter density slowly varying inside the layers and sharp density changes at the borders between the layers. The detailed density profile can be found in table II of reference [38]. Neutrino propagation within the layers is very close to adiabatic, whereas at the borders the adiabaticity is maximally violated [26,31,39]. Hence, U(t N , t D , p w ) can be computed by stacking up the time evolution operators for the propagation in the individual layers calculated using equation (B.4). The trajectory, and therefore also the number of the layers crossed, depends on the nadir angle η, see e.g. figure 1 in reference [40]. The night-day asymmetry [39,41], , computed numerically as a function of cos(η) for the neutrino energy E = 10 MeV is presented in the right panel of figure 2. Its shape can be qualitatively understood in the two-layer approximation, that distinguishes only between the mantle and the core. For a neutrino crossing only the mantle (0 ≤ cos(η) 0.84) where t D + = t D + is the moment shortly after the neutrino enters the Earth, and t N − = t N − is the moment shortly before it leaves it. To a good approximation, the matter density profile of the Earth is spherically symmetric, hence ϑ 12 (t N − , p) = ϑ 12 (t D + , p). The resulting expression for the difference of the survival probabilities reads [42] see appendix D for a few further details. As can be read off from the right panel of figure 2, equation (2.32) reasonably well reproduces the numerical results for 0.4 cos(η) 0.84. For a neutrino crossing also the core (0.84 cos(η) ≤ 1) where t 1 and t 2 (with t D < t 1 < t 2 < t N ) are the moments when the neutrino enters and leaves the core. The approximate spherical symmetry of the Earth implies ϑ 12 (t 1 − , p) = ϑ 12 (t 2 + , p) and The resulting expression for the difference of the survival probabilities reads is the jump of the diagonalization angle induced by the jump of the matter density at the border between the mantle and the core, see appendix D.
As can be read off from the right panel of figure 2, equation (2.34) acceptably well reproduces the numerical results.
Evolution equation for the single-particle Wigner function. As follows from the above considerations, the phenomenology of solar neutrinos is determined by an interplay of the phase and shape factors. The phase factor, equation (2.4), describes the neutrino flavor evolution and satisfies a Schrödinger-like equation, The shape factor, equation (2.9), describes position of the neutrino wave packet and satisfies the Liouville equation, As has been demonstrated above, supplemented by the mass-basis counterpart of the initial conditions ee (t P , x, p) = g ee (t P , x, p), where we have used hermiticity of the Hamiltonian, H * jk (t, p − ∆/2) = H kj (t, p − ∆/2), in the second line. Expanding the Hamiltonian around p, H(t, p ± ∆/2) = n (±∆) n 2 n n! ∂ n p H(t, p), and using ∆ n e i∆x = (−i∂ x ) n e i∆x we can recast equation (2.39) in the form In the approximation used above ∂ p H ≈ v p and therefore commutes with . Neglecting the n ≥ 2 terms, i.e. performing the first-order gradient expansion, we recover equation (2.37).
Evolution equation for propagation in an inhomogeneous potential. The analysis presented above is based on an approximation replacing neutrino propagation in an inhomogeneous potential by propagation in a homogeneous but time dependent potential. An evolution equation for the Wigner function of a neutrino propagating in a general potential can be derived directly from the Schrödinger equation as well. For H ij (t, x, p) = δ ij ω i (p) + V ij (t, x) the momentum-representation counterpart of the Schrödinger equation (2.1) reads where we have used hermiticity of the Hamiltonian in the third line. Using equation (2.42) we find Next, we expand the Hamiltonian, x, p) can be pulled out of the integral. Using ∆ n e i∆x = (−i∂ x ) n e i∆x as well as y m e −ipy = (i∂ p ) m e −ipy , and integrating over y and k we obtain Performing the first-order gradient expansion, i.e. keeping at most the first-order derivatives with respect to x or p, we can recast equation (2.46) in the familiar form [22], (2.47) In the approximation ∂ p H ≈ v p used above the first anticommutator on its left-hand side simplifies to v p ∂ x (t, x, p) thereby reproducing the velocity term of the Liouville operator in equation (2.37). A detailed analysis of the Liouville term can be found in reference [22].
Integration over neutrino spectrum and production points. So far we have studied evolution of individual neutrinos. In the detector contributions of neutrinos produced in the different production chains i and different distances R from the center add up. In an experiment detecting neutrinos via the neutral-current interactions the count rate in an energy bin [E, E + ∆E] is proportional to the total neutrino flux in this energy range, In an experiment detecting neutrinos via the charged-current interactions the differential neutrino flux is weighted by the survival probability P e and the count rate is proportional to Two sources of kinematic decoherence. The count rates can also be expressed directly in terms of the single-particle Wigner functions. In an experiment detecting neutrinos via the charged-current interactions the count rate is proportional to where l enumerates the solar neutrinos produced by the detection time t, V is the detector volume, and P is the momentum-space volume containing the momenta corresponding to the energy range [E, E + ∆E]. The volume integral selects from the sum over l only those neutrinos whose wave packets are located inside the detector at t. The wave packets corresponding to different l have different characteristic momenta p w , whose distribution reflects the neutrino energy spectrum, see figure 1. The order of the integrations and summations in equation (2.50) is chosen so as to compute first the contribution of each individual neutrino and then sum over the neutrino ensemble. In an experiment sensitive to the neutral-current interactions the count rate is proportional to As has been discussed above, if the momentum integration volume P , determined by the experimental setup, exceeds the size of the neutrino wave packet in the momentum space, characterized by σ p , then the wave packet separation leads to averaging out of the oscillations in the density matrix. As a result, flavor composition of the day-time neutrinos does not exhibit a seasonal variation. This applies even to a hypothetical experiment detecting one neutrino at a time. Crossing the Earth surface the neutrinos start to oscillate. Their path in the Earth is relatively short and insufficient for the wave packet separation to occur. Therefore, the hypothetical experiment detecting one neutrino at a time would observe an oscillation pattern. However, this does not guarantee that the oscillation pattern be observable in a realistic neutrino experiment. At the detection time t flavor of the neutrinos produced at the same time t P and at the same point x P but possessing different p w differs due to the momentum-dependence of U(t N , t D , p w ) and θ P 12 in equation (2.30). In the sum over l their contributions add-up destructively thereby "washing out" the oscillation pattern, unless the detector resolution ∆E is smaller than the inverse distance from the source to the detector [32]. That is, there are two sources of kinematic decoherence, both averaging out the oscillations: dephasing of different momentum modes of a single-neutrino wave-packet, related to the uncertainty principle, and dephasing of many neutrinos [43]. Technically both are related to the integration over the neutrino momentum and finite momentum resolution of the detector.
Furthermore, even the neutrinos characterized by the same p w but produced in different parts of the Sun arrive to the detector in different flavor states because the effective mixing angle θ P 12 depends on the matter density at the production point, see equation (2.15). This effect also contributes to the "washout" of the oscillation pattern [31].

N -particle quantum-mechanical approach
Given the low density of solar neutrinos and the resulting smallness of the neutrino-neutrino interactions, the individual neutrinos can be considered as independent. In this approximation wave function of the N -particle neutrino system is given by the Slater determinant. The anti-symmetry of the wave function under exchanges of any two particles allows the definition of a reduced singleparticle Wigner function, with the coordinates and momenta of N − 1 particles integrated over. In the limit of non-overlapping individual wave functions the latter is simply the sum of the individual single-particle Wigner functions. Hence, it automatically reproduces the summation prescription used in the single-particle quantum-mechanical approach. An evolution equation for the reduced single-particle Wigner function can be derived from the Schrödinger equation. Its form matches the form of the evolution equation for the single-particle Wigner function. Solutions of the evolution equation constructed from solutions of the Schrödinger equation are consistent with the Heisenberg and Pauli principles at any given time, in particular at t = t P . Turning this the over way around, a solution of the evolution equation with the initial conditions constructed from a solution of the Schrödinger equation is consistent with the Heisenberg and Pauli principles also for t > t P .
N -particle wave function. Due to the low density of solar neutrinos, their self-interactions can be neglected. Therefore, the individual neutrinos can be considered independent. In this approximation wave function of the N -particle neutrino system, that by the Pauli exclusion principle must be antisymmetric under a permutation of any two particles, is given by the Slater determinant. In the mass basis . (3.1) For non-overlapping individual wave functions the normalization factor is given by N = (N !) − 1 2 . For overlapping wave functions, i d 3 p/(2π) 3 ψ a,i (t, p)ψ * b,i (t, p) = 0 for a = b, the overlap contribution must be included into the normalization. One can show that for identical particles the overlap contribution is time-independent and so is the normalization factor, as expected. Note that a set of orthogonal wave functions can be constructed from ψ a,i using the Gram-Schmidt method. The flavor-basis counterpart of equation (3.1) is obtained by replacing the mass-basis wave functions by the flavor-basis ones, ψ α 1 ... α N (t, p 1 , . . . , p N N (t, p 1 , . . . , p N ).
Reduced single-particle Wigner function. The Wigner function of the N -particle system is defined analogously to the single-particle one [44], The antisymmetry of the wave function under a permutation of any two particles allows the definition of a reduced single-particle Wigner function with N − 1 coordinates and momenta integrated over [44], x, x 2 , . . . , x N , p, p 2 , . . . , p N ) . (3.5) Expressed in terms of the momentum-representation, N -particle wave function it reads, Using the flavor-basis N -particle wave function we find that, similarly to the single-particle Wigner function, the flavor-basis counterpart of the reduced single-particle Wigner function is given by (3.6) and integrating over the momenta p 2 to p N , we find that in the limit of non-overlapping individual wave functions the reduced single-particle Wigner function is related to the single-particle Wigner functions of the individual neutrinos by In the flavor basis we obtain a similar relation. Equation (3.7) implies that in an experiment sensitive to the charged-current interactions the count rate is proportional to Thus, the reduced single-particle Wigner function naturally reproduces the prescription used in the single-particle quantum-mechanical approach to integrate over the neutrino production points and spectrum. Similarly, in an experiment detecting neutrinos via the neutral-current interactions the count rate is proportional to which reproduces equation (2.51). Equation (3.7) together with equations (3.9) and (3.8) implies that the N -particle approach to neutrino oscillations simultaneously accounts for both sources of kinematic decoherence: dephasing of different momentum modes of a single-neutrino wave-packet, related to the uncertainty principle, and dephasing of many neutrinos.
Evolution equation for the reduced single-particle Wigner function. An evolution equation for the reduced single-particle Wigner function can be derived directly from the Schrödinger equation. Differentiating equation (3.6) with respect to time, subsequently using equation (3.2), and taking into account hermiticity of the Hamiltonian we obtain Let us first consider those terms of the Hamiltonian (3.3) that contain H laka (t, p a , k a ) with a > 1, H laka (t, p a , k a ) . . . δ l N k N (2π) 3 δ(p N − k N ) + . . . . Their contribution to the right-hand side of equation(3.10) reads, , k a , p a ) .
Renaming l a ↔ k a and p l ↔ k l in the second term, we observe that the right-hand side of equation (3.11) vanishes. That is, the H laka (t, p a , k a ) terms do not contribute to the evolution equation. The contribution of the remaining term is given by, where we have renamed k 1 → k and k 1 → k to emphasize its similarity to equation (2.43). Using equation (2.42) we can recast it in the form, (3.13) Applying the same transformations that led from equation (2.44) to equation (2.45) we arrive at (3.14) Performing the same steps that led from equation (2.45) to equation (2.46), we obtain an evolution equation for the reduced single-particle Wigner function. Its form matches the form of equation (2.46). Performing the first-order gradient expansion, we can recast the evolution equation in the familiar form, H(t, x, p), x (t, x, p) . (3.15) Supplemented by the initial conditions constructed by combining equations (3.1) and (3.6), equation (3.15) accounts for the Heisenberg uncertainty principle and the Pauli exclusion principle. This has long been known in the context of non-relativistic quantum-mechanical systems, see e.g. references [45,46] and references therein. On the other hand, unlike the Schrödinger equation, the evolution equation also admits classical solutions describing an ensemble of particles with definite coordinates and momenta, if supplemented by the classical initial conditions.

Kinetic approach
The kinetic equation for the matrix of densities is derived in the framework of perturbative quantum field theory in the Heisenberg representation. Its form matches the form of the evolution equation for the (reduced) single-particle Wigner function. For a single neutrino system the initial matrix of densities is given by the single-particle Wigner function evaluated at t = t P . For an N -particle neutrino system the initial matrix of densities is given by the reduced single-particle Wigner-function evaluated at t = t P . Because the form of and the initial conditions for the evolution and kinetic equations match, their solutions are also identical. This implies in particular, that the uncertainty and exclusion principles can be accounted for also in the kinetic approach to oscillations by considering initial conditions consistent with these fundamental quantum principles.
Kinetic equation for the matrix of densities. In astrophysical objects such as core-collapse supernovae or neutron-star mergers neutrino oscillations and collisions, including processes of neutrino emission and absorption, might be equally important. Processes changing the particle number are naturally accounted for in the second-quantized formalism. In a homogeneous system the particle number density is given by the expectation value of the operatorâ † j (t, p)â i (t, p) constructed from the ladder operators of the neutrino field. In a system with weak inhomogeneities one should consider the expectation value of the operator instead [11,43]. Let us emphasize that its arguments x and p are not coordinates and momenta of the field excitations. To derive an equation of motion for the operatorˆ ij (t, x, p) we differentiate equation (4.1) with respect to time. The time derivative of the ladder operators is determined by the Heisenberg equation, e.g.
As has been emphasized in reference [43], in the relativistic limit the first anticommutator on its left-hand contains the full flavor-dependent velocity structure of the Liouville term.
Single-particle initial conditions. At the production time t = t P the state vector in the Schrödinger picture, used in the single-and N -particle quantum-mechanical approaches, coincides with the state vector in the Heisenberg picture, used in the kinetic approach. This allows us to construct initial conditions for the matrix of densities corresponding to those for the neutrino wave function.
Equation (4.8) is valid for arbitrary initial conditions, including the single-particle ones. The state vector of the i'th mass eigenstate of momentum p readsâ † i (t P , p)|0 . The respective wave packet is constructed as p) is the i'th component of the neutrino wave function. Using orthogonality of the state vectors of different mass eigenstates, 0|â j (t P , k)â † i (t P , q)|0 = (2π) 3 δ(k − q) δ ij , we obtain the following expression for the one-particle state vector in the Heisenberg picture 'initiated' at t P , Taking the expectation value of equation (4.1) with respect to the state vector equation (4.10) we obtain an expression for the matrix of densities identical to the single-particle Wigner function, In other words, the constructed initial conditions for the kinetic equation match those for the evolution equation for the single-particle Wigner function. Given that the two equations have the same form they produce identical solutions. Hence, solutions of the kinetic equation account for the Heisenberg uncertainty principle if the initial conditions are consistent with this fundamental quantum principle. To provide an explicit example let us return to the Gaussian wave packet specified in equation (2.5). The resulting initial matrix of densities reads in the flavor basis (4.12) Note that αβ (t P , x, p) is localized neither in coordinate nor in momentum space, thereby reflecting the Heisenberg uncertainty principle. The dephasing of its momentum modes in the course of the neutrino propagation is the source of kinematic decoherence for single neutrinos. The authors of references [76,77] argued that, due to very high temperatures and densities, in supernovae the neutrino wave packet is very small, σ x ∼ 10 −11 cm, which corresponds to σ p ∼ 1 MeV. Note that in general the matrix of densities can attain negative values, just like the Wigner function.
N -particle initial conditions. The kinetic equation has originally [11] been formulated to study flavor neutrino evolution of multiparticle systems in the mean-field approximation. The state vector of an N -particle neutrino system is constructed as where ψ i 1 ...i N (t P , p 1 . . . p N ) is the N -particle wave function. For a system of independent neutrinos it is given by equation (3.1). The wave function of a system with strong interactions does not factorize into a product of the individual wave functions. In either case the Pauli principle requires the wave function to change sign under permutations of any two particles. For fermions the ladder operators anticommute andâ † i 1 (t P , p 1 ) . . .â † i N (t P , p N )|0 changes its sign as well. Therefore, the state vector remains invariant. Taking the expectation value of equation (4.1) with respect to equation (4.13), we obtain an expression for the matrix of densities identical to the reduced singleparticle Wigner function, ..l N (p + ∆/2, p 2 , . . . , p N )ψ * jl 2 ...l N (p − ∆/2, p 2 , . . . , p N ) . (4.14) In other words, the constructed initial conditions for the kinetic equation match those for the evolution equation for the reduced single-particle Wigner function. Given that the two equations have the same form they produce identical solutions. Hence, solutions of the kinetic equation account for the Heisenberg uncertainty principle as well as the Pauli exclusion principle if the initial conditions are consistent with these fundamental quantum principles. To provide an explicit example, let us consider the wave function specified in equation (3.1) in the limit of non-overlapping individual wave functions. The exchange terms become sizeable only at extremely high densities which are far above those present in supernovae [77], which justifies this approximation. The resulting initial matrix of densities is given by  ij (t P , x, p) are given by equation (4.12) with (in general) different p w and x P for different l. For a small σ p the spectrum of the neutrino ensemble is determined by the distribution of the characteristic momentum p w . On the other hand, a large σ p "superimposed" on the distribution of p w may leave an imprint on the observed neutrino spectrum. Equation (4.15) implies that the kinetic equation is capable of accounting for both sources of kinematic decoherence: dephasing of different momentum modes of a single-neutrino wave-packet, related to the uncertainty principle, and dephasing of many neutrinos. Their influence on collective neutrino oscillations has been addressed in references [78] and [79] respectively. The authors of reference [78] argued that collective neutrino oscillations can only take place if the product of the neutrino-neutrino interaction potential µ and the initial length P 0 of the global flavor "spin" vector exceeds the neutrino momentum uncertainty σ p . Therefore, the existence of the additional energy scale, related to the Heisenberg principle, may have an impact on the phenomenology of supernovae neutrinos.
Let us note that in the kinetic approach it is relatively straightforward to consider more complex initial states, for example the Glauber state addressed in reference [80] or states with nonzero particle-antiparticle correlations addressed in references [47,[81][82][83][84]. In the latter case the kinetic equation for the matrix of densities is supplemented by evolution equations for the pair correlators.

Summary and conclusion
The present work shows that for collisionless neutrino propagation (in a background medium) the quantum-mechanical and kinetic approaches to neutrino oscillations produce equivalent results. Based on this observation, we also argue that solutions of the kinetic equation account for the Heisenberg uncertainty principle and the related effect of wave packet separation (for single neutrinos), as well as the Pauli exclusion principle (for N -particle systems) if the initial conditions for the matrix of densities are consistent with these fundamental quantum principles.
To this end, we consider single-and N -particle neutrino systems and 'translate' the Schrödinger equation for their wave functions, ψ i (t, p) and ψ i 1 ...i N (t, p 1 , . . . , p N ), into an evolution equation for the (reduced) single-particle Wigner function ij (t, x, p). The term 'reduced' refers to the Nparticle systems and means that coordinates and momenta of N − 1 particles are integrated over. Both the single-particle and the reduced single-particle Wigner function satisfy the same evolution equation. Solutions of the evolution equation, supplemented by initial conditions constructed from the neutrino wave function, at any given time t reproduce the Wigner function calculated using ψ i (t, p) or ψ i 1 ...i N (t, p 1 , . . . , p N ) respectively. Because the wave function accounts for the Heisenberg uncertainty principle and, for N -particle systems, also for the Pauli exclusion principle, this applies also to the (reduced) single-particle Wigner function. This implies that the evolution equation for the neutrino Wigner function is capable of reproducing the standard results for the day-time and night-time neutrinos, including the impact of wave packet separation on the neutrino flavor composition. Because (in the collisionless limit) the form of the kinetic equation for the matrix of densities matches the form of the evolution equation for the Wigner function, this applies also to the kinetic equation, if the same initial conditions are used in the latter. Constructing the neutrino state vector using the wave function of the single-or N -particle system, we find that this requirement is automatically fulfilled.
The kinetic equation is a powerful tool for analysis of supernovae neutrinos, capable of describing neutrino propagation from the neutrinosphere, where collisions dominate, to the outer layers, dominated by oscillations. The Pauli-blocking factors in its collision term ensure that the exclusion principle is fulfilled in the scattering processes. As shown here, its Liouville and oscillation terms ensure that the matrix of densities remains consistent with the Pauli and Heisenberg principles if it was initially consistent with these principles. The Liouville and oscillation terms are responsible for both sources of kinematic decoherence: dephasing of many neutrinos, and dephasing of different momentum modes of a single-neutrino wave-packet. While the former is related to the energy spectrum characterizing the entire neutrino ensemble, the latter is related to the additional energy scale characterizing the size of the neutrino wave packets and thus to the uncertainty principle. The existence of this additional energy scale affects collective neutrino oscillations. Further analysis of the interplay between the many energy scales involved -the spectrum, the vacuum oscillation frequency, the matter potential, the neutrino-neutrino refraction potential, and the neutrino momentum uncertainty -may unravel new facets of the already very rich and complex phenomenology of neutrino oscillations in core collapse supernovae. An important prerequisite for this analysis is a detailed study of the shape and size of the neutrino wave packets in supernovae, for which only relatively crude estimates exist at present.
(A.9b) Although written in a slightly different form, equation (A.9) is identical to the one obtained in the flavor basis [16,77]. Note that ∆E (considered as a function of p) reaches its minimum at the resonance.

B Time evolution operator in the adiabatic limit
The time evolution operator satisfies a Schrödinger-like equation, As has been discussed in appendix A, in the approximation used here the Hamiltonian is diagonalized by a rotation in the 12-plane, E(t, p) = O † 12 (ϑ 12 (t, p))H(t, p)O 12 (ϑ 12 (t, p)). Using equation (B.1) we obtain the following evolution equation for the matrix O † 12 (ϑ 12 (t, p))U(t, t P , p): The term ∂ t O † 12 (ϑ 12 (t, p)) · O 12 (ϑ 12 (t, p)) ∝ ∂ t ϑ 12 (t, p) can be neglected in the adiabatic limit. Solution of the resulting approximate equation reads Taking into account the initial conditions, U(t P , t P , p) = 1, we finally arrive at

C Propagation in the Sun
Applying equation (B.4) to the neutrino propagation in the Sun we obtain where components of the diagonal matrix E can be read off from equations (A.7) and (A.9). As an overall phase cancels out in the Wigner function, it is convenient to subtract t S t PĒ (τ, p)dτ from the integral. The (1, 1) and (2, 2) terms of the phase exponent are then determined by ∆E. Hence, ∆φ 12 (t S , t, P , p) = − t S t P ∆E(τ, p)dτ and ∂ p ∆φ 12 (t S , t, P , p) = − t S t P ∂ p ∆E(τ, p)dτ . Using equation (A.9) we obtain [16,77,86] ∂ p ∆φ 12 (t S , t, P , p) = − t S t P dτ ∆ω 21 (p) − V e c 2 13 · cos 2θ 12 (∆ω 21 (p) − V e (t)c 2 13 · cos 2θ 12 ) 2 + (V e (t)c 2 13 · sin 2θ 12 ) 2 . (C. 2) The integrand vanishes at the MSW resonance and changes its sign after crossing the resonance. As has been emphasized in reference [16], this effect can partially alleviate the build-up of kinematical decoherence for neutrinos produced above the resonance.
The largest density jumps, from zero to ∼ 3.4 g/cm 3 and from ∼ 5.5 g/cm 3 to ∼ 9.9 g/cm 3 , occur within the first three layers (radius from 6371 km to 6346 km) and at the border between the mantle and the outer core (radius 3480 km) respectively. The density jumps between the other layers are considerably smaller. For the trajectories crossing more than three layers, the contribution of the short span occupied by the first three layers can be neglected. Therefore, the Earth can be approximated by a sphere of radius 6346 km with the density ∼ 3.4 g/cm 3 at the 'surface' and a density jump from ∼ 5.5 g/cm 3 to ∼ 9.9 g/cm 3 at the radius of 3480 km. In this approximation equation (D.8) can be used for all trajectories crossing the mantle, while equation (D.10) can be used for all trajectories crossing also the core. Numerically, this approximation yields acceptably accurate results, especially for the trajectories crossing many layers, see right panel of figure 2.