The Liouville equation for flavour evolution of neutrinos and neutrino wave packets

We consider several aspects related to the form, derivation and applications of the Liouville equation (LE) for flavour evolution of neutrinos. To take into account the quantum nature of neutrinos we derive the evolution equation for the matrix of densities using wave packets instead of Wigner functions. The obtained equation differs from the standard LE by an additional term which is proportional to the difference of group velocities. We show that this term describes loss of the propagation coherence in the system. In absence of inelastic collisions, the LE can be reduced to a single derivative equation over a trajectory coordinate. Additional time and spacial dependence may steam from initial (production) conditions. The transition from single neutrino evolution to the evolution of a neutrino gas is considered.


Introduction
A field of propagating neutrinos inside a supernova can be described by the matrix of densities in the flavour space, n(t, x, v, E), which depends on the space-time coordinates (x, t), velocity v and energy E. The evolution of the matrix of densities in space and time is governed by the Liouville equation (see e.g. [1]) where is the Liouville operator 1 and H(x, t, v, E) is the Hamiltonian that depends on properties of the background. The velocity matrix v is usually taken to be proportional to the unit matrix although its flavour dependence has previously been noticed [1,2]. Eq. (1.1) has been derived in a number of papers [1][2][3][4][5][6][7][8] and most of the derivations are based on a presentation of the density matrix in terms of the Wigner function [9]. The latter is a specific Fourier transform of the neutrino field bilinear. Furthermore, most of the derivations rely on formalism of statistical mechanics [10].
There is a number of questions related to the form, derivation and applications of eq. (1.1) for supernova neutrinos. 1) Using the Wigner function for the derivation of the LE from first principles. This function allows to reconcile the classical nature of the Liouville equation obtained from conservation of phase space with the quantum nature of propagating neutrinos. Indeed, in order to derive the Liouville equation, the simultaneous treatment of the neutrino momentum and the neutrino coordinate is needed [10]. Such a treatment of a quantum system can be obtained using the Wigner functions. However, the Wigner functions come with the price that in an inhomogeneous background, they can attain negative values, and this would also propagate to quantities such as the matrix of densities.
2) Origin of two separate derivatives over time and spatial coordinates in the equation.
In some recent studies the derivatives in space and time in the Liouville equation are treated independently [11,12]. One concern regarding these studies is that the locality of the solution could be violated and that information could be transferred superluminally from one space point to another [13].
3) The correspondence of eq. (1.1) to the evolution equation of the individual neutrino. Indeed, in the absence of inelastic collisions we can in principle consider evolution of the individual neutrino moving along a certain trajectory (and not an ensemble). For a single neutrino, two derivatives in eq. (1.1) are reduced to the one over the coordinate along the trajectory. Then a gas of neutrinos is just the system of independently evolved neutrinos, and the density is the sum over neutrinos in a given unit volume. In general, for a gas, x and t are independent since different values of (x, t) correspond to different neutrinos. For a single neutrino x and t are connected by x = vt up to uncertainty principle.
4) The form and nature of the matrix v in eq. (1.1). Notice that in quantum theory the velocity v follows from the solution of the equation rather than appears in the equation itself. However if the inelastic collisions are neglected and one can consider propagation of single neutrinos, then v should have an interpretation in terms of velocities of the individual neutrinos. In the case of two mixed neutrinos we can speak about two group velocities. Using the averaged velocity in the LE, so that v = v 0 1, should be justified. Indeed, the difference of velocities in v produces the term ∆vp ≈ ∆m 2 /2|p| (in vacuum), which is comparable to the term with the Hamiltonian (after subtraction of the part proportional to the unit matrix).
In this paper we will address these four groups of questions. Instead of Wigner functions we propose to use wave packets which explicitly encode the uncertainty relation, and so, allow to reconcile the momentum and coordinate description. Our method is the following: We reconstruct the solutions of the evolution equation for the wave function using physics arguments and known asymptotic behaviour for certain known cases. Then using this solution we find the equation which it satisfies, and finally we obtain the corresponding equation for the density matrix. In a sense we solve the inverted problem here. We confront the obtained equation with the standard Liouville equation. Our analysis is applied to the cases of the uniform medium and a medium with adiabatically varying parameters where the eigenstates of the evolution equation are well defined. The case of a non-adiabatic density change and neutrino self-interactions will be considered elsewhere [14]. This allows us to clarify the origin of two (time and coordinate) derivatives in LE. We study how and under which conditions the evolution equation is reduced from the full Liouville form with three spacial and one temporal derivative to a differential equation with only one derivative. Our treatment allows also to clarify form and meaning of the matrix of velocities. The transition from the single particle description to the description of a gas is important for this discussion. Using a modified evolution equation for a single neutrino we derive the equation for the matrix of densities.
The paper is organised as follows. In section 2 we will show that the standard Liouville equation is satisfied for a system of point-like particles. In section 3 we will find the differential equation that corresponds to the wave packet description of propagating neutrinos and confront it with the usual Liouville equation. The obtained equation contains an additional term, as compared to the LE, proportional to the difference of the group velocities. In section 4 we study the physical meaning of this term and consider explicit examples of Gaussian and exponential forms of the wave packet. In section 5 we derive the evolution equations for a neutrino gas. Our conclusions are presented in section 6.

Liouville equation for point-like particles
Let us consider a gas of freely propagating point-like neutrinos. A neutrino emitted at (x 0 , t 0 ) is described by a wave function where τ is the time elapsed since emission, and ψ is a vector in flavour space, ψ T = (ψ 1 , ψ 2 , . . . ). Treated as a point-like particle, the neutrino has a well determined velocity v = (v x , v y , v z ). We assume here the same v for all eigenstates of propagation thus ignoring decoherence. For point-like particles with different velocities the coherence will be lost immediately.
The neutrino flavour evolution is described by a Schrödinger-like equation where the Hamiltonian H is in general a function of the local coordinates (x, t).
For the individual neutrino, we define the density matrix in flavour space as ρ ≡ ψψ † .
Using (2.1) it is straightforward to check that ρ satisfies the equation Let us show that being expressed in terms of (x, t), the density matrix ρ satisfies the Liouville equation. For an individual neutrino, the local coordinates are related to the coordinates of emission: Therefore expressed in terms of the local coordinates, the density matrix can be presented as ρ = ρ (x(τ ), t(τ )) .
Inserting this derivative into (2.2) we obtain which coincides with the standard form of the Liouville equation. The reason of this result is that due to the point-like nature, the local variables x and t are related according to (2.3), and essentially depend on a single variable τ , so that The Liouville operator is nothing but the total derivative over time.

Wave packets
Essentially, the consideration in section 2 was a classical one. In particular, because it uses exact velocities and spatial coordinates simultaneously for the description. Here, instead of the Wigner function, we will use the wave packet description to account for the quantum nature of the neutrinos.
In what follows we will reconstruct the solution of the evolution equation for the density matrix using wave packets. As for the point-like neutrinos, the wave function in position space, ψ T = (ψ 1 , ψ 2 ...) with ψ j being the components in flavour space is given by the Fourier transform: ψ(x, t) = dp (2π) 3 e ip · x φ(p, t).
Here the evolution of φ T = (φ 1 , φ 2 , . . . ) is determined by the Schrödinger equation In a non-uniform medium, the Hamiltonian depends on time t along the trajectory. It is convenient to solve the evolution equation in the basis of eigenstates of the instantaneous Hamiltonian. In this basis where H j are the eigenvalues of the Hamiltonian. In a uniform medium or a medium with adiabatically changing density the evolution equation splits into separate equations for the eigenstates φ j (p, t): The integration of (3.1) gives for the jth component: where f j (p) is the spectral function in momentum space.
In what follows we will neglect the difference in shape of f j (p) for different j: f j (p) = A j f (p), where A j are the numerical coefficients given by mixing parameters at the production. (The shape difference is of the order ∆m 2 /E 2 and it does not increase with time, in contrast to the phases.) Additionally, we assume that in momentum space the spectral function peaks around the characteristic momentum in the packet, Thus, for a uniform medium, the jth component equals We can further specify ψ j by performing an expansion of the eigenvalues H j (p, t) around the characteristic momentum p w : Inserting the expression (3.4) into (3.3), we obtain (see also [15]) is the shape factor which propagates with velocity v wj . This shows that v wj defined in (3.5) can be identified with the group velocity of the j'th eigenstate. In the following the dependence of H on p w will be implicit.
To obtain a form of the evolution equation as close as possible to the standard Liouville equation, we will use the average velocity in the Liouville operator. With this, the Liouville operator acting on (3.6) gives Here so that Y is the traceless matrix proportional to the differences of group velocities and . (3.9) It is important to notice that H, Y and G are diagonal matrices, and the matrix Y can be written as Furthermore, Y is much smaller than |v 0 |. Explicitly, in the two neutrino case and with propagation in vacuum, we find is the third Pauli matrix. Y should be identified with the second flavour asymmetric and traceless part in eq. (3.11).
Computing the time derivative of ρ ≡ ψψ † and using eq. (3.7), we obtain the evolution equation for the density matrix: where is the new term in comparison to the standard Liouville equation (1.1). Being scalar in the position space, C is a matrix in the space of eigenstates.
Since Y and G are diagonal and G is real in realistic examples (see section 4), the correction term (3.13) can be written as the anticommutator (3.14) Consequently, the full evolution equation for the density matrix (3.12) becomes where the dependencies on x and t have been left out for brevity. Thus, in terms of G we obtain a closed equation for the density matrix. G is determined by the production process and therefore can be considered as an external factor (see below). Since both G and Y are real, the correction C is imaginary for real values of ρ, which means that it describes effects of absorption and/or loss of coherence. Eq. (3.15) is one of the main results of our paper, and it takes into account the non-trivial matrix character of V w .
In what follows, we will study properties of C and clarify its physical meaning. Notice that C is proportional to the flavour asymmetric matrix of velocities Y. If velocities are equal for all eigenstates, then Y = 0, and consequently C = 0. So, deviations from the standard Liouville equation is related to the presence of different group velocities in the system.
If C is small compared to [H, ρ], the standard Liouville equation provides a good description of the evolution. Notice that the dominant flavour invariant component of H gives a vanishing commutator and does not contribute to the evolution. The flavour non-symmetric component of H is due to the difference of eigenvalues which depends on masses and background potentials. Therefore the smallness of Y in comparison to H cannot be related to the smallness of the neutrino mass.
Let us estimate C. The general expression in (3.10) leads to the simple estimate |Y| ∼ H traceless /|p w |. In turn, G can be estimated using the size of the shape factor σ x : Thus, the correction C is suppressed in comparison to the main term in (3.12) by the ratio σ p /|p w | ≪ 1, and for a wave packet which is narrow in momentum space, the Liouville equation is a good approximation. However, even in this case the correction from (3.13) can become significant over time, as we will see in section 4.
The matrix nature of V was considered in various derivations of the Liouville equation. In particular, in [2] the following equation has been obtained: 2 Using the expression for the velocity matrix (3.8), it can be rewritten as Let us compare (3.17) with our result in (3.15). For the wave function in eq. (3.6) we find ∂ x ρ = Gρ + ρG * . 2 We still ignore the term i 2 {∂pρ, ∂xH} in [2].
With this, the last term in (3.17) (analogy of our C) becomes which is still different from our expression for C. The two expressions coincide provided that 1) G is real, G = G * , which, indeed, is satisfied for the examples we will discuss.
Under these conditions According to our consideration, condition 2), and consequently eq. (3.18), are not satisfied due to different group velocities even with the assumption of universal spectral functions f . As we will show in section 4, the difference of group velocities leads to separation of the wave packets and loss of propagation coherence.
Finally, we remark that there is also an alternative route to the evolution eq. (3.12). Using eq. (3.6) we can reconstruct the density matrix in the position space immediately: The diagonal elements are 19) and the off-diagonal elements equal Applying the Liouville operator on the matrix formed by (3.19) and (3.20), we obtain (3.12).

Decoherence
Here we study the physical effects of the correction C. To simplify calculations we consider first propagation in a uniform medium or in vacuum. In this case H does not depend on t and consequently Considering two neutrinos allows Y to be written where ∆v ≡ v w1 − v w2 . In vacuum, ∆v =p w ∆m 2 21 /(2|p w | 2 ), which can also be seen from (3.11). Using the definition (3.14) and Y (4.1), we obtain In terms of components of the density matrix, C can be rewritten as In what follows we find C for specific shape factors. For the Gaussian shape factor (which also corresponds to a Gaussian spectral function f (p)) using the definition (3.9), we obtain As we assumed in section 3 G is, indeed, real. Since |x − v wj t| ∼ σ x ∼ 1/σ p , we have |G j | ∼ σ p . This agrees well with the simple estimate used to obtain (3.16).
Inserting expression (4.4) for G j into (4.2) we find According to (4.5), 1) C ∝ σ 2 p ∼ 1/σ 2 x , that is, the shorter the wave packet in the position space, the stronger the effect.
2) C ∝ ∆v -the correction is proportional to the difference of group velocities.
From this, one can conclude that C describes the effect of loss of coherence due to separation of the wave packets. This can also be inferred from the matrix structure of C, where in (4.5) the diagonal part separates the densities (wave packets) of the two eigenstates in space. Indeed, x = v wj t determines the centres of the wave packet, so that x > v wj t is the front and x < v wj t is the back parts of the wave packet. For ρ 11 and the front part, we have −iC > 0. Therefore the 11-term of the matrix increases the density distribution ρ 11 in front and decreases it in the back for v w1 > v w2 , while the 22-term does the opposite. These changes are equivalent to a relative shift of the packets.
The second matrix in (4.5) affects the off-diagonal part of the density matrix. Due to the sign and the factor of (∆v) 2 , it leads to a decrease of the off-diagonal parts of ρ in the uniform medium, thus producing decoherence.
The physical meaning of C can also be seen performing an integration of eq. (3.12) over x: Let us define ρ(t) ≡ d 3 xρ(x, t), and chose the integration volume large enough, so that ρ(x, t) = 0 at the boundary of the volume. Then the integral of the second term on the LHS (which can be reduced by the Stock theorem to the surface integral) vanishes. In a uniform medium H does not depend on x, and consequently eq. (4.6) gives Inserting expression (4.5) into (4.7), we find that the diagonal matrix in (4.5) integrates to zero because the density distribution of the j'th eigenstate is symmetric around x = v wj t. Consequently (4.7) becomes Here L coh = 1 |∆v|σ p is the coherence length and ρ 21 (t) and ρ 12 (t) are again the x-integrated elements. This result can also be obtained by taking the explicit form of the wave packets from (3.6) and (4.3) and calculating ∂ t dxψψ † .
As a second example, let us consider the case of exponential wave packets in one dimension which correspond to neutrinos produced in decays [16]. To make the packets continuous in whole space, we parametrise them as the limit σ f → ∞ of for v wj t < x, (4.9) j = 1, 2. From this expression we obtain

Then according to (4.2) the correction C equals in assumption of
Notice that the off-diagonal terms of C are non-zero only in the second region between the two maxima of the wave packets. The coefficient in these terms, 1 2 (σ p + σ f ), grows unrestrictedly when σ f → ∞. This is not a problem since according to eq. (4.9) ρ 12 and ρ 21 decrease exponentially in the same limit. However, it highlights that the evolution equations only make sense for a continuous g.
The evolution equation for the x-integrated density matrix is given in (4.7). The correction integral can be computed using the expression for ρ 12 in (3.20) with the explicit form of the shape factors in (4.9). In the second region, which gives the correction, we have

Its integration leads to
and in the limit σ f → ∞, we obtain Integration of ρ 12 over the whole space (over all three regions) gives And in the limit σ f → ∞: Using (4.11) we can express the correction (4.10) in terms of ρ 12 (t) in the limit σ f → ∞: which is also valid for v w1 < v w2 . As in the Gaussian case, we see that the diagonal terms disappear as they must preserve probability. The result in eq. (4.12) differs from the result for the Gaussian wave packet case in eq. (4.8) by a factor 2t/|L coh |. Since the effect of wave packet separation is significant at t ∼ L coh , the two expressions give similar results, although the decoherence is less pronounced at early times for the Gaussian wave packet. The transition region is sharper in time in the Gaussian case. The correction term found in (4.12) agrees with the one used in [17], the only difference is a factor due to the definitions of σ p .
The decoherence can also be seen in the explicit form of the density matrix. The offdiagonal elements from (3.20) integrated over x: disappears with time due to decrease of the overlap of g(x − v wj t) and g(x − v wk t), as the consequence of different group velocities .
In a non-uniform medium, ∆v varies and furthermore changes sign close to the MSW resonance [18]. Consequently, the eigenstates might decohere at high densities and later restore partially or completely the coherence at lower density, where the difference in group velocity ∆v = ∂(H 1 − H 2 )/∂p has changed sign. To take such effects into account, we need to take into account the dependence of H on time, and therefore restore the integral over time of H and ∆v. Carrying this change through the calculations for Gaussian wave packets gives the space integrated evolution equation corresponding to eq. (4.8): The presence of both ∆v(t) and t 0 dt ′ ∆v(t ′ ) allows the equation to describe both decoherence and the restoration of coherence when ∆v changes sign, but t 0 dt ′ ∆v(t ′ ) does not.

The neutrino gas. From a single neutrino to an ensemble of neutrinos
The consideration in the previous sections has been done essentially for a single neutrino and the density matrix ρ is simply the matrix representation of the neutrino flavour polarisation vector (see e.g. [19]). Here we present a generalisation of the analysis to the case of a field of propagating neutrinos characterised by the matrix of densities in flavour space n(x, t, v, E). (We can use p instead of v, E). In the absence of inelastic collisions, the matrix of densities n can be reconstructed from the density matrices of individual neutrinos propagating along certain trajectories determined by v. After obtaining n, we will find the equation it satisfies.
Let us start with the simplest case of a flux of point-like neutrinos F (t) emitted from a source at x = 0. We consider first one spacial dimension and relatively small distances, x ≪ L coh , so that the difference of group velocities can be neglected. For a point-like neutrino ρ(x, t) = ρ(x) depends on the distance from the production point, and satisfies the single derivative equation 3 iv∂ It also satisfies the Liouville equation, so that v∂ x can be substituted by D since ∂ t ρ(x) = 0. The density of neutrinos in the point (x, t) is given by It is straightforward to show that n(x, t) satisfies the standard Liouville equation. Acting on n(x, t) by the Liouville operator with single velocity we obtain Here in the second equality we used eq. (5.1) for the derivative ∂ x ρ(x) and in the last equality -the definition (5.2). Notice that in (5.3) the derivatives of F cancel: due to the dependence on coordinates in the combination (t − x/v).
For an ensemble of point-like neutrinos, the time dependence follows from the production condition and it is independent of the x-derivative. As a consequence of the point-like character of neutrinos, there is no restriction on the speed of time-variations. The flux factors out from the evolution.
Notice that in more than 1D the divergency of the flux in space is taken into account by the LE which implies conservation of the flux in the complete phase space which includes both coordinates and momenta.
In the three dimensional case neutrinos are emitted from a surface fixed by the coordinates x 0 . In general, the flux can depend on the emission point on this surface and on the direction determined by v. So, F = F (x 0 , t 0 , v), and in terms of coordinates (x, t): where v = |v| and L(x) = |x − x 0 (x)|. 3 Formally we could introduce the time dependence in ρ(x) as This density, indeed, satisfies the LE and leads to the same equation (5.3) for the matrix of densities n.
Acting by the 3D Liouville operator (with universal v) on the matrix of densities where ρ(L) = ρ(L(x)) is the density matrix of an individual neutrino propagating along the trajectory with a baseline L and satisfies the equation Here we have taken into account that The action of the Liouville operator on F gives where differentiation of t 0 = t − L/v vanishes, as in (5.4). The term in (5.7) is due to dependence of the flux on coordinate of the emitting surface. Let us show that this term vanishes too. Given coordinates (x, v) determine the emission point x 0 on the emission surface. If the surface is continuous and differentiable around x 0 , we can consider a plane tangent to the surface at x 0 , and take a small patch of the plane around x 0 . Then the surface is well approximated by the tangent plane. We can consider the system of coordinates x = (x (1) , x (2) , x (3) ) in which the tangent plane is at x ) (that is, the trajectory of the neutrino is in the plane with x (2) = 0). In this coordinate system and then immediately, Consequently, the RHS of (5.7) is zero, DF = 0, and eq. (5.5) becomes Finally using eq. (5.6) we again arrive at the standard Liouville equation for n: Let us consider the case of individual neutrinos described by wave packets. Now three scales are relevant for the problem: The width of the wave packet in position space, σ x , the effective distance between the centres of the wave packets of neighbouring neutrinos, λ, and the coarse graining scale, ∆x, in which we will compute the number of neutrinos. Then the density is obtained dividing the number of neutrinos by ∆x. For simplicity we will consider one spatial dimension, so that the distance between the centres of the wave packet equals λ = v∆t and ∆t is the averaged interval between emission of two sequential neutrinos. 4 Then the matrix of densities is where t j is the time of emission of the j'th neutrino. We take ∆x ≫ λ, so that in the interval ∆x there are many neutrinos, otherwise, the description of the gas is not different from the description of every particle independently. Depending on σ x , we have two different physics cases: 1) If ∆x ≪ σ x , the density matrix is almost constant on the interval [x, x + ∆x], it can be put out of the integral in (5.8), and the integration is trivial leading to The number of wave packets per coarse graining scale equals n x = ∆x/λ, so that ∆x = n x λ = n x v∆t. Inserting this expression into (5.9) and substituting summation by integration Since the integral over ρ is in an integral over shape factors with the form g(x − v wj t), it is very well approximated by where it was used that the shape factors give a significant contribution to the integral when t ′ ≈ t − x/v. The number of neutrinos per ∆x can also be expressed in terms of the flux as F/v. If the flux changes slowly with time compared to ∆x/v this does not affect the integral, and We notice that the expression in (5.11) has the same form as the n for point-like particles in (5.2). As a consequence, the derivation is very similar to the derivation for point-like particles leading to (5.3) with the difference that (5.1) should be replaced by Eq. (5.12) is the 1D version of eq. (4.7) with the substitution t → x/v. The evolution equation for n(x, t) in (5.11) is then where C is given in (3.14). The last term describes the wave packet decoherence as it has already been discussed in section 4.
2) In the case of a narrow wave packet, λ, σ x ≪ ∆x, the matrix of densities can be obtained from (5.8): (5.14) Here it was used that almost all of the wave packets are fully contained in the interval, so that the integration limits can be extended to infinity. The number of wave packets in the interval ∆x is n x , and the emission time for wave packet at (x, t) is given by t j = t − x/v as we also used in (5.10). Since (5.14) gives the same expression for n(x, t) as in (5.10), the evolution equation for the matrix of densities also coincides with eq. (5.13). Thus, we have demonstrated that the basic form of the Liouville equation for the density matrix can be obtained without referring to the Wigner functions, and hence it is be independent of the associated peculiarities.
6 Conclusions 1) We study the validity and applications of the Liouville equation for oscillating neutrinos in vacuum and in matter. Inelastic interactions and neutrino self-interactions were neglected here. In contrast to the previous studies we have accounted for the quantum nature of neutrinos using explicitly the wave packet description of the neutrino state avoiding the use of Wigner functions.
2) We have reconstructed explicit solutions for the wave functions and then the density matrix using known exact solutions. Then we have found equations which such a density matrix satisfies and confronted them with the standard Liouville equation. The obtained equation differs from the Liouville equation by an additional term C which is related to the presence of components with different group velocities (essentially, with different eigenvalues of the Hamiltonian). We show that C describes the effect of loss of propagation coherence due to separation of the wave packets. The correction C depends on the specific shape of the wave packet which, in turn, is determined by production conditions. It is proportional to t/2L 2 coh for Gaussian and to 1/4L coh for exponential shape factors. With C we put certain information about the initial conditions into the evolution equation. This is consistent with the fact that the density matrix itself is constructed using the wave packets.
3) In the wave packet approach the matrix v equals the unit matrix multiplied by the average group velocity of the eigenstates in the system. The difference of the group velocities appears in the additional term C in the equation. 4) We show that the density matrix constructed with wave packets, and thus explicitly accounting for the quantum nature of neutrinos, satisfies the modified Liouville equation with only one additional term which describes separation of the wave packets. The latter is unavoidable in the wave packet picture. This also implies that one can simply consider the standard Liouville equation for the plane waves with equal momenta for all eigenstates and then perform integration over the production spectrum and the detection energy (momentum) resolution. In this case the density matrix depends on time only; derivatives over coordinates and velocity are irrelevant. 5) Our approach allows to clarify the meaning of derivatives over time and coordinates in the Liouville operator. For an individual point-like neutrino the density matrix depends essentially on a single variable (time or coordinate along the trajectory related by the equality L = vt). The Liouville equation is satisfied trivially since it can be reduced to the Schrödinger-like equation.
6) If inelastic collisions are absent, the field of neutrinos can be described as the sum over individual freely propagating neutrinos. In this case additional dependence of the matrix of densities on time and coordinate may follow from dependence of the produced neutrino flux on time and the point on the emission surface. We show that in this case the matrix of densities still satisfies the modified LE. The additional dependence does not affect the evolution equation and can be taken in to account as an external factor. 7) Additional dependence of time may follow from dependence of the matter potential on time. If this dependence is not very fast, it simply corresponds to an adiabatically changing H. The special case of very fast time dependence or non-adiabatic variation of density can be realised in the presence of neutrino self-interactions. This case will be considered elsewhere.