Inhomogeneous Kinetic Equation for Mixed Neutrinos: Tracing the Missing Energy

Flavor-dependent neutrino transport is described by a well-known kinetic equation for occupation-number matrices in flavor space. However, as an overlooked theoretical problem, we show that in the inhomogeneous case, neutrino-neutrino refractive energy is not conserved. We derive the missing gradient terms in the fast flavor limit (vanishing neutrino masses), and prove that the missing refractive energy is traded with the huge reservoir of neutrino kinetic energy through gradients of the weak interaction potential. Even small changes of the kinetic energy accommodate the refractive energy gained or lost. Flavor evolution alone is negligibly affected by the new terms.

Predicting the consequences of fast conversions is thus a central question.In tackling this problem, a key role must be played by conserved quantities, which are the only exact guide into the final state reached after conversions.One such quantity must be energy.However,we show that the usual equations of motion (EOMs) do not conserve energy.In their often-used form, they are where H p,r,t = √ 2G F p ′ ρ p ′ ,r,t (1 − cos θ p,p ′ ) is a matrix driving flavor evolution by neutrino-neutrino refraction.The advection term proportional to v = p/|p| quantifies a drift in coordinate space caused by inhomogeneities.This term allows for p-dependent instabilities, so even an initially homogeneous system can develop flavor disturbances growing exponentially and drifting to ever smaller scales [15,[29][30][31][32][33][34].
As detailed later, even an elementary example of two colliding beams consisting initially of ν e and ν x reveals that the refractive interaction energy changes dramatically as soon as initial seeds of inhomogeneity grow.The solution of this puzzle is that one cannot look at flavor evolution in isolation.Inhomogeneities in the weak potential exert a force and trade refractive interaction en-ergy with the huge reservoir of neutrino kinetic energy.We here derive the missing terms in Eq. ( 1) and show that to lowest order in a gradient expansion, the energy exchange can be exactly accounted for.
Interaction energy.-Wefirst define the ν-ν interaction energy.The starting point is the many-body Hamiltonian H = H 0 + U, where H 0 = α,p ϵ p a † α,p a α,p is the kinetic energy operator with ϵ p = |p| in the massless limit.Moreover, the interaction energy is where G F is Fermi's constant, {p} is performed over all momenta such that p 1 +p 3 = p 2 +p 4 , and u p is the spinor of a left-handed massless neutrino with momentum p, normalized such that u p γ 0 u p = 2.The energy of the system is the expectation value of H over the quantum state.In particular, the kinetic energy is K = ⟨H 0 ⟩ = d 3 r p ϵ p Tr(ρ p ).For the interaction energy, we need to determine the average value of a string of four creation and annihilation operators.We use mean field approximation, where the quantum state is assumed to be the product of single-particle states.The expectation value factorizes as In this way we define a potential energy of interaction U = ⟨U⟩.For a homogeneous system, it is where we neglect an irrelevant (infinite) renormalization of the neutrino chemical potential.
In the two-flavor case, the density matrices are often decomposed into their trace and trace-free part as

arXiv:2401.05278v2 [hep-ph] 30 Jul 2024
where 1 is a unit matrix, P 0 p = Tr(ρ p ), ⃗ σ a vector of Pauli matrices, and ⃗ P p a polarization vector.( ⃗ P is a vector in flavor space, p one in phase space.)Exactly one particle in mode p implies P 0 p = | ⃗ P p | = 1.Only ⃗ P p evolves by oscillations, whereas the trace (total particle number in a p-mode) is conserved.The oscillatory part of the interaction energy is corresponding to that part of Eq. ( 3) that depends on the trace-free parts of the ρ matrices.
For simplicity we assume in the following a system that remains homogeneous in all but the spatial z-direction and use v = v z .The standard EOM Eq. ( 1) falls into pieces for particle number and flavor polarization as where we use the angular moments ⃗ P 0 = p ⃗ P p and ⃗ P 1 = p v ⃗ P p .Thus, the usual equations decouple the neutrino number from their flavor evolution.
Two-beam example.-Thenonconservation of U osc is a large effect as we show in a simple example of two opposite beams along the z-axis (v = ±1).The energies and number densities are equal, one beam initially occupied with ν e , the other with ν x .As usual, the beams represent many neutrinos with p's close enough that the small spread is unimportant on the relevant time scales.Therefore, we represent the flavor polarization of each beam by a highly occupied ρ matrix and concomitant polarization vector ⃗ P ± (z, t).We thus need to solve where the interaction strength was absorbed in the units of time and space.For ∂ z ⃗ P ± = 0 there is no instability, but the system has unstable inhomogeneous modes.
We can show energy nonconservation analytically in the initial linear regime.In this limit, it is convenient to express the x-y part of the polarization vector (the off-diagonal element of the density matrix) in the form ψ ± = P x ± + iP y ± and initially we take P z ± (z, 0) = ζ ± to be constant.The linearized EOMs (|ψ| ≪ |ζ|) are These EOMs are most easily solved for the spatial Fourier modes for a periodic box of length L where k = 2πn/L with integer n is a discrete wave vector.Equation (9) shows that a k mode is unstable if The oscillation energy ), where in the continuous limit k = dk/2π.This piece can grow exponentially if there are small seeds for the unstable Fourier modes.If we decompose the v = ±1 modes in terms of the eigenmodes of the linear analysis, the k modes evolve as Matching this expression with the initial conditions ψ ±,k (0) = ψ 0 ±,k reveals for the growing piece (11) Substituting in the transverse part of the oscillation energy, we finally find where we restrict the summation to the range of unstable eigenmodes.Notice that the most unstable eigenmode k = k does not contribute to energy nonconservation due to the vanishing prefactor.Otherwise, an exponential growth appears already at the linear level.
For more than two beams, from Eq. ( 7) it follows that , and thus it is still true that the energy change grows with e 2γt , provided that k ̸ = 0; homogeneous modes, with k = 0, do not lead to energy change due to the prefactor kv.
For a numerical solution that extends to the nonlinear regime, we use periodic boundary conditions on a box of length L = 100, divided in N spatial bins.As initial conditions, we use P z ± (0) = ζ ± = ±1/2.These parameters imply that the range of unstable modes is delimited by k 1 = −2 and k 2 = 0, the maximum growth rate obtains for k = −1 and is γ = 1, and the initial oscillation energy is ± are seeded with randomly sampled functions with a spatial dependence where the amplitudes are sampled from a normal distribution with a variance σ 2 = 10 −8 and the phases are uniformly sampled from 0 to 2π; the functions are real and so c n = c −n and ϕ n = −ϕ −n .We use arbitrarily N max = 100 to avoid seeds at too small scales.Numerical stability requires keeping fluctuations at the scale of the grid spacing as small as possible.
Figure 1 shows the solution for a single realization of initial seeds.The upper panel shows contours of the ν xcontent of the beam initially occupied with ν e , as a function of z and t.The lower panel shows U osc (t) in units of U osc (0).On the shown timescale, the range of unstable modes −2 < k < 0 grows nonlinear and then keeps oscillating, with beats between different modes causing an irregular pattern that in detail depends on the choice of seeds.In the longer run, these nonlinear modes feed higher-k modes and flavor variations will obtain on ever smaller scales.If we had used a larger box, or equivalently averaged over more than one realization of the initial conditions, the amplitude of the oscillations in the final state would shrink.Nevertheless, the bulk of the nonconservation of U osc happens in the initial phase, and is clearly a large effect.
Inhomogeneous kinetic equations.-Thelack of refractive-energy conservation questions the validity of the traditional EOMs.The strategy for their derivation is a perturbative expansion in three small parameters: (i) The mass-to-energy ratio r m = m ν /ϵ p , where m ν is the neutrino mass, but in our fast flavor limit we do not worry about it.(ii) The fractional refractive energy shift r µ = µ/ϵ p , where µ = √ 2G F n ν is a scale for the ν-ν interaction energy.For T ≃ 5 MeV in the decoupling region of a supernova, n ν ≃ 10 33 cm −3 , µ ≃ 0.2 meV, and r µ ≃ 10 −11 and thus indeed very small.(iii) The scale of density variations ℓ in units of the neutrino de Broglie wavelength r ℓ = (ℓ|p|) −1 ≪ 1. Density variations are on scales much larger than |p| −1 .The EOMs are typically derived to lowest order in all of these small ratios, notably neglecting terms of the order of r ℓ r µ , which causes the nonconservation of energy.
To augment Eq. ( 1) with the missing terms, we note that in an inhomogeneous setting, the space-dependent occupation-number matrix is defined by a Wigner distribution [49], the expectation value of Here, the Fourier wavevector k measures the typical length scale over which the density matrix is changing, and is therefore |k| ∼ ℓ −1 ≪ |p|.The time derivative is obtained using Heisenberg's equations.The commutator with the kinetic energy is easily established as recovering the usual advection term.
For the interaction part, we here outline the general strategy and report details in the Supplemental Material (SM).The commutator [D, U] produces strings of four creation and annihilation operators, or four-point correlation functions.The average of these combinations are evaluated again using the mean-field approximation.In the homogeneous limit, when taking the expectation value of ⟨a † 1 a 2 ⟩, the momenta of state 1 and 2 must be equal.We go one step further and include the inhomogeneity to first order in k.The final result, fully derived in the SM, is where we introduce that coincides with the standard H p,r stated after Eq. ( 1) except for the additional trace term that drops out in commutators, but not in the anticommutators on the left-hand side.Moreover, where p = |p| and p ′ = |p ′ |.
Ω p,r is the renormalized quasiparticle energy, in general a matrix in flavor space, analogous to the renormalized quasiparticle energy in Landau's theory of a Fermi liquid [50].This induces a renormalization in the group velocity ∇ p Ω (0) p,r (only the zero order terms must be kept which already contain one gradient) and a weak force field −∇ r Ω (0) p,r .All additional terms are small, of the order of r µ r l , and thus will not produce quantitatively large flavor-conversion effects.However, as a qualitatively new effect, the renormalized group velocity causes a slow spatial drift of neutrinos; the order of magnitude of the velocity is |∇ p Ω (0) p,r | ∼ r µ .Since r µ ∼ 10 −11 , this is completely negligible compared to the standard neutrino velocity.
The most significant impact comes from the weak force −∇ r Ω (0) p,r , which can change the neutrino kinetic energy.Part of this force originates from gradients in the neutrino number density.Our new insight is that an additional component originates from gradients in flavor composition.To see this clearly, we rewrite the EOMs as where is the number flux of neutrinos passing through a surface element of coordinate space at phase-space location {p, r}, as seen from the structure analogous to a continuity equation, whereas F p,r is the same in momentum space, given by the same expression with ∇ p → ∇ r .This term couples the trace of the density matrix with the polarization vector; taking the trace, we find Thus, spatial gradients of the polarization vectors can feed ν kinetic energy.This term implies a net rate of energy gain or loss of order dϵ p /dt ≃ µ/ℓ.Over a timescale ℓ, the energy that a ν can accumulate is of order µ ≪ ϵ p , so we recover that this is a small effect relative to ϵ p , but large relative to the interaction energy.It provides the missing channel by which kinetic and refractive energy can be traded.Even a small amount of energy µ lost (gained) from the large ν kinetic energy explains the large change in the refractive energy, since U osc /K ≃ µ/ϵ p .
With the new kinetic equations, we can directly prove that to first order in r µ and r ℓ , the total energy is conserved.After integrating by parts, we find Substituting the polarization vector part of Eq. ( 21) in Eq. ( 22), we find that it precisely balances the rate of change in the oscillation energy found in Eq. ( 7).
Discussion.-We have shown that the usual kinetic equation for mixed neutrinos Eq. ( 1) is non-conservative in the inhomogeneous case.We have illustrated this point with a simple two-beam example, where the ν-ν interaction energy strongly changes, and we have also shown this effect analytically in the linear regime.Deriving missing gradient terms beyond Eq.( 1), we have shown that the refractive energy gained or lost is precisely traded with neutrino kinetic energy that usually is not followed.In our two-beam example, the monochromatic initial energy distribution develops small space-time dependent shifts that account for the missing energy.We have not tried to study this effect numerically because it requires many p-modes around the original one, but we have proven energy conservation analytically.
The usual EOMs correctly account for flavor evolution, meaning the trace-free part of the density matrices often described by polarization vectors.On the other hand, the particle number in a given p-mode is conserved, corresponding to the trace part of the EOM, where the lefthand side of Eq. ( 1) is a continuity equation.Therefore, the trace part of the EOM must be expanded to higher order in the gradients to capture nontrivial evolution.The reshuffling of neutrinos among p-modes is a small effect relative to the large kinetic energies, yet precisely absorbs the missing refractive energy.
These novel terms come from the gradients of the neutrino self-energy.With hindsight, that they would affect neutrino evolution, is obvious from the viewpoint of physical kinetics, and implicitly present in previous formal derivations [2,3,[5][6][7], but usually assumed to be a small effect.Our new insight is that fast conversions spontaneously break homogeneity, magnifying the gradients of the trace-free density matrix, making these terms large enough to explain completely the apparent nonconservation of energy.
We note that Eq. ( 1) conserves entropy and including the new gradient terms, this is also the case, i.e., entropy is conserved order by order in the gradient expansion (see Supplemental Material).On the practical level, the conservation of entropy, but not of energy, could be used to test numerical stability.Even more importantly, in making predictions on the final state induced by conversions, energy conservation cannot be used in practice, unless the new terms are kept.
Our finding may shed new light on a recent proposal that the final outcome of fast conversions may be some sort of thermalized state [51].The quasi-steady state generically observed in numerical simulations of fast conversions [52][53][54][55][56][57] is only determined by the oscillatory part of the density matrix, yet its dynamics does not admit a separately conserved energy.Thus, it cannot separately thermalize, since it exchanges energy with the kinetic energy of neutrinos.

Supplemental Material for the Letter Inhomogeneous Kinetic Equation for Mixed Neutrinos: Tracing the Missing Energy
In this Supplemental Material, we provide further details on the derivation of the equation of motion (EOM) discussed in the main text.

A. Detailed derivation of the equations of motion
For the reader's convenience, we reproduce some of the equations from the main text.The starting point of the derivation is the many-body Hamiltonian describing the neutrino kinetic energy and interactions.This is written H = H 0 + U, where H 0 = α,p ϵ p a † α,p a α,p is the neutrino kinetic energy operator, and We will henceforth use the notation (p 1 , p 2 )(p 3 , p 4 ) = u p1 γ µ u p2 u p3 γ µ u p4 .We want to determine the time evolution, under this Hamiltonian, of the density operator or more conveniently of its spatial Fourier transform (S4) Expanding to first order in k, and using ∇ p ϵ p = v, we find which Fourier transformed leads to the usual advective term quoted in the main text.
Let us now consider the commutator of the interaction Hamiltonian U, which gives where δ p1,p2 is a Dirac delta forcing the equality p 1 = p 2 .We can now take the expectation value of this operator over the mean-field state; the corresponding terms can be grouped in pairs which are identical after the exchange p 1 → p 3 , p 2 → p 4 , so we finally obtain In each term, it will prove convenient to reparametrize the momenta in such a way as to separate the large momenta (associated with the neutrino energies) from the small momenta (associated with the spatial gradients of the density matrix).In particular, in terms of the form ⟨a † β,p1 a α,p2 ⟩, because of the slowly-varying nature of the density matrix, we find a non-vanishing result only if p 1 and p 2 differ by terms of order ∼ ℓ −1 .Thus, for example, it proves convenient to reparameterize the first term ⟨a † γ,p1 a γ,p2 ⟩⟨a † α,p3 a β,p−k/2 ⟩ choosing p 1 = p ′ + (k − q)/4, p 2 = p ′ − (k − q)/4, p 3 = p + q/2, p 4 = p + k/2, which makes it clear that p and p ′ are large momenta, while q and k are small.Adopting a similar parameterization for all the terms we find We can now take the inverse Fourier transform by multiplying by e −ik•r1 and summing over all k.We perform the full manipulations only on the first of the four terms; the other ones are treated analogously.Introducing the notation a = k − q, b = k + q, we can write this term as In this expression, we must now perform the expansion for |k, q| ≪ |p, p ′ |.For the spinor part, this requires using the explicit form of the spinors expanded to first order in the momenta a and b.The result depends on the mutual orientation of the vectors p and p ′ .Thus, we introduce a set of three orthonormal directions The expansion in powers of k and q involves three terms.The zeroth-order term is then simply

4 FIG. 1 .
FIG.1.Solution of our two-beam example for one realization of initial seeds.Top: Contours of νx content of v = +1 beam as a function of z and t.Bottom: Evolution of oscillation energy Uosc(t) in units of Uosc(0).

)
Following the notation in the main text, we denote the expectation value over the quantum state ⟨D αβ (p, r)⟩ = ρ αβ (p, r) and analogously for the Fourier transform ⟨ Dαβ (p, k)⟩ = ραβ (p, k).We start by obtaining the commutator of the non-interacting Hamiltonian H 0 with this operator; using the standard commutation rules we easily find [H 0 , Dαβ (k, p)] = ϵ p+k/2 − ϵ p−k/2 Dαβ (k, p).