Relating quantum mechanics and kinetics of neutrino oscillations

Simultaneous treatment of neutrino oscillations and collisions in astrophysical environments requires the use of (quantum) kinetic equations. Despite major advances in the field of quantum kinetics, the structure of the kinetic equations and their consistency with the uncertainty principle are still debated. The goals of the present work are threefold. First, it clarifies the structure of the Liouville term in the presence of mixing. Second, we derive evolution equation for neutrinos propagating in vacuum or matter from the Schrödinger equation and show that in the relativistic limit its form matches the form of the (collisionless part of the) kinetic equation derived by Sigl and Raffelt. Third, by constructing solutions of the evolution equation from the known solutions of the Schrödinger equation, we show that the former also admits solutions consistent with the uncertainty principle and accounts for neutrino wave packet separation. The obtained results speak in favor of a (quantum) kinetic approach to the analysis of neutrino propagation in exploding supernovae where neutrino oscillations and collisions, as well as the effect of wave packet separation, might be equally important.


Introduction
The disappearance of solar neutrinos reported by Ray Davis Jr and his collaborators [1] was the first hint of the existence of flavor neutrino oscillations. Shortly thereafter an expression for the survival probability of ν e was obtained by Gribov and Pontecorvo [2]. While propagation of the low-energy fraction of solar neutrinos detected by the Gallium experiments [3] is almost unaffected by the matter effects [4], forward scattering off the background matter strongly affects propagation of the high-energy fraction observed by the water Cherenkov experiments [5]. However, propagation of even the high energy neutrinos in the sun's interior is almost collisionless and can be described by the Schrödinger equation for the neutrino wave function ψ i (t, x) or by the equivalent Liouville-von Neumann equation for the respective density matrix, see reference [6] for a comprehensive review. As has been emphasized in references [7,8], this approach is consistent with the uncertainty principle and convenient for analysis of wave packet separation.

JHEP01(2020)138
On the other hand, neutrino collisions with particles of the ambient medium can play a dominant role in certain phases of supernovae evolution [9]. Because particle number is conserved in quantum mechanics, description of neutrino production and absorption processes typically relies on the Boltzmann kinetic equation for the neutrino occupation numbers. Their quantum-mechanical counterpart is the Wigner function: ̺ ij (t, x, p) = d 3 ye −ipy ψ i (t, x + y/2)ψ * j (t, x − y/2) . (1. 2) The Boltzmann equation in its classical form is not capable of describing neutrino oscillations. Thus, neither the Liouville-von Neumann nor the Boltzmann equation can consistently describe neutrino propagation from the neutrinosphere, where collisions dominate, to the outer layers, dominated by oscillations. The difficulty of "marrying" oscillations and collisions is evident already from the fact that the Liouville-von Neumann and Boltzmann equations operate with different objects, the density matrix and occupation numbers respectively. One of the first attempts to incorporate oscillations into the kinetic equations was made by Dolgov [10] who replaced, for each momentum mode, the neutrino occupation numbers by matrices in flavor space. Whereas diagonal entries of these matrices represent the usual occupation numbers, the off-diagonals correspond to coherence between the neutrino flavor eigenstates. In reference [9] Sigl and Raffelt derived kinetic equation for the matrices of occupation numbers featuring correct scattering, oscillation, and, last but not least, Liouville term, applicable in the physically relevant relativistic limit. However, their result is sometimes viewed with a degree of skepticism, mainly because the derived kinetic equation is believed to be in conflict with the uncertainty principle as well as unable to account for the effect of wave packet separation. A technically different approach to the analysis of neutrino oscillations was employed by Yamada in reference [11] and by Vlasenko, Fuller, and Cirigliano in references [12][13][14][15][16]. Instead of analyzing the density matrix or the matrix of densities they used Kadanoff-Baym equations for the neutrino Wightman propagator. Its quantum-mechanical counterpart is the two-point correlator: ̺ ij (t, x, ǫ, p) = dτ e iǫτ d 3 ye −ipy ψ i (t + τ /2, x + y/2)ψ * j (t − τ /2, x − y/2) . (1. 3) The obtained quantum kinetic equation is capable of describing neutrino flavor and spin oscillations, as well as neutrino collisions with the ambient medium. However, their results have not yet received the deserved attention in the community, mainly because of the high technical complexity of the used formalism. The goal of the present work is to establish connection between these three approaches to neutrino oscillations and present arguments in favor of the (quantum) kinetic approach to the analysis of neutrino propagation in supernovae, where oscillations and collisions, as well as the effect of wave packet separation, might be equally important. To this end we first consider neutrino oscillations in vacuum in section 2. In subsection 2.1 we briefly review the quantum-mechanical approach. In subsection 2.2 we derive evolution equation for the neutrino two-point correlator from the Schrödinger equation. Its form matches JHEP01(2020)138 the form of the quantum kinetic equation in vacuum. We also construct its solutions from the known solutions of the Schrödinger equation in vacuum and analyze the shell structure of the resulting two-point correlator. In addition to the mass shells ǫ = ω i we identify middle shells ǫ = (ω i + ω j )/2, that are responsible for neutrino oscillations. In subsection 2.3 we derive evolution equation for the neutrino Wigner function from that for the two-point correlator in the on-shell approximation (that necessarily includes the middle shells) and clarify the structure of its Liouville term in the presence of mixing. For relativistic neutrinos its form matches the (vacuum limit of the) kinetic equation of Sigl and Raffelt. Considering Gaussian and Lorentzian initial conditions we show that it admits solutions that are consistent with the uncertainty principle and account for the wave packet separation. In section 3 these results are generalized to collisionless neutrino propagation in matter. In subsection 3.1 we derive, to the first order in the gradient expansion, evolution equation for the neutrino two-point correlator valid also for nonrelativistic neutrinos. Studying it term by term we single out contributions relevant in the physically interesting relativistic limit. In subsection 3.2 we derive evolution equation for the Wigner function from the relativistic limit of the evolution equation for the twopoint correlator. Its form matches the (collisionless limit of the) kinetic equation of Sigl and Raffelt. Constructing the Wigner function from the formal solution of the Schrödinger equation we show that it admits solutions that are consistent with the uncertainty principle and account for neutrino wave packet separation also for neutrino propagation in matter. In section 4 we summarize our results. Finally, appendices A to E contain supplementary technical material.

Neutrino propagation in vacuum
In this section we consider neutrino propagation in vacuum and establish a connection between three technically different approaches to description of flavor neutrino oscillations. In subsection 2.1 we briefly review the quantum-mechanical approach operating with the neutrino wave function ψ i (t, x) or density matrix ρ ij (t, x). In subsection 2.2 we derive evolution equation for the two-point correlator ̺ ij (t, x, ǫ, p) from the Schrödinger equation. Its form matches (the vacuum limit of) the quantum kinetic equation. In subsection 2.3 we derive evolution equation for the Wigner function ̺ ij (t, x, p) from the evolution for the two-point correlator using on-shell approximation. In the physically relevant relativistic limit its form matches (the vacuum limit of) the kinetic equation of Sigl and Raffelt.

Description in terms of the density matrix
We first briefly review the standard quantum mechanical wave packet approach to neutrino oscillations. A comprehensive review and analysis of subtleties related to the neutrino production, propagation, and detection can be found in reference [6].
Schrödinger equation for neutrino wave function. Neutrinos produced in chargedcurrent weak interaction processes are considered to be in a flavor eigenstate ν α ∈ (ν e , ν µ , ν τ ) [17]. Up to small corrections [18] the latter is a linear superposition of the JHEP01(2020)138 neutrino mass eigenstates ν i ∈ (ν 1 , ν 2 , ν 3 ). In terms of the wave functions ψ α = i U αi ψ i . Evolution of the mass eigenstates is governed by the Schrödinger equation, (2.1) In vacuum H ij (t, x) = K ij (x) with the kinetic energy operator given by 2 (note that throughout this paper we work in the mass basis). As can be verified by substitution, resulting solution of the Schrödinger equation is given by [7] The initial conditions are encoded in ψ i (0, p). For the reason of computational simplicity it is common to assume that the mass eigenstates constituting the produced flavor state are described in the momentum space by Gaussian wave packets [7,8], where σ p is the momentum uncertainty of the produced neutrino state and p w -its characteristic momentum. Following reference [7] we neglect the difference in the shape of ψ i (0, p) for different mass eigenstates. As has been shown in reference [17] by analyzing neutrino oscillations in the QFT framework (see also references [19][20][21][22]), expression (2.3) is oversimplified and accounts neither for the difference of the momentum uncertainties in different directions, nor for the contribution of the energy uncertainty. However, it is sufficient for the purposes of this work, where it is used solely for illustration. For small σ p the corresponding coordinate-space wave function is well approximated by with v i (p) ≡ ∂ p ω i (p) being velocity of the respective mass eigenstate. Equation (2.3) states that the momentum-space wave function is localized in the vicinity of p w within the ∼ σ p range. Its coordinate-space counterpart is localized in the vicinity of x = vt within the ∼ σ x range, where the latter is defined by the uncertainty relation, σ x σ p = 1 2 . That is, in agreement with the uncertainty principle, delocalization in the momentum space results in localization in the coordinate space, and vice versa.
Liouville-von Neumann equation for neutrino density matrix. Since one is typically interested in the probability of detecting a neutrino of a particular flavor α, a convenient quantity for tracking the neutrino evolution is the density matrix defined in equation (1.1). Using the latter we find that the Schrödinger equation translates, without any approximations, into a closed-form Liouville-von Neumann equation for the density matrix, (2.5)

JHEP01(2020)138
Note that as in reference [23] ρ and H without generation indices denote matrices in the mass basis. Using the approximate solution equation (2.4) we obtain for the density matrix [7,8] (2.6) states that components of the density matrix are localized in the vicinity of x =v ij t within the ∼ σ x range, thereby reflecting the uncertainty principle. In addition, it implies that the off-diagonal elements of the density matrix experience a suppression proportional to the difference of the respective group velocities. This is a manifestation of kinematical decoherence [6,23] caused by wave packet separation: because wave packets of the neutrino mass eigenstates have finite spatial size ∼ σ x and propagate with different group velocities they separate in the course of neutrino propagation.

Description in terms of the two-point correlator
The definition of the density matrix ρ ij (t, x) requires the two underlying wave functions to be computed at the same point of the space-time, thereby reducing the number of 'degrees of freedom' from four (two time and two space arguments) to two (one time and one space argument). This reduction is avoided in the two-point correlator ̺ ij (t, x, ǫ, p) that trades the four arguments of the wave functions for another four -the central time, central coordinate, energy, and momentum -thereby providing a more granular description of the neutrino state than the density matrix.
As we demonstrate below, in vacuum the Schrödinger equation for the neutrino wave function translates, without any approximations, into a closed-form evolution equation for the two-point correlator. Its form matches (the vacuum limit of) the quantum kinetic equation.
The energy argument ǫ is independent of the momentum argument p. This means in particular, that the two-point correlator is not constrained to be on-shell. Of course, for typical initial conditions its diagonals do peak in the vicinity of the respective mass shells, ǫ ∼ ω i . On the other hand, the off-diagonals live on the intermediate shell(s) Given the arguments p and ǫ one can construct the usual relativistic velocity v = p/ǫ entering the Liouville term ∂ t + v∂ x . The velocity does not carry any generation indices because it is determined by the energy spectrum of the two-point correlator. Since the diagonals peak on the mass shells, their propagation velocity matches that of the respective mass eigenstates. Because the off-diagonals live on the intermediate shell, their propagation velocity matches that of the center of momentum. As this is in no way obvious if one deals with the Wigner function or the density matrix, the question of the propagation velocity in the Liouville operator was debated in the literature recently [7,23].

JHEP01(2020)138
Evolution equation for neutrino two-point correlator. The two-point correlator corresponding to equation (2.2) is derived by substituting the latter into equation (1.3) and integrating over τ and y, Using the identity as well as relations k + q = 2p and ω i (k) + ω j (q) = 2ǫ implied by the delta-functions, we can recast equation (2.7) as a product of a shape and a phase factor, The shape factor reads with the velocity v defined as v ≡ p/ǫ. Because the shape factor is a function of vt − x, it satisfies the Liouville equation, Action of the Liouville operator on the two-point correlator, equation (2.8), therefore results in Let us emphasize that the derivation of equations (2.8) and (2.11) does not rely on any approximations, and that the form of equation (2.11) is independent of the initial conditions.
Propagation velocity. Results similar to equation (2.11) have been obtained in a number of previous works [11,12,24,25]. However, the authors quickly turned to the ultrarelativistic approximation setting ǫ = |p| and leaving the meaning of ǫ unspecified.
To clarify the role of ǫ we consider the plane wave limit, ψ i (0, p) ∝ δ(p − p w ). In this limit the shape factor simplifies to

JHEP01(2020)138
where we have omitted the overall normalization for brevity. The second delta-function in equation (2.12) forces the diagonals to the mass shell of the respective mass eigenstate, ǫ = ω i (p w ). On the other hand, the off-diagonals live on the intermediate shell, ǫ = 1 2 (ω i (p w ) + ω j (p w )), see references [26,27] for a related analysis in the context of leptogenesis.
The non-trivial shell structure of the two-point correlator explains why the velocity v in the Liouville operator, see equation (2.11), does not carry any generation indices. The velocity is determined by the value of ǫ on the respective shell and is different for the diagonal and off-diagonal components of ̺ ij (t, x, ǫ, p). In the plane wave approximation the diagonals propagate with the velocity of the respective mass eigenstate, p w /ω i (p w ), whereas the off-diagonals propagate with the velocity of the center of momentum, 2p w /(ω i (p w ) + ω j (p w )).
Derivation from the Schrödinger equation. Because the dynamics of the underlying wave functions is governed by the Schrödinger equation, one can expect that, similarly to the Liouville-von Neumann equation for the neutrino density matrix, the evolution equation for the neutrino two-point correlator (2.11) can be derived from the Schrödinger equation without using an explicit solution of the latter. To demonstrate this we first recast equation (2.11) in the form Next, we differentiate equation (1.3) with respect to the central time t, multiply the result of differentiation by ǫ, and further use ǫe iǫτ = −i∂ τ e iǫτ on the right-hand side. Integrating by parts we can replace −i∂ τ e iǫτ by ie iǫτ ∂ τ . Using the chain rule we finally arrive at (2.14) Proceeding similarly, i.e. differentiating equation (1.3) with respect to the central coordinate x, multiplying the result of differentiation by p, using e −ipy p = i∂ y e −ipy , integrating by parts to replace i∂ y e −ipy by −ie −ipy ∂ y , and finally using the chain rule we obtain Combining equations (2.14) and (2.15) we arrive at x). Action of the d'Alembert operator hence results in ψ i = −m 2 i ψ i , as expected. Comparing the right-hand side of equation (2.16) to the definition of the two-point correlator we conclude that the former x, ǫ, p). In other words, we recover equation (2.13).
Quantum-kinetic equation. Equation (2.13) can be rewritten in a manifestly covariant form, where x µ = (t, x) and p µ = (ǫ, p), see appendix A for more details. Its form matches the form of the quantum kinetic equation for the left-handed neutrino component (considered in the vacuum limit) discussed in reference [11]. In the quantum kinetic approach to oscillations the two-point correlator ̺(x, p) is replaced by the Wightman propagator ̺ < (x, p), defined as the expectation value of the neutrino field operators [11][12][13][14][15][16]. A few further details are presented in appendix E. Apart from the quantum kinetic equation, the Wightman propagator must also satisfy the complementary constraint equation that determines its spectrum.
In addition to the neutrino and antineutrino correlators, the formalism of the nonequilibrium quantum field theory also naturally incorporates the particle-antiparticle correlators [28][29][30][31]. The authors of references [32][33][34][35] have shown that the particle-antiparticle correlations also affect neutrino oscillations in vacuum, albeit their contribution is strongly suppressed in the physically interesting relativistic limit. As has been argued in reference [36], they are nevertheless important to ensure that neutrino detection probability exactly vanishes outside of the light cone. An approach, in which faster-than-light signaling is automatically and manifestly prevented, has been developed in reference [37]. This approach works directly with probabilities instead of squared matrix elements. The probability-level calculation implicitly sums over the unobserved emissions that are crucial to ensuring that the measurement is local and consistent with causality. The need to sum inclusively over unobserved emissions was first described in reference [38] and made explicit in reference [39].

Description in terms of the Wigner function
The Wigner function ̺ ij (t, x, p) is obtained from the two-point correlator by integration over ǫ. Since the Liouville term also depends on ǫ, strictly speaking, the integration does not yield a closed-form evolution equation for the Wigner function. On the other hand, for typical initial conditions the width of the two-point correlator in the ǫ space is negligibly small and the off-shell contributions can be neglected. In the relativistic limit this results in an evolution equation for the Wigner function whose form matches the form of (the vacuum limit of) the kinetic equation of Sigl and Raffelt.
The latter is believed to be inconsistent with the quantum mechanical uncertainty principle. We argue below that this is not the case. Applied to the Wigner function, the uncertainty principle states that the stronger the Wigner function is peaked around JHEP01(2020)138 a characteristic momentum p w , the weaker is its dependence on vt − x, and vice versa. If the Wigner function initially satisfies the uncertainty principle, the evolution equation ensures that this holds true in the course of the neutrino propagation. However, unlike the Schrödinger equation, the evolution equation itself does not enforce the inverse relation between localization in the coordinate and momentum spaces.
Each individual momentum mode of the Wigner function does not experience any suppression in the course of the neutrino propagation. On the other hand, because of nonzero width of the Wigner function in the momentum space related to the uncertainty principle, off-diagonal components of the density matrix do get suppressed at late times due to dephasing of the momentum modes.
Evolution equation for Wigner function. In the plane wave limit the Liouville equation (2.11) can be trivially integrated over ǫ. On the left-hand side the integration yields a matrix of velocities, .
On the right-hand side the integration yields the familiar Let us emphasize that equation (2.19) is a direct consequence of the fact that the offdiagonals of the two-point correlator live on the intermediate shell [27]. The resulting generalization of the classical Liouville equation to the case of mixing neutrinos reads Also for initial conditions other than the plane wave ones, the two-point correlator is typically strongly peaked in the vicinity of the respective shells. In this case the small variations of ǫ around the on-shell value can be safely neglected and integration of equation (2.11) again yields equation (2.20).
The uncertainty principle. In the plane wave approximation the neutrino is completely localized in the momentum space and completely delocalized in the coordinate space. Let us now return to the case of Gaussian initial conditions, see equation (2.3).
Trading the integration variables k and q in equation (2.9) for the center and relative momentum, s ≡ (k + q)/2 and ∆ ≡ k − q, and integrating over the center momentum we obtain for the shape factor

JHEP01(2020)138
Because contributions of ∆ larger than the neutrino momentum uncertainty are strongly suppressed, for σ p ≪ p we can safely approximate ω i (p + ∆/2) + ω j (p − ∆/2) by ω i (p) + ω j (p), thereby approximating the propagation velocity by v ij (p). The integration over ǫ is then trivial. The resulting approximate Wigner function also factorizes into a shape and a phase factor, with the shape factor given by As can be verified by substitution, equation (2.22) satisfies the evolution equation (2.20) for arbitrary initial conditions. For Gaussian initial conditions the remaining integration over ∆ in equation (2.23) can be done analytically and yields The form of the shape factor for Lorentzian initial conditions is presented in appendix B. Equation (2.24) implies that the evolution equation admits solutions consistent with the uncertainty principle. Note, however, that equation (2.22) remains a solution of the evolution equation even if σ p and σ x are considered as independent. In particular, one can simultaneously take the limits σ p → 0 and σ x → 0, thereby 'creating' a state that is localized both in the momentum and coordinate spaces. That is, the evolution equation itself does not enforce the inverse relation between localization in the coordinate and momentum spaces. The latter is encoded in the initial conditions for ̺ ij (t, x, p) instead.
Kinematical decoherence. In the approximation used to derive equation (2.22), individual momentum modes of the Wigner function do not experience any suppression in the course of the neutrino propagation. On the other hand, their growing dephasing results in kinematical decoherence and related suppression of the off-diagonals of the density matrix, known as wave packet separation [6,23].
To demonstrate this we derive the respective density matrix by integrating equation (2.22) over p. For a wave function strongly peaked around the characteristic momentum p w the Wigner function can be expanded in the vicinity of p w . In the relativistic limit v ij (p) ∼ 1 and is not sensitive to small, of the order of ∼ σ p , variations of p around p w . Therefore, we can approximate v ij (p) by v ij (p w ) in the shape factor, see equation (2.23). On the other hand, expanding the argument of the phase factor, , we find that it is sensitive to the difference of the respective group velocities. In this approximation the Wigner function reads

JHEP01(2020)138
Integrating over p we recover equation (2.6), Therefore, kinematical decoherence in the density matrix is caused by dephasing of the individual momentum modes, which grows in the course of the neutrino propagation. This picture is equivalent to the one of the wave packet separation [6]. The smaller is the considered momentum range, the less pronounced is the suppression. This corresponds to coherence restoration in a detector with a very good energy resolution. Let us emphasize, that the finite energy resolution of a detector and an emission process can just as well be taken into account by post-processing the results through an integral over energy, instead of considering solutions of finite width in the momentum space. The nonzero width of the wave packets in the momentum space, crucial for the onset of kinematical decoherence at late times, can be traced back to the uncertainty principle and a nonzero spatial size of the wave packets [8]. It is interesting to note, however, that if we were to treat σ p and σ x as independent and set the latter to zero in equation (2.23), we would still get the same suppression factor in equation (2.26). This (inconsistent from the viewpoint of quantum mechanics) approximation still solves equation (2.20), as has been discussed above. If we considered an ensemble of neutrinos, this approximation would correspond to describing neutrinos as classical particles with a particular momentum distribution. This picture has been used in reference [40] to study collective oscillations of a two-flavor neutrino system. Hence, the results obtained there do take into account the effect of wave packet separation, at least for low neutrino densities. For recent developments in the field of collective neutrino oscillations see e.g. references [41][42][43][44][45][46][47][48][49].
As has been noted in references [7,8], because ̺ ij (t, x, p w ) satisfies the evolution equation (2.20) with p = p w the density matrix also satisfies an effective Liouville equation, where C ij (t, x) stems from the last term on the right-hand side of equation (2.26) and describes kinematical decoherence. The explicit form of C(t, x) depends on the initial conditions and can be calculated a posteriori, once the solution for the density matrix is known [7].
Spread of the wave packet. Because different momentum modes propagate with different velocities, the width of the density matrix increases with time. To estimate this effect we that the expansion is understood to be made in the direction of the neutrino propagation). Integrating over p we recover equation (2.6) with σ 2 x replaced by an effective generation-and time-dependent width [50], The phase factor also receives a tiny, of the order of ∆v ij (p w )v ′ ij (p w ), correction that is proportional to v ij (p w )t − x and hence vanishes at the center of the wave packet.

JHEP01(2020)138
Off-shell effects. As can be read off from equation (2.21), for each p there is a range of allowed values of ǫ. In other words, the two-point correlator possesses an off-shell component. Therefore, even for a fixed value of the momentum, ̺ ij (t, x, p) is a superposition of plane waves propagating with different velocities v = p/ǫ centered around v ij with a tiny spread ∼ v ij (p)∆v ij (p) · σ p /(ω i (p) + ω j (p)). By analogy with the wave packet separation one could expect kinematical decoherence in the Wigner function caused by a gradual dephasing of the individual modes at late times. Analysis of this subtle effect is beyond the scope of the present work.
Kinetic equation. In the relativistic limit v ij is well approximated by the average velocityv ij introduced earlier. In this approximation the evolution equation (2.20) can be recast in the form that matches the form of the kinetic equation discussed in references [9,23], where [. , .] denotes a commutator, and {. , .} an anticommutator. In the kinetic equation ̺(t, x, p) is replaced by the matrix of densities, defined as the expectation value of bilinears of the creation and annihilation operators of the neutrino field [9].

Neutrino propagation in matter
In this section we generalize the results of section 2 to neutrino propagation in an external potential. In subsection 3.1 we derive, to the first order in the gradient expansion, an evolution equation for the neutrino two-point correlator, valid also for non-relativistic neutrinos. Up to small corrections related to a somewhat different expansion scheme its form matches (the collisionless limit of) the quantum kinetic equation. Analyzing the energy dependence of its terms we derive a simplified evolution equation valid in the relativistic limit.
In subsection 3.2 we use the latter to derive an evolution equation for the Wigner function. Its form matches (the collisionless limit of) the kinetic equation of Sigl and Raffelt. Considering Gaussian initial conditions we show that it admits solutions consistent with the uncertainty principle and accounts for neutrino wave packet separation. We limit our analysis to forward scattering and do not consider the Lindblad terms [51][52][53] describing dynamical decoherence.

Description in terms of the two-point correlator
Because the Liouville operator on the left-hand side of the evolution equation for the two-point correlator, ǫ∂ t + p∂ x , has the dimension of energy squared, its right-hand side contains (anti) commutators of the two-point correlator with the Hamiltonian squared. For non-commuting kinetic and potential energy operators this yields a new contribution proportional to their commutator. Shifts of the on-shell value of ǫ induced by a timedependent potential also result in a new contribution proportional to the time-derivative of the potential. In the relativistic limit, which we consider first to make contact with results of references [9,23], these new terms can be neglected and the evolution equation reverts to the conventional one.
In order to keep the derivation as simple as possible we assume that the background matter is unpolarized and also free of convective currents. With these assumption the matter is characterized by a potential V ij (t, x) determined by occupation numbers of the charged leptons. Because equation contributions from the kinetic and potential terms can be treated separately.
The contribution of the potential term can be conveniently analyzed in the coordinate representation employed in equations (1.3) and (3.1). Expanding V ij (t ± τ /2, x ± y/2) to the first order in τ and y, using e iǫτ τ = −i∂ ǫ e iǫτ and e −ipy y = i∂ p e −ipy in equation (3.1), and comparing the resulting expression to the definition of the two-point correlator, see equation (1.3), we find for contribution of the potential term Contribution of the kinetic term is more easily analyzed in the momentum representation where the action of the kinetic operator amounts to multiplication of the wave functions by ω, To leading order in the gradients ω i (p + ∆/2) − ω j (p − ∆/2) ≈v ij (p)∆ + ω i (p) − ω j (p) in the relativistic limit. Using e i∆x ∆ = −i∂ x e i∆x and comparing the resulting expression to the definition of the two-point correlator in the momentum representation, see appendix C, we find that in the relativistic limit contribution of the kinetic term can be approximated by

5)
JHEP01(2020)138 valid in the relativistic limit. The first two terms on the left-hand side generalize the flux conservation equation. The third one describes momentum changes by coherent external forces [23]. Finally, the fourth one describes spectrum shifts caused by a changing external potential.
As we have seen in section 2, in order to derive a closed-form evolution equation without resorting to the relativistic limit one has to work with the Liouville operator of the form ǫ∂ t + p∂ x . The form of the resulting evolution equation in matter is discussed in the remainder of this subsection.
Neutrino propagation in a constant potential. Before generalizing equation (2.17) to space-or time-dependent potentials let us first consider the technically simple case of a constant potential V ij . By analogy with equation (2.17) we multiply equation (3.1) by ǫ and further use ǫe iǫτ = −i∂ τ e iǫτ on the right-hand side. Integrating by parts we can replace −i∂ τ e iǫτ by ie iǫτ ∂ τ . Because the Hamiltonian H ij (x) = K ij (x) + V ij is time-independent, subsequently using the chain rule and the Schrödinger equation (2.1) we arrive at where we employ a compact notation H 2 ij ≡ H in H nj . In the momentum representation action of the Hamiltonian H ij (x) = K ij (x) + V ij amounts to multiplication of the wave function by H ij (p) = δ ij ω i (p) + V ij , see appendix C for more details. Expressing the coordinate-representation wave functions in equation (3.6) in terms of their momentumrepresentation counterparts and subsequently introducing the center and relative momenta we arrive at To the first order in the gradient expansion H 2 In vacuum H 2 ij (p) = δ ij ω 2 j (p), ∂ p H 2 ij (p) = p and equation (3.8) reverts to equation (2.17). In appendix C equation (3.8) is re-derived using the solution of the Schrödinger equation in a constant potential. As can be read off from the explicit expression for the two-point correlator, its spectrum is determined by eigenvalues of the Hamiltonian and is shifted with respect to the vacuum spectrum by terms proportional to the potential. Thus, the evolution equation (3.8) must be supplemented by initial conditions that consistently account JHEP01(2020)138 for modifications of the neutrino spectrum in an external potential. As has been mentioned above, in the quantum kinetic approach to neutrino oscillations the spectrum can be obtained from the supplementary constraint equation.
Neutrino propagation in a spatially inhomogeneous potential. For a spatially inhomogeneous but time-independent potential equation (3.6) still holds. However, because the operators K ij (x) and V ij (x) no longer commute, it does not reduce to equation (3.7) in the momentum representation. To compute the right-hand side of equation (3.6) we reorder the operators in H 2 such that V is left to K, plus a commutator, (3.9) Expanding the kinetic energy operator in powers of ∂ x /m we find to the first order in the gradients where a ℓ are coefficients of the Taylor expansion. Substituting the first line of equation (3.9) into equation (3.6) and trading the momentum integration variables for the center and relative momenta we obtain where H ij (x, p) = δ ij ω j (p) + V ij (x). Next, we expand the square of the Hamiltonian to the first order in the gradients, H 2 (x ± y/2, s ± ∆/2) ≈ H 2 (x, s) + ∂ s H 2 (x, s) ∆/2 + ∂ x H 2 (x, s) y/2. The first term of the expansion yields the commutator which generalizes the right-hand side of equation (3.8) to spatially inhomogeneous potentials. Using e i∆x ∆ = −i∂ x e i∆x and the product rule we obtain for contribution of the second expansion term which generalizes the second term on the left-hand side of equation (3.8) to spatially inhomogeneous potentials. Using e −ipy y = i∂ p e −ipy we find for contribution of the third JHEP01(2020)138 (3.14) Substituting the second line of equation (3.9) into equation (3.6), using the approximation equation (3.10), and trading the momentum integration variables for the center and relative momenta we find: Because the right-hand side of equation (3.15) already contains a derivative of the potential, to the first order in the gradient expansion it is sufficient to approximate . In this approximation we obtain whose structure is similar to that of the second term on the right-hand side of equation (3.14). Combining equations (3.12), (3.13), (3.14) and (3.16) we finally arrive at Equation (3.17) implies that in a spatially inhomogeneous potential neutrino oscillations are governed by H 2 (x, p) As this term is a sum of the Hamiltonian and the commutator of the kinetic and potential energy operators, its structure replicates the structure of equation (3.9).
Neutrino propagation in a time-dependent potential. To generalize equation (3.6) to homogeneous but time-dependent potentials we again multiply equation (3.1) by ǫ and use ǫe iǫτ = −i∂ τ e iǫτ on the right-hand side, subsequently replacing −i∂ τ e iǫτ by ie iǫτ ∂ τ upon integration by parts. Because the Hamiltonian H ij (t, x) = K ij (x) + V ij (t) is timedependent, using the chain rule and the Schrödinger equation (2.1) we find that in addition to an expression similar to the right-hand side of equation (3.6), Written in this form, the evolution equation for the two-point correlator is valid for generic potentials. For a time-dependent but homogeneous potential the kinetic and potential energy operators commute and the action of the Hamiltonian amounts, in the momentum representation, to multiplication of the wave function by H ij (t, p) = δ ij ω i (p) + V ij (t). Expressing the coordinate-representation wave functions in equations (3.18) and (3.19) in terms of their momentum-representation counterparts and subsequently introducing the center and relative momenta we obtain respectively. Proceeding as above we expand the Hamiltonian to the first order in the gradients in equation (3.20), ij (t, s)τ /2. Subsequently using e i∆x ∆ = −i∂ x e i∆x and e iǫτ τ = −i∂ ǫ e iǫτ we obtain from equation (3.20) Because equation (3.21) already contains a derivative of the potential, to the first order in the gradient expansion it is sufficient to approximate ∂ t V in (t ± τ /2) by ∂ t V in (t). In this approximation we obtain from equation (3.21) Combining equations (3.22) and (3.23) we finally arrive at

JHEP01(2020)138
In appendix D we rederive the left-hand side of equation (3.24) for a single neutrino generation using the gradient expansion to obtain an explicit expression for the two-point correlator. This expression illustrates that a homogeneous but time-dependent potential affects primarily the spectrum of the two-point correlator. Furthermore, the terms x, ǫ, x)} have a common origin and stem from time-derivative of the energy-conserving delta-function. Therefore, even though at first sight the {∂ t V(t), ̺(t, x, ǫ, x)} term seems to describe effects of absorption or loss of coherence, it actually describes shifts of the spectrum induced by the time-dependent potential, just like the term {∂ t H 2 (t, p), ∂ ǫ ̺(t, x, ǫ, p)}.
Recovering relativistic limit. Since we work to linear order in the gradient expansion here, the quantum kinetic equation for neutrinos propagating in a general potential can be deduced by combining equations (3.17) and (3.24), To conclude our analysis of the evolution equation for the two-point correlator we demonstrate how to recover equation (3.5), derived in the relativistic limit, from equation (3.25), which is valid for any value of the momentum.
To handle the first anticommutator on the left-hand side of equation ( To the first order in the potential the second term can be approximated by [∂ x V, [ω, ∂ p ̺]] which vanishes in the limit of ω ∝ 1. Hence this term is proportional to ∆ω and is again negligibly small compared to the first one. The first term can be approximated by {ω, {∂ t H, ∂ ǫ ̺}} ij = (ω i + ω j ) · {∂ t H, ∂ ǫ ̺} ij , and the second, to the first order in the potential, by The second term is proportional to ∆ω and, in the relativistic limit, negligibly small compared to the first.
On the right-hand side of equation ( Common to all the (leading) terms discussed so far is the overall factor ω i + ω j . The two new terms, {∂ t V(t, x), ̺(t, x, ǫ, x)} one the left-hand side and [∂ x V(t, x), ∂ p ω(p)] on the right-hand side, do not share this feature. Therefore, they can be neglected in the relativistic limit.

JHEP01(2020)138
As has been discussed in section 2, barring small off-shell contributions 2ǫ ≈ ω i + ω j . Dividing the left-and right-hand side of equation (3.25) by ω i + ω j we therefore recover equation (3.5). It is also worth mentioning that the subleading contribution ∆v ij (ω i − ω j )/(ω i + ω j )/2 that stems from the anticommutator {∂ p H 2 (t, x, p), ∂ x ̺(t, x, ǫ, p)} is equal to v ij −v ij . This explains why the center of momentum velocity is replaced by the average velocity in the relativistic approximation.
Quantum-kinetic equation. Depending on its nature the neutrino possesses either two (Majorana particle) or four (Dirac particle) components. On the other hand, the wave function in the quantum-mechanical approach to neutrino oscillation describes a single component, either a left-handed neutrino or a right-handed antineutrino. For typical experimental setups this is well motivated because the spin-flipping (or active-sterile) conversions are strongly suppressed by the smallness of the m/ω ratio. Nevertheless, strictly speaking, this means that the quantum-field theoretical counterpart of the conventional quantum-mechanical approach to neutrino oscillations is the theory of mixing scalar fields. To the first order in the gradient expansion the quantum kinetic equation for the respective Wightman propagator reads with the effective Hamiltonian defined as see appendix E for more details. The local self energy Σ(t, x) models the lowest-order contribution to the neutrino self-energy evaluated in the limit m W,Z → ∞ [12]. The form of equation (3.26) is almost identical to the form of equation (3.25), except for the absence of the two new terms, {∂ t V(t, x), ̺(t, x, ǫ, x)} and [[∂ x V(t, x), ∂ p ω(p)], ̺(t, x, ǫ, p)]. Their absence can be traced back to the slightly different gradient expansion schemes used to derive equations (3.25) and (3.26). Whereas to derive the former we expanded the potential V, to derive the latter we expanded the self-energy Σ. The self-energy can be approximately expressed through V as Σ ∼ ωV + V ω. Therefore, appearance of terms containing a different combination of ω and V in the quantum kinetic equation is not possible. As has been argued above, in the relativistic limit the two new terms are subdominant and can be safely neglected anyway.

Description in terms of the Wigner function
The (relativistic) evolution equation for the two-point correlator (3.5) can be trivially integrated over ǫ. The resulting evolution equation for the Wigner function matches (the collisionless limit of) the kinetic equation of Sigl and Raffelt and admits solutions consistent with the uncertainty principle also for neutrinos propagating in matter.

JHEP01(2020)138
Kinetic equation. To derive evolution equation for the Wigner function we integrate equation (3.5) over ǫ. Because ̺(t, x, ǫ, p) vanishes for ǫ → ±∞, the last term on the left-hand side vanishes upon the integration. Integration of the remaining terms yields valid in the physically relevant relativistic limit. The form of equation (3.28) matches the (collisionless limit of the) kinetic equation of Sigl and Raffelt [9]. As has been emphasized in reference [23], in the relativistic limit the first anticommutator on the left-hand side of equation (3.28) contains the full flavor-dependent velocity structure of the Liouville term.
Neutrino propagation in a constant potential. For a constant potential the second anticommutator on the right-hand side of equation (3.28) vanishes and it simplifies to We can construct its solution by substituting the (well known) solution of the Schrödinger equation in a constant potential, ψ(t, p) = e −iH(p)t ψ(0, p), into equation (1.2). This results in

JHEP01(2020)138
with the effective velocity υ(t, p) given by Performing a gradient expansion of the Hamiltonian in equation (3.30) and subsequently using the Feynman disentanglement theorem we obtain for the shape factor with T c and T a being the time-ordering and the anti-time-ordering operators respectively. As can be verified by substitution, equation (3.34) is consistent with equations (3.32) and (3.33).
The uncertainty principle. The matrices H(p) and ∂ p H(p) do not commute in general. The lowest order approximation to the effective velocity consists in neglecting [H(p), ∂ p H(p)] in equation (3.33). Parametrically this commutator is of the order of sin 2θV∆v ij , where θ is the vacuum mixing angle, and is therefore strongly suppressed with respect to the leading term ∂ p H(p). It becomes comparable to ∂ p H(p) only at late times, sin 2θV|∆v ij |t ∼ 1.
For V ∼ 1 eV, ω ∼ 1 MeV, and the solar mass difference this amounts to a distance that is an order of magnitude larger than the solar radius. This approximation results in and shape factor of the form For Gaussian initial conditions the approximate shape factor equation (3.35) takes the form identical to equation (2.24), up to the (negligible in the relativistic limit) replacement v ij →v ij , As has been argued in section 2.3, this shape factor is consistent with the uncertainty principle. The same applies to Lorentzian initial conditions, see appendix B. In other words, also for neutrino propagating in matter the evolution equation admits solutions consistent with the uncertainty principle.
Kinematical decoherence. As follows from equations (3.31) and (3.36), also in matter each individual momentum mode of the Wigner function does not experience any suppression in the course of collisionless neutrino propagation. On the other hand, their growing dephasing results in kinematical decoherence in the density matrix at late times.
To demonstrate this we derive the latter by integrating equation (3.31) over p. For Gaussian initial conditions the integral can be estimated using the saddle point approximation.

JHEP01(2020)138
To this end we use e ∓iH(p)t = U(p)e ∓iE(p)t U † (p), where U(p) is a unitary matrix diagonalizing the Hamiltonian, U † (p)H(p)U(p) = E(p). Expanding E(p) in the vicinity of p w , E(p) ≈ E(p w ) + u(p w )(p − p w ), where u(p) ≡ ∂ p E(p) denotes group velocity of the propagation eigenstates, and taking into account that E(p) and u(p) commute, we find for the density matrix x)e iH(pw)t , (3.37) with the shape factor given (for the Gaussian initial conditions) by Hence, for ∆u = 0 oscillations are suppressed at late times. In vacuum equation (3.37) reverts to equation (2.26). Let us emphasize that even though the evolution equation (3.29) neglects medium corrections to the neutrino propagation velocity, it does take into account medium corrections to the group velocity of the propagation eigenstates. Close to the Mikheyev-Smirnov-Wolfenstein resonance [54] the latter, though small in the absolute value, are of the same order of magnitude as the difference of the vacuum neutrino velocities. As a result ∆u mn can vanish. Therefore, matter effects can suppress development of the kinematical decoherence in the density matrix under certain circumstances, see e.g. [50,54] and references therein. One should however keep in mind, that scattering on the matter particles induces dynamical decoherence in the Wigner function and density matrix [16].
Neutrino propagation in a time-dependent potential. As the width of the neutrino wave packet in coordinate space is typically 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 time-dependent potential, which is governed by The value of this effective matter potential at a time t is given by the potential at the position the neutrino wave packet reaches by that time. This approximation is widely used in the literature, see e.g. references [8,55].
Because the Hamiltonian at different times does not in general commute with itself, the solution of the Schrödinger equation contains the time-ordered exponent, ψ(t, p) = T c e −i t 0 H(t, p)dt ψ(0, p). Therefore, the evolution operators e −iH(p+∆/2)t and e iH(p−∆/2)t in equation (3.30) are replaced by T c e −i t 0 H(t, p+∆/2)dt and T a e i t 0 H(t, p−∆/2)dt ,

JHEP01(2020)138
Equation (3.40) only approximately solves equation (3.39). Differentiating the former with respect to t we obtain −iH(t, p+∆/2)T c e −i t 0 H(t, p+∆/2)dt and T a e i t 0 H(t, p−∆/2)dt iH(t, p− ∆/2) under the integral sign. Similar to the case of a constant potential, to recover equation (3.39) we need to perform a gradient expansion of the pre-exponential factors, H(t, p ± ∆/2) ≈ H(t, p) ± ∂ p H(t, p)∆/2 (and subsequently use ∆e i∆x = −i∂ x e i∆x ). This implies that the solution of equation (3.39) is obtained from equation (3.40) by performing a gradient expansion of the Hamiltonian in the (anti-) time-ordered exponents.
Proceeding similarly to the case of a constant potential we tentatively write the Wigner function in the form From equations (3.39) and (3.41) it follows that the shape factor satisfies equation (3.32) with the effective velocity now given by Performing a gradient expansion of the Hamiltonian in equation (3.40) and further using the Feynman disentanglement theorem we obtain for the shape factor an expression identical to equation (3.34) but with the effective velocity given by equation (3.42).

Summary and conclusion
The present work establishes a connection between three approaches to the description of neutrino oscillations: (i) a quantum-mechanical approach operating with the neutrino wave function or density matrix, (ii) a quantum kinetic approach operating with the Wightman propagator, and (iii) a kinetic approach operating with the matrix of densities. The quantum-mechanical approach is well suited to analysis of neutrino oscillations in vacuum and its collisionless propagation in matter. In particular, it is automatically consistent with the uncertainty principle and, as a result, takes into account the effect of wave packet separation. On the other hand, it is not capable of accounting for scattering processes. This limits its applicability to analysis of neutrino propagation in exploding supernovae, where oscillations and collisions may be equally important.
The quantum kinetic equation for the Wightmann propagator is derived using methods of non-equilibrium quantum field theory. An evolution equation for its quantum-mechanical counterpart -the two-point correlator -can be derived from the Schrödinger equation. For neutrino propagation in vacuum these two (quite distinct) approaches lead to identical equations. For collisionless neutrino propagation in matter the resulting equations differ only by subleading terms which can be traced back to slightly different gradient expansion schemes. Analysis of the two-point correlator for typical initial conditions reveals that in addition to the usual mass shells ǫ = ω i there are ǫ = (ω i + ω j )/2 shells responsible for neutrino oscillations. Their existence explains why the Liouville term of the evolution (as well as of the quantum kinetic) equation, ∂ t + v∂ x , does not carry any generation indices. The propagation velocity v = p/ǫ is determined by the spectrum of the twopoint correlator. In the propagation basis its diagonal elements peak in the vicinity of the JHEP01(2020)138 respective mass shells and propagate with velocity v ≈ p/ω i . The off diagonal elements peak in the vicinity of the oscillation shells and propagate with the center of momentum velocity v ≈ 2p/(ω i + ω j ).
For typical initial conditions the off-shell contributions are negligible and it is sufficient to keep track of the evolution of the two-point correlator in the vicinity of the mass and oscillation shells. In this approximation the evolution equation for the two-point correlator can be reduced to a closed-form evolution equation for the Wigner function. In its Liouville term the velocity v is replaced by the matrix of velocities, v ij = 2p/(ω i + ω j ). In the physically relevant relativistic limit v ij is well approximated by the average velocityv ij = (v i + v j )/2. In this approximation the form of the evolution equation for the Wigner function is identical to the form of the (collisionless limit of the) kinetic equation for the matrix of densities. The latter is derived using methods of perturbative quantum field theory. Solutions of the evolution equation for the Wigner function, derived from the Schrödinger equation, can be obtained from known solutions of the latter. The resulting Wigner function is localized neither in the momentum, nor in the coordinate spaces. The stronger it peaks around a characteristic momentum p w , the weaker is its dependence onvt − x, and vice versa. In other words, the evolution equation for the Wigner function admits solutions consistent with the uncertainty principle. Let us emphasize that unlike the Schrödinger equation, it does not enforce the uncertainty principle and admits also classical solutions describing a particle with definite coordinate and momentum. On the other hand, it consistently evolves solutions satisfying the uncertainty principle if supplemented by respective initial conditions. For a neutrino propagating in vacuum or collisionless matter each individual momentum mode of the Wigner function does not experience any suppression in the course of the neutrino propagation. On the other hand, their growing dephasing results in kinematical decoherence in the density matrix at late times.
One can expect that the three approaches to neutrino oscillations produce identical results for neutrino propagation in vacuum or an external potential (up to small differences related to different approximation schemes). Results of the present work indicate that this is indeed the case. Because the evolution equations for the two-point correlator and Wigner function, derived from the Schrödinger equation, match the quantum kinetic and kinetic equation respectively, the two-point correlator and Wigner function, derived from the solution of the Schrödinger equation, solve the quantum kinetic and kinetic equation respectively. This observation hints, that the effect of wave packet separation can be accounted for also in the (quantum) kinetic approach to neutrino oscillations by considering initial conditions consistent with the uncertainty principle. Given that the (quantum) kinetic approach naturally incorporates scattering processes, this finding speaks in favor of using the (quantum) kinetic description for the analysis of neutrino propagation in exploding supernovae where neutrino oscillations and collisions, as well as the effect of wave packet separation, might be equally important.

JHEP01(2020)138
As can be inferred from equation (B.2), localization in the momentum space is inversely proportional to localization in the coordinate space, in agreement with the uncertainty principle.

C Propagation in a constant potential
In the momentum representation equation (1.3) can be recast in the form Fourier-transforming equation (2.1) we obtain the momentum-representation Schrödinger equation The momentum-representation Hamiltonian is related to its coordinate-representation counterpart by The Hamiltonian is diagonalized by a unitary transformation, U † in (p)H nm (p)U mj (p) = δ ij Ω i (p). The wave function in the propagation basis, Ψ i (t, p) = U † ij (t, p)ψ j (t, p), satisfies the Schrödinger equation i∂ t Ψ i (t, p) = Ω i (p)Ψ i (t, p). Its solution reads Ψ i (t, p) = e −iΩ i (p)t Ψ i (0, p). Transforming the solution back to the neutrino mass basis we obtain ψ i (t, p) = U in (p)e −iΩn(p)t Ψ n (0, p) . (C.5) Substituting equation (C.5) into equation (C.1) and integrating over τ we find for the two-point correlator As expected, the spectrum of the two-point correlator, 2ǫ = Ω n (p + ∆/2) + Ω m (p − ∆/2), is determined by eigenvalues of the Hamiltonian. Proceeding as in section 2 we represent the phase difference as . with H(t, p) = ω(p) + V(t). Its solution is given by [7] ψ(t, p) = e −i t 0 dt ′ H(t ′ , p) ψ(0, p) .

E Quantum kinetic equation
In the collisionless approximation (i.e. including only forward scattering) the quantum kinetic equation for the Wightman propagators, ̺ ≶ (u, v), of a system of mixing scalar fields reads [56,57]: Further expanding the self-energy to the first-order in the gradients we obtain: Note that following the standard convention we use the same symbol for the Wightman propagator considered as a function of the center and relative coordinates. Expressing ̺ < (x, s) in terms of its Wigner-transform, where x = (t, x) and p = (ǫ, p). The Wightman propagator ̺ < (t, x, ǫ, p) is the quantum kinetic counterpart of the quantum-mechanical two-point correlator ̺(t, x, ǫ, p).
To rewrite equation (E.5) in a form analogous to equation (3.25) we separate the time and space components of ∂ x . Introducing an effective Hamiltonian, H 2 (t, x, p) ≡ p 2 + m 2 + Σ(t, x) = ω 2 (p) + Σ(t, x) , (E.6) and using the identity 1 4 {∂ p H 2 (t, x, p), ∂ x ̺ < (t, x, ǫ, p)} = p∂ x ̺ < (t, x, ǫ, p) we can finally recast equation (E.5) in the form: Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.