Quantum transport theory for neutrinos with flavor and particle-antiparticle mixing

We derive quantum kinetic equations for mixing neutrinos including consistent forward scattering terms and collision integrals for coherent neutrino states. In practice, we reduce the general Kadanoff--Baym equations in a few clearly justified steps to a generalized density matrix equation that describes both the flavour- and particle-antiparticle coherences and is valid for arbitrary neutrino masses and kinematics. We then reduce this equation to a simpler particle-antiparticle diagonal limit and eventually to the ultra-relativistic limit. Our derivation includes simple Feynman rules for computing collision integrals with the coherence information. We also expose a novel spectral shell structure underlying the mixing phenomenon and quantify how the prior information on the system impacts on the QKE's, leading to a direct effect on its evolution. Our results can be used for example to accurately model neutrino distributions in hot and dense environments and to study the production and decay of mixing heavy neutrinos in colliders.


Introduction
One of the key long term goals in neutrino physics is to obtain practically solvable quantum kinetic equations (QKEs) that can accurately and consistently model the coherent neutrino flavor evolution including also decohering collisions.Traditionally neutrino oscillations have been described either in the quantum-mechanical (QM) wave-packet approach or using some partial quantum-field theory (QFT) methods.In the QM treatment [1][2][3][4][5][6][7][8][9][10][11] neutrinos are described by wave packets whose widths are related to the uncertainty of neutrino momentum at the production process.In the partial QFT approaches [11][12][13][14][15][16][17][18][19][20][21][22][23] neutrino production, detection and propagation are considered via compound Feynman diagrams that treat neutrinos internally as QFT propagators and externally as wave packets.In vacuum, when decoherence effects related to neutrino production, detection and propagation can be neglected, both approaches lead to the standard formula for neutrino oscillations probabilities.However, in hot and dense environments, where interactions modify neutrino flavor evolution significantly, more fundamental approaches are needed.
Other efforts to derive quantum kinetic equations in the specific neutrino physics context include [48][49][50][51], but these treatments were less general than that of [40][41][42][43][44][45][46][47], relying on the ultra relativistic (UR) limit and the mean field limit [48] or expansions in small perturbative quantities [49][50][51].These articles found some effects of the particle-antiparticle mixing dubbing them as spin coherence [49] or helicity coherence [48], and proposed that they may be relevant in hot and dense environments (see also [52]).This is usually not the case however and there indeed seems to be some confusion in the literature concerning the notion of particle-antiparticle mixing.This term seems to be used when actually referring to the CP-violating flavour mixing induced effects on the evolution of the particleantiparticle asymmetry, already included in the early treatments discussed above.Such CP-violating effects are of course at the heart of the leptogenesis [53] and the electroweak baryogenesis [54] problems.They are also relevant for the production and decay of heavy Majorana neutrinos in high energy physics experiments [55] as well as for the dynamics of the supernovae explosions [52].The true particle-antiparticle mixing on the other hand, is relevant e.g. for particle production in the early universe during the (p)reheating phase after inflation [56,57].
Articles [41][42][43][44][45][46] discussed dispersive corrections only at a generic level.A more detailed formulation, but still with no explicit expressions for dispersive corrections was given in [53].These corrections are important as they eventually give rise to neutrino forward-scattering potentials.A fully general and self-consistent derivation of the QKEs from fundamental principles, which include both the forward scattering potentials and the decohering collision integrals and encompass both flavour and particle-antiparticle mixing coherences, has still been missing until now.This work fills this gap responding to a demand sometimes voiced in the literature [52].Our formalism is not restricted to treating just neutrino mixing and coherences and it is indeed much more general than is necessary for most neutrino physics applications.We will still tune our presentation at all times keeping the neutrino physics motivation in mind however.
We start from a SD-equation in the Closed Time Path (CTP) formulation [58][59][60].In a series of well motivated steps reduce the SD-equations to quantum kinetic equations which contain all coherence information essential for neutrino mixing and oscillations.Our derivation assumes only slowly (adiabatically) varying background fields, the validity of weak coupling expansion and eventually the spectral limit (although this is optional).We employ the local approximation that was recently developed and applied in the context of resonant Leptogenesis [53] and is closely related to the coherent quasiparticle approximation (cQPA) developed in [41][42][43][44][45][46][47].An essential element in our derivation is the introduction of a projective representation which reduces the original SD-equation to a set of Boltzmann-type equations for distribution functions classified according to their natural oscillation frequencies.In the weak coupling limit it also reveals a novel spectral shell-structure with new "coherence shells" that are recognized to carry information about the particle-antiparticle and flavour coherences.We generalize the work of refs.[40][41][42][43][44][45][46] in several ways.In addition to using the projective representation of [47,53] we allow for an adiabatic evolution in the spatial coordinates.We also include a derivation and evaluation of the neutrino forward scattering potentials, whose general structure is much more complex than that found in [48,49].
Our master equations take a very elegant and intuitive form of a set of Boltzmanntype equations (or a generalized density matrix equation) that fully contain the neutrinoantineutrino mixing and are straightforward to solve numerically.Our results are also valid for arbitrary neutrino masses and kinematics, while e.g.[48][49][50][51] assume UR-limit from the outset.However, we do present our results also in the limit when particle-antiparticle mixing can be neglected and eventually in the UR-limit.Moreover, we take into account some higher (infinite) order gradient corrections in the Kadanoff-Baym-equations, which seem to be neglected elsewhere [49,50], but are necessary for correct evaluation of the self-energies and collision terms.In addition (going also beyond the treatments in [40][41][42][43][44][45][46][47]53]) we present our results clearly indicating at each step how they can be generalized to include finite width corrections beyond the spectral limit.Lastly, our derivation comes with a comprehensive set of generalized Feynman rules which provide a straightforward and systematic way to compute collision integrals for our QKEs including all coherence effects.
Our most general equations are useful to model the production of flavour-mixing particles by background fields, for example at the reheating phase after the inflation or during some phase transitions.In most other cases the particle-antiparticle mixing can be neglected however, and the frequency-diagonal limit can be assumed.This is sufficient for the resonant leptogenesis problem as well as the heavy (Majorana) neutrino-antineutrino oscillations mentioned above.Finally, the very simple UR-limit equations are useful tool to set up numerical framework to study neutrino distributions in hot and dense astrophysical environments, such as neutron stars and compact object mergers.
The paper is organized as follows.In sections 2 and 3 we set up the Kadanoff-Baym equations, and show how the pole-and statistical equations can be decoupled.In sections 4-6 we reduce these equations to local QKEs which take the form of a density matrix equation with and without the particle-antiparticle oscillations.In section 7 we show how to compute collision integrals appearing in these QKEs, and derive simple Feynman rules to automatize this task.In section 8 we compute the general 1-loop forward scattering potentials.In section 9 we discuss how the localization is related to the (lack of) prior information on the system, and how in general a preparation of the system affects its evolution.Finally, section 10 contains our conclusions and outlook.

Kadanoff-Baym equations
In this section we briefly review the derivation of the Kadanoff-Baym equations for the neutrino two-point function from the CTP-formalism.Indeed, the quantity that holds the information about coherence for mixing neutrinos in out-of-equilibrium conditions is the 2-point correlation function: where ρ is an unknown density operator for the system, ψ is the fermion field, T C is time ordering operator and u 0 and v 0 are complex time arguments on the usual Keldysh-contour C [59].The path ordered 2-point function S C (u, v) obeys the Schwinger-Dyson equation [60][61][62][63] (we suppress the Dirac and flavor indices when there is no risk of confusion): where S −1 0 is the free inverse fermion propagator, The self-energy function Σ C depends on the model.It can be computed for example from the 2-PI effective action: where Γ 2 is the sum of the 2-PI vacuum graphs of the theory, truncated to a desired order in coupling constants.
The complex-time SD-equation (2.2) is equivalent to a coupled set of Kadanoff-Baym equations [64] for the real-time valued correlation functions: (2.4) Here p = r, a refers to the retarded and advanced (pole) functions and s =<, > to the statistical Wightman functions.For more details on precise definition of the real time functions and self-energies see [41,42,47,53].The convolution is now defined over real time variables, and one can separate the internal and external degrees of freedom in (2.4) by performing the Wigner transform: where x ≡ (u + v)/2 is the average coordinate and k is the conjugate momentum to the relative coordinate r = u − v.This leads to the following mixed space equations: where x and we defined a shorthand notation for the mixed space correlation function [47]: Here the superscript Σ indicates that the gradient ∂ Σ x acts only on the self-energy function, while ∂ k is a total derivative.We also defined , and absorbed the mass term into the singular part of the Hermitian self-energy: Σ H (k, x) = Σ H,sg (x) + Σ H,nsg (k, x).The two forms of the KB-equations: (2.4) and (2.6) are equivalent and exact to the given approximation for the self-energy function.Each for has its unique advantages that we shall use in what follows.

Decoupling of the pole and the statistical equations
In addition to their manifest non-locality the KB-equations (2.4) and (2.6) feature a direct coupling between the statistical and the pole functions.In order to reduce them to a single local quantum kinetic equation, we must both localize them and decouple the pole equations from the statistical ones.We will address the decoupling problem first.The key idea is to split the statistical function into a background part that is strongly coupled to the pole functions and to a perturbation, whose equation formally decouples.The formal decoupling then suggests a wide range of approximate solutions that make the decoupling exact, leading to a single self-consistent equation for the perturbation.

The general background solution
The pole functions S r,a can be expressed in terms of the Hermitian function S H and the spectral function A: S r,a = S H ∓ iA, and similarly Σ r,a = Σ H ∓ iΣ A .We can then write the pole equations in (2.4) symbolically as follows: while the statistical equation becomes The inverse Hermitian operator in these equations is defined as We chose to work explicitly with S < .The equation for S > is equivalent, but not needed because A = i 2 (S > + S < ).We chose to use the direct space notation here, but one may go to Wigner space by basically just replacing " * " by "⊗" everywhere.
We observe that while the pole functions appear explicitly in the equation for S < , the converse is not true; statistical functions affect the pole equations only indirectly through the self-energy functions.This suggests dividing the statistical functions as S < ≡ S < 0 + δS < , where the background solution S < 0 satisfies the equation This implies that the perturbation δS < satisfies the equation where the collision integral is given by The essential feature of this construction is removing the direct coupling term Σ < * S H from equation (3.5) for the perturbation δS < .We stress that S < 0 is not guaranteed to be the background solution in the sense that it would make the collision term vanish in the equation for δS < .For that to be true additional constraint needs to be imposed which is easiest to see by going to the Wigner space.From (3.1) and (3.4) one readily finds the solutions A = S H0 ⊗ Σ A ⊗ S H and S < 0 = S H0 ⊗ Σ < ⊗ S H . Using these results we can write the Wigner space expression for the collision integral (3.6) as follows: This expression vanishes identically if where g < (k) is an arbitrary 4-momentum dependent scalar function.Equation (3.8) defines a large class of consistent background solutions, including the vacuum: g < vac (k 0 ) = θ(−k 0 ) and the thermal background: 1 g < th (k 0 ) = 2f FD (k 0 /T ).We stress that while we call δS < a perturbation, it does not need to be small; we have only made a convenient division of the solutions, but no approximations yet.All equations written so far are as general as the full interacting field theory itself.

Adiabatic background solutions
The separation of equations (3.1) and (3.4) from (3.5) suggests a way to construct efficient approximations.First note that the pole functions S H and A are strongly constrained by the spectral sum rule, which prevents them from having rapidly varying coherence solutions [40,53].In many cases they can be solved in an adiabatic approximation which then must hold by construction also for the background solution S < 0 .Taking a cue from the exact solutions for A and S < 0 used in (3.7), we may define adiabatic solutions in the Wigner space as follows: A ad ≡ S H0,ad Σ A,ad S H,ad and S H,ad ≡ S H0,ad (1 The idea is that after the division S < = S < 0,ad + δS < , the term δS < should describe a transient around the adiabatic background solution S < 0,ad .Again, this is not guaranteed to hold automatically and for the consistency of the definition two additional conditions are needed.First, we interpret Σ < ⊗ S H as a coherence damping term that only gives the width to the background solution.That is, we set Σ < ⊗ S H → Σ < ad S H,ad in equation (3.2) as a part of the adiabatic approximation.Second, we require that the collision integral (3.6) vanishes for the adiabatic solution.Dropping all gradients in Wigner space convolutions (2.7) (note that when acting on adiabatic solution also K0 → k 0 ), we get from (3.7): This expression vanishes if which is the adiabatic generalization of (3.8).
Beyond the choice of g < ad (x, k), there is a lot of freedom in defining the adiabatic pole solutions.It is often sufficient to work in the spectral limit, where Σ A,ad ≡ 0, or even in the vacuum, where Σ p ad ≡ 0, but taking the finite width, dispersion and even backreaction from δS < could be accounted for, as long as the constraint (3.11) holds.There is yet more freedom in choosing the approximation for the operator S −1 H0 which can differ in the perturbation equation (3.5) from S −1 H0,ad in the pole and background equations (3.1) and (3.4).Moreover, the approximations for the convolutions involving the Hermitian self-energy functions can be different from the leading order result employed for the absorptive self-energy functions leading to (3.11).There are thus many ways to define consistent approximation schemes for the split equations, whose accuracy will obviously depend on the problem at hand.

Decoupled equation for the statistical function
After an approximation scheme is defined for the pole and the background functions, the equation (3.5) for the perturbation δS < acquires additional source terms.Given a specific scheme that satisfies (3.9) and (3.11) and another specific definition for S −1 H0 to be used in (3.5), we can write the decoupled equation for the δS < as follows: where the collision term in (3.12) is written in terms of full S < including the vanishing adiabatic part, and the source term S < ad can be written as: To get this form we used the adiabatic pole equation for the spectral function along with the result (3.11).We remind that the approximation used for Σ H does not need to be the same as the one used for the adiabatic self-energy function Σ H,ad and moreover, the convolutions associated with the self-energies can be computed also to higher order in gradients.Indeed, note that the last two terms cancel to the lowest order in gradients.
The gradient source i 2 ∂ /S < 0,ad is relevant in applications with rapidly changing backgrounds, such as resonant leptogenesis [53].On the other hand, in the transport equations for the electroweak baryogenesis the entire CP-violating source term comes from the two last terms in equation (3.13) [54].However, both terms can be dropped in very slowly varying backgrounds which is usually the case in light neutrino physics.The term involving Σ H on the other hand, will allow accounting for the forward scattering corrections to the evolution equations even when Σ H,ad was set to zero in the pole and the background equations.Conversely, if Σ H,ad = Σ H , then the self-energy corrections are already resummed to the quasiparticle dispersion relation and also the explicit self-energy corrections drop from the source term.We will return to this issue in section 6.3.
It will be useful to rewrite equation (3.12) in an alternative form, in terms of the full S < -function: The first term on the right hand side drops in the spectral limit, leaving no explicit source terms in (3.14).The gradient source and the sources arising from the different approximations imposed on the Hermitian self-energy functions in (3.12) remain however, hidden in the notation.

Local equations
The quintessential feature of the KB-equations (2.4) and (2.6), and of (3.12) and (3.14) is their non-locality.Our quest to reduce (3.14) to a local density matrix equation then clearly requires some further approximations.From the Wigner-space point of view the task is to curtail the infinite expansions in gradients.This can be justified by a further adiabaticity assumption, now concerning the perturbation δS < , or it can be enforced by integration over some of the momentum variables (encoding the lack of information on them).In what follows, we shall use both methods.

Adiabaticity in space coordinates
In essentially all problems relevant for neutrino physics, backgrounds are changing slowly in microscopic scales.This means that we can drop all spatial gradients acting on self-energies in the Wigner space evolution equations.In a generalization from the purely homogeneous case studied in [53], we keep the gradients acting on the correlation function however.To this end we Wigner transform equation (3.12)only over the spatial coordinates and then work in the adiabatic limit.The result is the following 2-time equation: where C< kx ≡ γ 0 C < kx γ 0 .The convolution is now only over time and we identified the vacuum Hamiltonian with α ≡ γ 0 γ and the correlation functions and self-energies with over-bars are defined as Ss ≡ iS s γ 0 , Sp ≡ S p γ 0 and Σs ≡ iγ 0 Σ s , Σp ≡ γ 0 Σ p .Also, to keep our notation compact and to highlight the essential features related to the time variable, we defined a shorthand notation S< kx (t 1 , t 2 ) ≡ S< (k, x; t 1 , t 2 ).Finally, we dropped the explicit source term appearing in (3.14), because we are eventually going to the spectral limit with the adiabatic solutions.This source could be simply added at any point of the derivation however.

Localization in time
Our next step is to localize equation (4.1) in time.From the 2-time perspective the motivation for this is obvious as discussed above.From the Wigner space point of view the localization corresponds to an integration over the frequency variable 3 .In other words, localization corresponds to working to the lowest order in a moment expansion in the frequency variable in Wigner space.We will discuss localization from a more general point of view in section 9, relating it to a definite statement of the prior information on the system.
Deriving local equation is simple.We start by taking the total derivative of S < kx (t, t) and using the chain rule to get In the second step we used (2.7), dropped the total derivatives under the integral and used the adiabatic assumption to replace Σout → Σ.From the 2-time perspective the problem is that the convolution depends on the unknown δ Skx (w 0 , t) for arbitrary w 0 , while equation (4.5)only yields the solution for w 0 = t.The loss of closure is a generic problem in deriving Boltzmann type equations from the SD-equations, and the usual solution is to reduce the Wigner space correlation functions to spectral limit [41,53].We will also go to spectral limit eventually, but for now we continue to develop the more general formalism including finite width corrections.To facilitate this we next introduce a novel homogeneous decomposition [53] for the perturbation δS < .

Homogeneous Ansatz
Remember that equation (4.5) is really an equation for δ S< , even though we wrote it in terms of the full S< for the simplicity of notation.Equations for δ S< admit a broader class of solutions than do the pole equations (3.1) and the equation (3.4) for the background.While the latter two admit only inhomogeneous solutions (in the sense of classifying the solutions to differential equations), equation (3.5) and its descendants have also homogeneous solutions which can always be written in the following form [53]: Homogeneous solutions are often used to model transients set up at some initial time t in .Indeed, the spectral function is the unitary time evolution operator in the free theory limit: , whose correct normalization is ensured by the spectral sum rule Here we follow the idea of [44][45][46]53] and use (4.7) as an ansatz for the non-local terms in the convolution integrals.This makes sense in dissipative systems with a finite damping rate γ k where only points with |u 0 − v 0 | < ∼ γ k are strongly correlated; for more discussion see [53].Remarkably, the structure (4.7) transforms all convolutions in (3.5) into simple matrix products: Here we used the spectral sum rule and we also adopted the more convenient notation setting e.g.δS < kx (t, t) ≡ δS < k (t, x) ≡ δS < k (x).The reverse ordered convolutions are just the Hermitian conjugates of the above:

The Master equation
The result (4.8) allows a tremendous simplification, reducing the original integro-differential equations to a set of ordinary differential equations.Indeed we can now immediately recast our local equation (4.5) simply as: where the forward scattering term is defined as and the Hermitian part of the local collision integral is We suppressed the x-dependence in all correlation and self-energy functions for clarity and we remind that we dropped the source term appearing in (3.14).Equation (4.9) is our final master equation.While the effective self-energies are still convolutions with the spectral function, they no longer depend explicitly on the perturbation.Also the convolution ΣH,k * S < 0k is just some known function.In the spectral limit it can be absorbed into the δ S< k -term in (4.10), which then allows writing . The master equation is formally closed, assuming that self-energy functions are defined externally.We define the Hermitian self-energy function in sections 5.3 and 8 and in section 7 we show that (4.7) allows expressing all self-energies including those in the collision integral in terms of the local correlation function.Note that no specific approximation scheme for the pole and the background solutions was needed in the above derivation.Different choices would lead to slightly different effective self-energy functions without changing the form of the master equation.

Density matrix equations
We now proceed to derive the quantum kinetic (density matrix) equations from the master equation (4.9).These equations take the simplest form in a projective representation onto the helicity and the Hamiltonian eigenbases.This classifies solutions according to their eigenfrequencies which simplifies equations and helps the analysis and interpretation.We will work explicitly in the vacuum Hamiltonian basis, as this is sufficient in most, if not all problems in neutrino physics.A generalization to the quasiparticle basis incorporating a resummation of thermal corrections to the Hamiltonian is discussed in section 6.3.

Vacuum representation
In a homogeneous and isotropic system helicity is conserved and the whole Dirac algebra is spanned by eight primitive structures that we list in (A.1).However, these structures can be more conveniently chosen [47,53] in the following combinations: where the helicity 4 and the vacuum energy projection operators are defined as: Here h = ±1 is the helicity, a, b = ±1 are the energy sign indices and is the vacuum energy of the neutrino eigenstate.Projection operators satisfy completeness relations P + ki +P − ki = P k+ +P k− = ½, the orthogonality and idempotence relations P a ki P b ki = δ ab P a ki and P kh P kh ′ = δ hh ′ P kh , and the eigenequations H ki P a ki = P a ki H ki = aω ki P a ki as well as α • kγ 5 P kh = P kh α • kγ 5 = hP a ki .The normalization factors N ab kij are most conveniently chosen as follows: This construction generalizes directly to the adiabatic case, so that we can parametrize our Wightman functions without any loss of generality as follows: where f <aa ′ khij (t, x) are some yet unspecified functions.Given the normalization (5.3) it is easy to show that Tr[ S< kij P e ′ e khji ] = f <ee ′ khij .This normalization also simplifies the dynamical equations and leads to the standard normalization of the distribution functions in the thermal limit.

Projected master equation
Using the projective representation (5.3) we can easily derive a generalized density matrix form for the master equation (4.9).We insert (5.3) into (4.9),multiply with P e ′ e khji (note the order of the flavour and the energy sign indices) and take a trace over the Dirac indices to extract scalar equations for the eigenfunctions f <ee ′ khij (t, x) 5 : 4 The non-covariant operator P kh corresponds to helicity only in the rest frame of the energy eigenstate.Otherwise it measures the spin along the momentum in a given frame.The covariant helicity operator be used, because its definition would require an exact definition of the energy shell, and this information is not available to us.Indeed, if it were, then each state would be defined precisely and no oscillations would take place. 5We implicitly assumed that in the forward scattering terms in (5.5) the background solution has the same form as the perturbation.To be consistent with our most general assumptions, one should replace e.g.
khji ( ΣH ⊗ S < 0 ) .While the trace-term is a known function, it has a different form than the first term, except in the spectral limit.Remember that we also dropped the source term appearing in (3.14), which would add a further term: where k = k/|k|, a sum over the repeated indices a and l is understood and e defined the oscillation frequency: with ω e ki ≡ eω ki .The forward scattering coefficient tensor is given by: khji ΣH effkil P ae ′ khlj . (5.7) In deriving (5.5) we used (f <ae khil ) * = f <ea khli , which follows from the Hermiticity of S< k (t, t), as well as the orthogonality relation Tr[P aa ′ khij P e ′ e khji ] = δ a ′ e ′ δ ae and we defined where the velocity tensor then is with and v ki ≡ |k|/ω ki .The master equation (5.5) is written in terms of frequency states rather than particle and antiparticle solutions.The positive frequency solutions directly correspond to particles of course, while the antiparticles correspond to the negative frequency solutions with inverted 3-momenta.At the level of distribution functions this implies the following relation: where functions in the l.h.s.refer to antiparticles.It is indeed notationally much simple to work with frequencies without an explicit identification of antiparticles at the level of the evolution equation (5.5).One can always convert the initial conditions and the final results to the particle-antiparticle language using (5.11), however.
The first term on the right hand side of (5.5) comes from the commutator with H k and it induces the leading time-dependence of solutions according to (5.6).The left-hand side of (5.5) displays a modified Liouville term where the tensor (V e ′ e khij ) aa ′ encodes the effect of different group velocities on the coherence evolution.The interaction terms are cleanly separated into a collision integral and the forward scattering terms.All terms in (5.5) thus have a clear physical meaning and the apparent complexity of the equation merely ad,kh S H,ad,kh )ij + h.c.)P e ′ e khji ] to the r.h.s. of (5.5).Instead of writing equation (5.5) in a complete but cumbersome form for δf <ee ′ khij , we prefer the simpler, slightly inaccurate notation, keeping these caveats in mind.These issues are eventually not relevant for the neutrino physics applications where the spectral limit can be taken. 6The minus sign in this definition is a convention that arises from our definition of the spectral representation (5.4).Indeed, from S> + S< = 2 Ā we get For the negative frequencies however, reflects its wide generality; equation (5.5) describes both flavour and particle-antiparticle oscillations for arbitrary neutrino masses with arbitrary interactions in backgrounds that are only constrained to be adiabatic in space, with a large freedom in the choice of the adiabatic approximation scheme.
To proceed we must finally specify our approximation scheme(s) and define how to compute the self-energy functions and collision integrals.We will first consider the simplest approximation, i.e. the (vacuum) spectral limit.After that we will define forward scattering terms involving known Hermitian self-energy structures.The definition and evaluation of generic self-energies and the collision integral will be discussed in section 7.

Spectral limit
Spectral solutions are easy to find directly using the projective parametrization (5.4).One can start with the Hermitian part of the statistical KB-equation (2.6) in the collisionless limit to the zeroth order in gradients: (5.12) Inserting here the equivalent of (5.4): . This is a spectral equation whose solutions are distributions where f see ′ khij are some shell-functions that parametrize the correlation functions S s il (k, x).The cQPA shell solutions (5.13) were found and used to derive QKEs for fermions and bosons in the spectral limit in [41,[44][45][46][47].The frequencies ωee ′ kij are given by ωee ′ kij ≡ 1 2 (ω e ki + ω e ′ khj ). (5.14) We display these shells for two-neutrino mixing in figure 1.A solution corresponding to a spectral shell k 0 = ωee ′ kij has the oscillation frequency 2∆ω ee ′ kij , given by (5.6).The particleantiparticle coherence solutions with e = e ′ reside at k 0 ≈ 0 and oscillate very rapidly in comparison to flavour coherence solutions with e = e ′ but i = j that form tight bundles with the usual mass-shell solutions with e = e ′ and i = j.We will use the vast difference in the oscillation frequencies between the particle-antiparticle and flavour mixing to derive a much simpler evolution equation limited to only flavour mixing in section 6.1 below.Homogeneous Ansatz with the spectral limit The Ansatz (4.7) is a more general construction than the spectral cQPA solution, but it reduces to (5.13) when one uses spectral free theory solutions for the pole functions.The free spectral function in the Wigner space is given by: Here it is simpler to use the direct space representation of this function, which is just the time-evolution operator: (5.16) where s =<, > and δ f sab khij ≡ exp(2i∆ω ab kij t)δf sab khij and we defined with Moving to the Wigner space in frequency and assuming t = 1 2 (u 0 +v 0 ), equation (5.17) becomes (now written for δS s ij without the bar) This is the just the spectral cQPA-result (5.13).We shall see in section 7 that the nonvanishing phase factors for t = 1 2 (u 0 +v 0 ) have a crucial role in ensuring correct 4-momentum conservation over the internal vertices in loops contributing to collision integrals [65].
In the diagonal limit a = b and i = j the structure D ab khij reduces to the standard form: We then get: This solution has the same form as the most general spectral adiabatic background solution depending on flavour, helicity and frequency, proving that in the spectral limit the background solutions can indeed be merged into diagonal transient solutions.One can then write the full solution in the same form as (5.21): where f sab khij = δ ab δ ij f sa eq,kii + δf sab khij .We stress that the spectral limit for adiabatic solutions is convenient and often sufficient, but not obligatory assumption for the analysis of (5.5) and the self-energy terms that appear in the equation.One could add non-trivial gradient corrections and finite widths to (5.5) and to spectral function (5.15), although this would come with a considerable amount of additional tedium.Our goal has been to balance between the full generality and the simplicity of notation for the benefit of both approaches.

Hermitian self-energy corrections
In essentially all problems in neutrino physics one can neglect dispersive and finite width corrections to background solutions and in the definition of the projective basis.However, we have included forward scattering corrections to (5.5).As we discussed at the end of the section 3, we can use different Hermitian self-energy functions in the pole-and the background equations and in equation (5.5).This freedom induces forward scattering corrections to (5.5) even when the pole equations are treated in the free theory limit.
Moreover, we can evaluate the effective Hermitian self-energies (4.8) using the vacuum spectral function (5.15).This immediately reduces the effective self-energies into simple products: The forward scattering coefficients then become khji ΣHkil (aω kl )P ae ′ khlj . (5.23) This expression involves ordinary self-energy function instead of the effective one.One should note that the self-energy function inherits its energy sign from the nearest energy projector to the right from the self-energy function.
In general ΣHk can depend on the dynamical quantities we are set out to solve.An example is given by the first diagram in the figure 2 which contains an internal neutrino line.This correction can cause a strong back-reaction from local neutrino densities, which have been shown give rise to interesting new phenomena in the early universe [32,[36][37][38][39] as well as in the core collapse supernovae and in the accretion discs around compact objects [52].Often Σ Hk is dominated by the interactions with the background however, and in such cases one can use some adiabatic (thermal of finite density) approximation for Σ Hk .We give a simpler treatment for such cases below and postpone a direct evaluation of the diagrams involving neutrinos to section 8.
Simple 1-loop weak gauge corrections.Any homogeneous self-energy function can be expressed in terms of structures (A.1) and evaluated in the projective basis using table 1.The weak gauge interactions in particular induce a one-loop self-energy for light neutrinos corresponding to diagrams shown in fig. 2. This self-energy can be written in the following general form: where a ij (k, x) and b ij (k, x) are some space-time varying flavor matrices and u µ is the plasma 4-velocity.In the second line we went to the plasma rest frame u µ = (1, 0) and in the third line we used the fact that operators α • kP L and −hP L have the same projective representation in the vacuum eigenbasis.In this case the forward scattering coefficient tensor becomes where V khij (k 0 , x) is the matter potential experienced by propagating neutrinos, defined above in (5.24) and the tensor σ abc klji ≡ Tr[P ca khil P L P bc khji ] is given by (see (A.1)): where ) is yet another useful projector function.In the ultra relativistic (UR) limit P a khi → δ a,−h .Even with thermal assumption, the expression (5.26) is much more general than the ones presented in the literature so far, see e.g.[48].It will simplify significantly in the frequency diagonal and in the UR-limits to be considered below.
Renormalization.Let us briefly discuss the role of renormalization in our quantum transport formalism.The key point is that renormalization only affects the vacuum parts of the correlation functions arising from a need to regulate and properly define the vacuum selfenergy corrections.From the results of section 3 it is then clear that renormalization only affects the vacuum limit of the pole and the background statistical functions.Moreover, the SD-equations for these functions decouple in the spectral limit, after which renormalization procedure only concerns the Hermitian self-energy function and can be done using the standard field theory methods, augmented to account for the flavour mixing.We refer the reader for example to [53] for a recent explanation of the procedure.For our current purposes, we can assume that our Hermitian self-energy functions are properly renormalized and we are using renormalized masses and couplings to parametrize our theory.

Limiting cases and extensions
So far our results are very general and in particular valid for arbitrary masses.This results in complex tensor structures in the projected equations and in the projected self-energies.In some problems, such as the production of flavour mixing particles from background fields, this complexity is unavoidable.Considerable simplifications arise if one can neglect the particle-antiparticle mixing and even more if one can take the UR-limit.Indeed, if one was only interested in the flavour oscillations, the particle-antiparticle mixing is but an unnecessary complication.We now show how to remove these structures from (5.5) by a simple integration [53].

Integrating out the particle-antiparticle oscillations
As was already pointed out, the particle-antiparticle oscillations have very high frequencies and they should average out in the usually much slower flavour oscillation scale.Indeed from equations (5.5) and (5.6) the leading time-dependence of the particle-antiparticle coherence functions is f e−e khij (t) ∼ exp(−2ieω kij t) with 2ω kij = ω ki + ω kj , while for the flavour mixing functions f ee khij (t) ∼ exp(−2ie∆ω kij t) with 2∆ω kij = ω ki − ω kj .This suggests [53] to take a Weierstrass transform of equation (5.5): where Then all terms in the equation which vary in the flavor scale are essentially unchanged by the transform, while the terms proportional to the coherence functions get exponentially suppressed: where c e kh stands for any coefficient of δf e−e kh in the equation of motion which is assumed to vary only in the flavour scale.The Weierstrass transform thus induces a coarse-graining with a temporal resolution scale σ which effectively washes out the particle-antiparticle mixing from the master equation.
In the absence of the particle-antiparticle mixing the tensor structures in (5.5) simplify considerably.In particular the velocity tensor in the spatial gradient term becomes just where vkij is the average velocity of the flavour states i and j with momentum k.The projected density matrix equation restricted either to particle or antiparticle sector then becomes: No sum over the energy signs remains here, but the sum over the repeated flavour index l persists.The tensor σ eee khlji still has a complex mass-dependence for general kinematics.We continue to defer the evaluation of the collision integral to a later stage.Here we also made explicitly the Feynman-Stueckelberg interpretation and replaced everywhere: After this identification f + kh refers to particle and f − kh to antiparticle flavour density matrix.Equation (6.4) is relevant for studying for example the resonant leptogenesis.It presents a generalization from [53] in that it gives explicit expressions for the dispersive corrections that so far have not been included in any leptogenesis calculation.Because (6.4) also allows for (adiabatic) evolution in the spatial coordinate, it could be easily used to accurately model the heavy neutrino production and the associated lepton number violating processes in the collider experiments with arbitrary heavy state kinematics.

The ultra-relativistic limit
Equation (6.4) simplifies further in the ultra-relativistic limit.UR-limit can be always taken when dealing with light neutrinos and it is therefore useful to show the equation explicitly in this case.Using UR-expansion for energies it is easy to show that With this result we immediately get: The forward scattering terms now collapsed to the familiar light neutrino matter potentials.
If we further define an effective matter Hamiltonian: and diagonal velocity matrix v kij ≡ δ ij |k|/ω ki , we can write equation (6.4) in the UR-limit in the compact, familiar form of a density matrix evolution equation: where ( Ce kh ) ij ≡ Tr[ C< H,khij P ee khji ].When dealing with light neutrinos over relatively small propagation distances, one can further set v kij → δ ij , in which case the spatial flow term reduces to 1  2 {v k , k • ∇f e kh } → k • ∇f e kh .Equation (6.9) looks deceivingly simple, but it still contains all information of flavour coherences and forward scattering potentials in the UR-limit.It also has the same form as early UR-limit kinetic equations derived by the S-matrix or operator formalism techniques [27][28][29][30] and it is sufficient for almost all light-neutrino physics applications, from laboratory experiments to light neutrino interactions in the early universe and within high density astrophysical objects.Before we turn to crucial issue of computation of the collision integrals we still discuss two issues related to the matter Hamiltonian and the choice of the basis one uses for the pole funcctions.
Flavour and the matter eigenbases So far we have worked exclusively in the vacuum basis.However, given the effective Hamiltonian, we can perform rotations to the flavour or matter bases in the usual manner.In particular the transformation between the flavor and mass basis is just a constant rotation.
where in the standard case with three light neutrinos U would be the usual PMNS mixing matrix.A similar transformation could be performed to go to the matter basis where H e kh is diagonal.The matter basis, which always depends on k and possibly on h, in general also rotates along the neutrino path, because V e kh = V e kh (x).As a result U m ≡ U e kh (x) depends on the space-time coordinate x as well, which gives rise to the additional Liouville terms after the rotation into the matter basis:

Quasiparticle basis
Until now we have used vacuum solutions both for the projective representation and in the reduction of the effective self-energy functions.This is a very good assumption in all light neutrino physics applications and in the leptogenesis problem.However, in some cases dispersive corrections can change the phase space structure significantly.A well known example is the infrared region of a thermal plasma with gauge interactions, where new collective hole excitations appear [66][67][68].This situation is realized in some electroweak baryogenesis scenarios [69][70][71][72] due to strong flavour blind QCD interactions.For illustration we briefly consider this case where additional weak flavour mixing interactions can be treated perturbatively as discussed above.
The QCD-interactions are vector-like and induce a Hermitian self-energy correction, which in the hard thermal loop (HTL) approximation has the form Because Σ th H is flavour diagonal, we can work in the vacuum mass eigenbasis.Neglecting for the moment all other corrections, the inverse propagator for the system can be written as S −1 = rnp 0 γ 0 − rγ • k − m i , where r ≡ 1 − a and rnk 0 ≡ rk 0 − b.This propagator has two branches of poles given by: 12) The positive sign corresponds to particle and the negative sign to the hole solutions of [67,68].One can now derive the effective Hamiltonian near these quasiparticle shells [54]: Thermal corrections are manifested mainly via the nontrivial refractive index n i .One can proceed to construct the energy projectors using (6.13) in entirely analogous manner to the vacuum case.One should also include the thermal wave-function corrections for the quasistates, see e.g.[54].Similarly, one can use quasi-particle generalization of the spectral function in the evaluation of the various self-energy functions.Finally the perturbative flavour changing interactions can be added on top of this structure similarly to the vacuum case.
An even simpler generalization to the energy-shell projectors, relevant for the leptogenesis problem, is to replace the constant masses by a time dependent masses.In this case the energy projectors pick up a time dependence through masses that gives rise to additional terms proportional to ∂ t m i in the projected equations (5.5), (6.4) and (6.9).For explicit expressions see [53].

Collision integrals
In this section we show how to compute the collision integrals with coherent states in our transport equations.Our approach is different from section 5.3 where we assumed generic structures for the Hermitian self-energy function and worked out their projections.Here we construct the expansion of the collision integrals directly in the projected basis.That is we will compute the collision integral traces: where in the last step we used the orthogonality of the energy projection operators and the normalization (5.3).We now insert the expression (4.11) for C< H,khij into (7.1) and proceed similarly to what we did with the Hermitian self-energy function terms in section 5.1.We can then write the generic collision integral as follows where (W se ′ e khji ) l a -tensors are defined similarly to (W He ′ e khji ) l a -tensor in (5.23): for s =>, <.Here we preferred to wrote the projection operators in the D-tensor notation of equation (5.18).The key quantity in the expression (7.3) is the self-energy function: where again s =>, <.An essential element in the following analysis is the fact that the homogeneous Ansatz (4.7), which reduced temporal convolutions to simple products, also reduces an arbitrary function Σ sa kil , with any number of internal lines, computable in the local limit.From the reduction process we can infer simple rules for a diagrammatic construction of collision integrals.We give the main points of the derivation and the final results here.More details can be found in the companion paper [65].

General reduction process
Any diagram contributing to the self energy (7.4) consists of a number of vertices connected by propagators, with an integration over the spacetime coordinates in each internal vertex.Wigner transforming propagators then introduces a momentum integral for each internal line and a group of phase factors that after integration give rise to 4-momentum conservation over vertices.For the spatial coordinates and momenta this works out in the usual manner, but the time coordinate requires more attention.Moreover, the internal lines are divided to statistical (cut) and pole propagators, following the standard rules of the thermal field theory [46,65,73].Each cut propagator introduces a statistical f -factor that gets associated with an external state in the interaction process, while the pole propagators give rise to the internal resonances or generate loop corrections to collision rates.
It is simplest to work in the 2-time representation and we continue to assume the free theory limit for the pole functions.In this case the dynamical and the background solutions can be combined as explained in section 5.2.Following (5.19) the full statistical correlation function can then be written as: where s =<, >.In this article we consider only the gauge interactions and treat the gauge fields as non-coherent resonances.This means that we only need the 2-time representation of the standard gauge-field propagator: where D µν (q) is the usual momentum space propagator function.Extension to gauge fields, or additional scalar fields, that are a part of the non-equilibrium quantum plasma is straightforward, but to keep discussion simple we do not present it explicitly here.Figure 3 shows the propagator (7.5) associated with the internal time coordinates u 0 and v 0 in a generic self-energy diagram contributing to the self energy (7.4).Remarkably, the phase factors in (7.5) appear with different frequencies at coordinates u 0 and v 0 when either e = e ′ or i = j.This difference is crucial to ensure the correct energy conservation at both vertices even when the connecting propagator has support on a constant energy shell ωee ′ kij .Including also the phase factors from the gauge-field propagators (7.6), we see that the integral over u 0 results in 2πδ(aω pk − bω rl − q 0 ) and the integral over v 0 gives 2πδ(cω rm −dω p ′ n +q ′ 0 ).As usual, the delta-functions from vertices kill all frequency integrals associated with the propagators, leaving one extra delta-function that gives a generalized energy conservation for the process in question [65].Using the delta-functions one can also show that the sum of all explicit t-dependent phases, which appear in the definition (7.4) and in the cut propagators (7.5), vanishes exactly for any process [65].Instead of giving general proofs of these statements, we will show below how this works in specific examples.
Having gotten rid of all phase factors we see that each internal cut propagator effectively contributes the following factor to the diagram: (7.7) The differential fraction in the integral in the first line in (7.7) defines a natural phase space density factor and f sab khij (t, x) is the corresponding phase space distribution function.This leaves D ab khij as the sole contribution to the scattering matrix element.Alternatively, one can use the singular cQPA-form given on the second line as the Feynman rule for the cut propagators with the associated full four dimensional momentum integral.
These results can be summarized by the Feynman rules for evaluating self-energy diagrams shown in figure 4. In addition to these rules one should associate each propagator line with an integral over the four-momentum and sums over the helicity, the energy-sign and the flavour indices.For the pole propagators one should use just the usual momentum space Feynman rules.

Collision integral from two-loop gauge diagrams
In figure 5 we show all two-loop self-energy diagrams that give rise to the 2-2 scattering terms in collision integrals mediated by the weak gauge interactions.Only the four first Figure 4: The Feynman rules for the internal lines and weak interaction vertices associated with coherent neutrino propagators.In the W -boson vertex U iα reduces to the PMNSmatrix in the case of pure active-active mixing (α refers to lepton flavour).In the Z-boson vertex c w = cos θ w and the mixing matrix Ūij reduces to ½ for pure active-active mixing.
The rules generalize to arbitrary Lorentz and flavour structures in an obvious way.
diagrams are two-particle irreducible (2PI).The last two diagrams on the second row should not be included in the full SD approach, where they would be accounted for by the oneloop corrections to the gauge-boson equations.However, when gauge bosons are treated as non-dynamical resonances, the 2PI-hierarchy is partially broken and these diagrams must be added as perturbative corrections to gauge boson propagators.They eventually produce the squared matrix elements for the s, t and u-channels for a given scattering process while the 2PI-diagrams give rise to the interference terms between the different channels.
In addition to scattering terms these self-energies create a large number of 1-3 decay and inverse decay terms as well as vertex corrections to 1-2 decays.Different processes correspond to different cuts on the self-energy diagrams following standard rules of the finite temperature field theory [46,65,73].For illustration we show cuts that give rise to Z-mediated 2-2 scatterings and 1-3 decays between neutrinos and other fermions in the first and fifth diagrams and a cut producing W -boson 1-2 decay correction in the third diagram.The 2-2 scatterings and the 1-3 decay corrections are further distinguished by the frequency sign signatures in the overall momentum conservation function following the standard kinematic analysis which we shall elaborate more in [65].
The self energy function Σ <a ZZ,il We now work out the self-energy function Σ <a ZZ,il that describes scatterings between coherent neutrinos mediated by Z-gauge bosons.The relevant diagrams are the first in the upper and the second in the lower row of fig. 5. We show these diagrams again in figure 6 including labels for the momenta and the discrete indices for the intermediate states.To alleviate the complexity of notation we will be using the shorthand Note that the rightmost propagators marked red in figure 6 are not part of the self-energy function, but they do contribute to the collision integral according to (7.3).The light blue numbers at each vertex correspond to the Keldysh time-path indices [65,73].They allow where p ≡ d 4 p/(2π) 4 and the curly brackets indicate groups of indices to be summed or variables to be integrated over.It is easy to perform the q 1 and q 2 -integrals using the delta-functions from vertices, which leaves out the one overall momentum conserving delta function.After performing the integrations one finds: where the phase space integral is defined as: Figure 6: Two-loop graphs that contribute to the matrix element squared of neutrinoneutrino scattering through s and t channels are shown with explicit index structures and cuts.The black dot implies the starting point of the evaluation of the trace.The red propagator is the dependent momentum propagator, see section 7.3 for explanation.
where the energy sign indices tell whether a given term contributes to a 2-2 scattering or a 1-3 decay channel.The distribution functions in (7.9) were gathered into the factor and the part that eventually contributes to the matrix element was defined as: where the gauge-boson 4-momenta are q 1 = (ω a kl −ω ; k−p 3 ).Evaluation of the "direct" diagram shown in the right in figure 6 proceeds analogously.The phase space-integral (7.11) and the element (7.12) containing the distribution functions are the same as in the interference diagram.The only difference stems from the different way of connecting the fermion lines and from the value of the gauge boson momentum q 2 , which give rise to a different A-factor: where Note that both 3-momenta in the direct channel (7.14) are the same q 1 = q 2 = k − p 1 , while in the interference channel q2 = k − p 3 is different.This is as expected, since direct diagram accounts for the s-and t-channels, while interference diagram accounts for their interference.Identifying different channels by Mandelstam variables is not obvious in presence of particle-antiparticle mixing and was done above based on 3-momenta.In the frequency diagonal limit the situation is simple.For example the t-and u-channel particleparticle collisions correspond to a = e = a i = a ′ i = 1 for all i, while the t-and u-channel particle-antiparticle collision correspond to a = e = a 3 = a ′ 3 = 1 and a 1 = a ′ 1 = a 2 = a ′ 2 = −1.These are just examples from the different processes embedded in the generic collision integral.Identifying and classifying all different channels is straightforward but somewhat tedious task of standard kinematic analysis.For more details see [65].
Collision integral for Z-mediated processes We have now collected all the pieces needed to write down the full collision integral corresponding diagrams in figure 6. Inserting the self-energies computed above to equation (7.1) and doing some simple reorganization we find where we collected all summed indices into curly brackets Y = {X i , a, a ′ , l}.The phase space factor is the same as in (7.11) and the factor containing all distribution functions is Finally, the effective matrix element squared is defined as with where the quantities in the right hand side are given in equations (7.13) and (7.14).Remember that f >ab phij = aδ ij δ ab − f >ab phij in general and the rule f <,> khij = −f >,<−− (−k)hij for translation between negative frequency and antiparticle distributions.Note that, as emphasized by our notation, the frequency and flavour indices get flipped in the Hermitian conjugate term in (7.17), as is clear from (7.2).
Let us stress that in the presence of mixing the collision integral components are in general complex; indeed, we found already in section 5.2 that a function f ab khij (t, x) has a leading phase −2∆ω ab khij t.From (7.16) and the phase space constraint in (7.11) it is then easy to work out that the leading phase of the phase factor term Λ (and hence that of the whole integrand) is This means that the integrand is a rapidly oscillating function in the particle-antiparticle mixing scale in all cases except the fully diagonal limit, where e = e ′ and a i = a ′ i for all i and in the fully anti-diagonal limit, where e = −e ′ and a i = −a ′ i .This implies that also the collision integral averages to zero under the Weierstrass coarse-graining of the equation of motion in the limit discussed in section 6.1.
Despite the compact notation, expression (7.15) is a very complex object that encompasses all flavour and particle-antiparticle mixing effects for arbitrary neutrino masses and kinematics.It still displays only the familiar structures from the usual Boltzmann theory and all complexity is reduced to a bookkeeping of the indices labeling the states.In particular, the phase-space factor and the dependence on the distribution functions are universal features independent of the structure of the interactions.The interaction matrix element is more complex than the usual result for non-mixing states, but its evaluation is formally straightforward and can be easily done using symbolical routines.We will evaluate (7.17) in the UR-limit in section 7.4 below.For now, we comment that the simple, intuitive form of our collision integral, which clearly separates the different flavour and frequency structures, appears to be in stark contrast with other existing QKE computations in the literature [49,50].

Simple Feynman rules for the matrix element squared
The simple factorization of the collision integral to universal phase space elements and a dynamical matrix element squared allows us to deduce a very simple set of rules to evaluate the collision integrals directly without the need to repeat the integration procedure for each new diagram and set of interactions.The rules we now spell out should be obvious from the above derivation.
• Draw the loop diagrams that contribute to a given interaction process to the desired order in perturbation theory, and assign unique momentum variable and flavor and frequency indices for each internal propagator line in the graph, allowed by the interaction vertices.
• Assign the Keldysh-path indices to all vertices to isolate the cut that gives rise to the desired interaction process.You only need to evaluate Σ > = Σ 21 directly, so the first index is always 2 and the last 1.
• Read off the phase space functions contributing to the Λ-factor from all internal cut propagator lines.Add the external phase space factor f <aa ′ khlj /2ω aa ′ khlj , which is associated with the external, dependent momentum propagator (DMP), marked red in diagrams in figure 6.
• Determine the phase space density factor with the overall energy conserving delta function.This depends on the number of loops in the diagram and the cut one is interested in.At two loops with a cut leading to 2-2 scatterings this was simple: we isolated the internal 11-propagator and eliminated its frequency and momentum from one of its end-vertex delta-functions using the other one.
• Compute the matrix element squared using the Feynman rules shown in fig. 7. Start from the equivalent of the black dot in the diagrams in figure 6 and follow the direction of momentum in the graph.For each internal cut-line insert the standard propagator shown in the first diagram in 7. Add the DMP at the end of the fermion line it is connected to.For each ("22") "11" line use the (anti) Feynman propagator.Add the DMP at the end of the fermion line it is connected to.Take a trace over the Dirac indices.
• Divide the result by two and add the Hermitian conjugate, accounting for the index changes as indicated in (7.17).
We invite the reader to verify that these rules indeed allow writing (7.15)-(7.17)directly from the diagrams shown in figure 6.These rules can also be extended, in an obvious way, to any other interaction types.One particularly interesting application is the Leptogenesis problem where the gauge interactions are replaced by scalar Yukawa interactions.In that case one also encounters Majorana neutrinos which require some special rules that can be found for example in [53].

Explicit results in the UR-limit
We return now to evaluate the effective matrix element (7.17) explicitly in the UR-limit.This limit is interesting for practical applications and because it allows to make a direct contact with some known results in the literature.We begin by pointing out some general simplifications: first the orthogonality of the energy projectors immediately sets a ′ = e ′ independent of the type of interactions.Second, in the case of neutrino-neutrino interactions the trace calculation simplifies due to the orthogonality of the chirality operators which immediately gives: where the last result applies in the UR-limit.To arrive at this result we used the expansion Similarly, one can show that in the UR-limit.From now on we will also assume that the Z-boson is very heavy in the energy scale of interest, which means that we can set: D Zµν = g µν /M 2 Z .With these simplifications the matrix element (7.17) becomes simple to evaluate.The result is: where G F is the Fermi constant and Ū 4,dir ilX ≡ Ūil 1 Ūl ′ 1 l Ūl ′ 2 l 3 Ūl ′ 3 l 2 and Ū 4,int ilX ≡ Ūil 3 Ūl ′ 3 l 2 Ūl ′ 2 l 1 Ūl ′ 1 l .Note that in the UR-limit all collision terms with particle-antiparticle mixing drop out at the leading order.The mixing terms would appear as O(m/E)-corrections only, which we have dropped above.Such terms also would have a large mixing phase.

Standard Model (SM) limit
To make equations even more transparent, we now assume the SM-limit where the neutral current rotation matrices are trivial: Ūij = δ i,j .Making use of the large number of Kronecker delta functions in (7.22), one can show that the general collision term (7.15) now reduces simply to F δ e,e ′ δ e,−h where the phase space factor is defined as before, but now the flavour indices are no longer necessary in the kinematic factors: and all flavour dependence is in the phase space factor: The first combination of f -functions in (7.25) arises from the direct diagram and the second set from the interference diagram.Note that while pure neutral current processes do not give rise to source terms in (5.5), nontrivial source and therefore nontrivial mixing contributions would be produced by charged currents.Equations (7.23) and (7.25) give the correct neutral current collision integral for such a setup.
Neutrino-antineutrino scattering Next consider the special case of ν ν − ν ν-scattering in the UR and SM-limits.To get this process we must set7 e = 1 = a 3 = 1 and a 1 = a 2 = −1 in (7.23)- (7.25).To impose the Feynman-Stueckelberg relations we first redefine p i → a i p i and k → ek.This immediately sets ) with ordinary dot-products.In the present case moreover ea 1 a 2 a 3 = 1.Similarly, the delta-function in the phase space factor now becomes δ 4 (k + p 1 − p 2 − p 3 ) appropriate for the process we are considering.Finally, we identify e.g.
. Working similarly with the other distribution functions, we find that the particle distribution factor becomes: and similarly for the antiparticle distributions, this falls into the familar form from the usual Boltzmann theory, except that the expression still contains all information of the flavour mixing.
If we finally take the flavour diagonal limit, all terms in Λ become real.We also see that the second (interference) term in (7.26) gives just one term with i = l = l ′′ = j.In the first (direct) term however, we only get i = l = j and the term still contains a sum over one flavour: l ′ = l ′′ .That is, when l ′ = i only the first term survives and gives the collision integral for the usual ν i νl − ν i νl -scattering coming from the single t-channel diagram.Indeed, the matrix element for this process is easy to compute and the result ) is in perfect agreement with the above results.For l ′ = i both terms contribute and reproduce the result for ν i νi − ν i νi -scattering obtained from summing the s-and t-channel terms using the usual field theory methods.
Two flavour active-sterile mixing Let us finally consider the case of two flavour mixing between an active (a) and a sterile (s) neutrino in the UR-limit.We will cast the collision integral for this system into the familiar form in the flavour basis.As before, we label the flavour states by Greek and vacuum basis states by Latin letters.In this case the rotation matrix between the flavour and vacuum basis U and the neutral current mixing matrix Ū are given by: where e.g.c ≡ cos θ, where θ is the vacuum mixing angle.Now observe that the rotation matrices Ū within Ū 4,dir ilX and Ū 4,int ilX always sandwich the distribution functions f <,> X i p i in equation (7.12).Because the matrix element function (7.17) does not depend on flavour in the UR-limit, the rotation amounts to replacing the vacuum basis matrices Ūij by their flavour basis representations Ūαβ as well as setting f sab phij → U αi f sab phij U † βj ≡ f sab phαβ everywhere.The final result is where we used (f <e khsa ) * = f <e khas in the Hermitian conjugate term and defined the real valued purely active rate: Here the direct term and interference term give equal contributions and are summed in the rate (7.29).As expected, the sterile state has no collision integral in the flavour basis and the off-diagonal terms are damped by rate that is half of the (sum of the sterile and) active rate(s).Further identification of the particle and antiparticle channels proceeds as in the previous example with neutrino-antineutrino scattering.
We have now reduced our initial collision integral with all flavour-and antiparticle mixing effects down to the flavour diagonal limit, where the contact to usual field theoretical methods was easy to make.We also extracted the known structure of the damping terms affecting the coherent flavour evolution in the two-flavour active-sterile mixing case in the UR-limit.We hope that showing these explicit results make it easier for the reader to understand how apply our methods to study any given problem at hand.

General forward scattering potential term
We now show how to evaluate the Hermitian self-energy diagrams directly at one-loop level.To be specific, we consider the leftmost bubble diagram in figure 2, with an internal neutrino line (we label it by ZB).From this example it should be evident how the other diagrams are computed.We begin by observing that the Hermitian self-energy is equal to the Hermitian part of the 11-component of the self energy in the CTP indices: ΣH = He( Σ11 ).We shall evaluate the self-energy and the corresponding forward-scattering coefficient (5.23) in the spectral limit, discussed in section 5.2.We also first assume that Ūij = δ ij .It then is easy to see that the ZB-diagram gives: To obtain the most general result, we should derive a non-equilibrium CTP-propagator also for the gauge boson.In most cases of interest however, the gauge bosons appear only as intermediate resonances, which allows to replace D 11 Zµν with the standard vacuum Feynman propagator.For the S 11 -function we use the identity S 11 = S r − S < = S H − iA − S < .Using the spectral results (5.15) and (5.21) we then get where we combined part of the spectral function with the Hermitian propagator S H (the principal value propagator in the spectral limit) to extract the vacuum propagator term.It is easy to see that (8.2) reduces to the standard thermal Keldysh propagator in flavour diagonal thermal limit.The vacuum contribution from iS 11 il is absorbed by the standard renormalization procedure and we only need to consider the second term in (8.2).Assuming that the energy scales in the problem are small compared to the gauge-boson mass, we can further set iD 11  Zµν ≈ ig µν /M 2 Z .In this (tadpole) limit Σ ZB Hil is actually independent of k, and we can perform the q 0 integral trivially to get To get the forward scattering coefficient one inserts this function to the trace in the formula (5.23).Equation (8.3) is a very general result however, and to keep the discussion simple, we now specialize to the UR-limit, which is also frequency diagonal.Employing the form (7.3), using results (7.20) and (7.21) and computing the resulting simple trace term, one finally gets Since the expression in the right hand side of (8.4) is independent of j, we can define, similarly to (6.8), (W ZB,Hee

khij
) l e ≡ (V ZB,e kh ) il .Setting q → bq and using the Feynman-Stueckelberg relation (5.11) along with the equality f If the background is isotropic, then the directional term proportional to q • k vanishes, and we recover the familiar expression [27,30], written in the vacuum mass eigenbases.In the supernova application, the directional term cannot be neglected however, and the complete structure shown in (8.5) should be used.
It is easy to generalize the above calculation for a general mixing matrix Ūij .The result is simply that (V ZB,e kh ) il → ( Ū V ZB,e kh Ū ) il .This can be further rotated back to the flavour basis using the rotation matrix U iα .For example, in the two-flavour active-sterile mixing case discussed above, the flavour space effective potential gets the expected form: where (V ZB,e kh ) aa is as in equation (8.5) with f < q−aa = c 2 f < q−11 + s 2 f < q−22 + cs(f < q−21 + f < q−12 ) and similarly for the antiparticle term.Other diagrams can be computed similarly.In particular the tadpole diagram contributes a term (V ZT,e kh ) il = Ūil Tr[ Ū (V ZB,e kh )] with the same set of approximations.In the above 2-neutrino active-sterile mixing case this leads to the same result as for the ZB-diagram: (V ZT,e kh ) αβ = (V ZB,e kh ) αβ , which is again the expected result.

Weight functions
A central element in our reduction of the non-local Kadanoff-Baym equations to a density matrix equation, was temporal localization, or from the Wigner space point of view, the integration over the frequency variable.The resulting loss of closure required new assumptions about the full correlation function, elaborated in sections 4.3 and 5.2.While this approach is clearly successful, one might ask if and to what extend it was unique?It is indeed not, albeit the final equations often do not depend on precise details.The issue is related to how the prior information, or preparation of the system affects its evolution.
First consider the exact solutions to the KB-equation (2.6) at formal level.Qualitatively we know that, due to gradient corrections and finite widths as well as flavour and particle-antiparticle mixing, they manifest some intricate structures, which in general are well localized in the phase space.We have solved these structures explicitly in the spectral limit in section 5.2.However, we can never have complete information of any system we observe or describe, and sometimes the resolution of the setup is too poor to discern certain individual structures.That is, we never have access to but some coarse-grained version of the actual system.
Since by the "system" we basically mean the correlation function, we can formally express the relation between the exact and the observable systems as follows: where W is the weight function encoding the observational resolution.The weight functions can affect any or all variables relevant for the problem.The parameters relevant for this paper were helicity, frequency, 3-momentum, and spatial and temporal coordinates.In the electroweak baryogenesis problem one would be particularly interested in the momentum, the spatial coordinate and the spin perpendicular to the phase transition wall.We only display the weight functions for continuous variables for brevity.
In this language taking the local limit corresponds maximal coarse graining in the frequency variable, or to a statement that there is no prior information about the frequencies.This can be formally described by the weight function This setup is suitable for studying problems including particle-antiparticle mixing, such as particle production, where the particle-and antiparticle mixing shells are widely separated.
However, in other problems some different weight functions might be more appropriate.The weight function (9.2) is very simple, consisting of a flat distribution in frequency and delta-functions in 3-momentum and space-time coordinates.These are idealizations of more general functions.The flat weight could, to the same effect, be replaced by a very broad and the delta-functions by very narrow Gaussian distributions.In some cases it can be useful to use such weight functions to enter more detailed prior information on the setup.Indeed, we have seen an example of this already in section 6.1, where we integrated out the particle-antiparticle mixing from the general projected master equation (5.5).For the parameters we used this procedure is effectively equivalent to using the weight function8 3) The same effect could be obtained also by a Wigner space weight function does not have the resolution to erase the flavour structures but erases the information from the particleantiparticle mixing.For example: Given 2∆ω ab kij ≪ σ k ≪ ωab kij this weight function picks the frequency diagonal particle-and antiparticle solutions of all flavours, but suppresses all particle-antiparticle mixing solutions around k 0 = 2∆ω ab kij .Of course the use of a particular weight function should be motivated by physical arguments, but it is not difficult to see how such motivation would arise in specific setups.Consider for example a laboratory experiment with a neutrino beam.One can usually determine very well what particle or antiparticle flavour states are produced, but their energy and momentum resolution is insufficient to resolve the emitted mass eigenstates.This lack of knowledge imposes the need to integrate exact QKE's over some phase space patch with a weight function, possibly of the form (9.4), with parameters reflecting the experimental resolution.When the experimental resolution is poor compared to spacing of flavour structures, this procedure would lead to our equation (6.4), independent of the precise structure of the weight function.
These are just some simple examples of how the prior information imposed on system by integration with a weight function affects the resulting evolution equations.More general weight functions could turn out to be useful tools to study quantitatively how the specific details of neutrino production and detection processes, e.g. the observer-system interference, affects the evolution of the system.

Discussion and conclusions
In this article we derived quantum kinetic equations (QKEs) for neutrinos that encompass both the flavour and the particle-antiparticle mixing.Our results include explicit forward scattering terms and collision integrals for the coherent neutrino states.We started from the most general Kadanoff-Baym equations and reduced them, in a set of well defined steps, into a single local density matrix equation (5.5), which is valid for arbitrary neutrino masses and kinematics and contains all flavour and particle-antiparticle coherence effects.To our knowledge evolution equations of this generality have not been presented before.We then showed how to consistently integrate out the particle-antiparticle mixing and derived separate (but coupled) equations for particles and antiparticles (6.4) that are still valid for all kinematic variables.Finally we took the ultra-relativistic (UR) limit of these equations recovering the familiar form of a density matrix equation (6.9).
Our analysis is closely related to earlier work done in refs.[40][41][42][43][44][45][46][47], but extending it in many ways.Pivotal elements in our derivation were the careful separation of the pole-and statistical KB-equations and the introduction of the projective representation (5.4).This unveiled a novel shell structure underlying the mixing phenomenon and allowed expressing the master equation as a set of scalar-valued equations for distribution functions classified according to well defined oscillation frequencies.The only physical assumptions made during the derivation were the slowly (adiabatically) varying background fields, the validity of the weak coupling expansion and eventually the spectral limit.
The definition of the collision integrals with all information of the flavour and particleantiparticle coherences is a very delicate problem, whose importance has been recently emphasized [52].Although some computations in the absence of the particle-antiparticle mixing exist [51,[74][75][76][77][78] and even some others that do include them [45][46][47], a complete and comprehensive treatment of the forward scattering terms and collision integrals with the most general mixing structure and arbitrary neutrino masses and kinematics has not existed so far.In ref. [50] a formulation of the collision integral is presented in the relativistic limit, but it does not seem convenient for practical purposes.In contradistinction, our derivation includes a simple set of Feynman rules which provide a straightforward and systematic way to compute collision integrals for the flavor-and the particle-antiparticle mixing systems.We showcased the simplicity of our formalism by deriving several simple and/or known limiting cases for the collision integrals.In particular we identified the damping terms in the two-flavour active-sterile mixing case in the UR-limit.
In addition to explicit collision integrals, we also showed how to compute the forward scattering corrections coming from diagrams involving coherent states, i.e. the forward scattering effect arising from a coherent neutrino background.Again we first showed how to obtain the most general structures with both flavour and particle-antiparticle coherences and then took the frequency-diagonal UR-limit of this result, which eventually revealed the familiar structures found using other, less fundamental approaches.
Finally, we briefly discussed how the prior information on the system impacts on its evolution.We pointed out that our localization procedure corresponds to a complete ignorance of the frequency structure of the correlation function and that it is precisely this lack of information that allows for the non-trivial oscillation structure to emerge.However, the perfect localization is an idealization and we sketched how more detailed information of the system could be imposed by the use of specific weight functions.This seems like a promising way to study the observer-system interference in general, beyond the astrophysical and early universe applications.
It has been speculated that neutrino-antineutrino coherence could be relevant in some astrophysical environments or in experiments involving the decay of heavy neutrinos [52,55], but we believe that this is not the case.Because due to their very fast rate the particleantiparticle oscillations average out completely in the time scales of interest.An exception to the rule could be the leptogenesis problem in the non-resonant regime.In contrast, the particle antiparticle mixing is essential for the particle production problem, e.g. during the (p)reheating phase after inflation, and our most general QKE's (6.4) are the right tool for studying this problem in the presence of flavour mixing.However, in the astrophysics applications the relevant physics is the CP-violating flavour mixing in the separate, but possibly strongly coupled (via the forward scattering term) particle and antiparticle sectors.For these systems, that include the core collapse supernovae, nascent neutron stars, compact object mergers and the primordial nucleosynthesis, our frequency diagonal equations (6.4), and quite often their UR-limit (6.9), are sufficient.
The same critique applies to some investigations of the CP-violating decays of heavy Majorana neutrinos in collider experiments.Indeed, there has been some interest in the possibility of measuring the heavy neutrino mixing parameters in colliders, because this could give insight into the neutrino mass generation mechanism and eventually into the question if the baryon asymmetry could be explained by the leptogenesis mechanism.To be specific, it has been argued [55] that the heavy neutrino-antineutrino oscillations could lead to oscillation in the lepton number conserving (LNC) and the lepton number violating (LNV) decay rates, which could be testable at the LHC.Again, the physics behind the phenomenon is not particle-antiparticle coherence, but the CP-violating flavour mixing of the heavy neutrinos.For these systems our frequency diagonal equations (6.4) are the ideal tool to use, valid for arbitrary neutrino masses and kinematics.
Our QKE's were derived with only a very few approximations, but there are still some limitations to their applicability.In the local limit, which corresponds to working to the lowest order in the frequency moment expansion in the Wigner space, all truly non-local effects, such as quantum entanglement, are lost.However, we have a rather clear idea about the size of the neglected effects, as they are all encoded in the gradient expansion in the Wigner space KB-equations.If the adiabatic limit in the spatial variables is a reasonable approximation, which often is the case, then all such effects are very small.This type of non-locality could be relevant in some early universe applications however, and it has been studied in the context of the QKE's e.g. in [47].Also, when computing the collision integrals and the forward scattering terms we eventually assumed the spectral limit, which neglects the finite width corrections.For the light neutrinos this should be an excellent approximation, but in the case of heavy unstable neutrinos the finite width effects could be interesting.We tried to address this issue by structuring our derivation such that the extension of our analysis to include finite width corrections should be at least in principle evident and it would be interesting to study these issues more in future.Our results should be useful in the study of many interesting problems related to astrophysics, particle physics and physics of the early universe.k, the projection can be written in the form:  where we defined the velocity v ki ≡ |k|/ω ki .

Figure 1 :
Figure1: Shown is the cQPA-shell structure for two-neutrino mixing.(The continuation to negative k is a simplification of the full 3-dimensional rotation symmetry.)Blue and red lines denote the mass shells, purple lines show the flavor coherence shells, and green dashed lines are the particle-antiparticle coherence shells.Black ellipses illustrate the uncertainty on momentum and frequency in the preparation of the system, which would prevent determining exact flavour shells (mass-eigenstates), but would eliminate the particle-antiparticle mixing.The black line on the negative momentum side illustrates the use of a flat weight in frequency and ideal accuracy in momentum (see section 9). ν

Figure 2 :
Figure 2: One-loop diagrams which contribute to light neutrino self-energy.

Figure 3 :
Figure 3: A section of a generic diagram contributing to the projected local collision integral (7.1).The dashed outline shows the extent of the self-energy contribution and the detail at the center displays insertion of nontrivial coherent propagators (7.5).

Figure 7 :
Figure 7: Simplified Feynman rules for computing the squared matrix element directly.The first propagator function is to be used in all internal lines and the special DMP propagator (see text for discussion), shown by red color applies for the outgoing line in the diagram.For the definition of quantities in the vertex factors see figure 4.
factor γ k = m/γ k and the coefficients A a khl , B b khj , C c khi and D c khl are listed in table hγ • kγ 5 a/γ kl b/γ kj c/γ ki γ 5 ; hα • k ahv kl bhv kj −chv ki

Table 1 :
Listed are the coefficients in reducted tensor in equation (A.3) for the spinor structures.Finally, for the Dirac structures γ 0 γ 5 and hγ • k the projection tensor is given by