Flavour mixing transport theory and resonant leptogenesis

We derive non-equilibrium quantum transport equations for flavour-mixing fermions. We develop the formalism mostly in the context of resonant leptogenesis with two mixing Majorana fermions and one lepton flavour, but our master equations are valid more generally in homogeneous and isotropic systems. We give a hierarchy of quantum kinetic equations, valid at different approximations, that can accommodate helicity and arbitrary mass differences. In the mass-degenerate limit the equations take the familiar form of density matrix equations. We also derive the semiclassical Boltzmann limit of our equations, including the CP-violating source, whose regulator corresponds to the flavour coherence damping rate. Boltzmann equations are accurate and insensitive to the particular form of the regulator in the weakly resonant case $\Delta m \gg \Gamma$, but for $\Delta m \lesssim \Gamma$ they are qualitatively correct at best, and their accuracy crucially depends on the form of the CP-violating source.


Introduction
Quantum coherence, and the related mixing and oscillation of quantum states, plays an important role in many interesting phenomena in particle physics and in the early universe. Examples include neutrino mixing and oscillations [1][2][3][4][5][6][7][8][9], particle scattering from phase transition walls in baryogenesis [10][11][12][13][14][15][16][17][18], mixing of Majorana neutrinos in leptogenesis [19][20][21], and particle production after inflation or in phase transitions [22][23][24]. In these applications a quantum field theoretical treatment of coherence is essential. Several different types of coherence may be relevant depending on the problem: particle production is driven by the particle-antiparticle coherence, and coherent mixing of left-and right-moving particles can be the engine for creating the particle-antiparticle asymmetry during the electroweak phase transition. Finally, coherence between different flavour states powers the familiar phenomenon of neutrino oscillations as well as the particle-antiparticle asymmetry generation in leptogenesis, which is the main topic of this paper. The primary goal of this work is to develop a general and transparent formalism for treating problems involving quantum flavour coherence. While our results are more general, a large part of the development will be done in the context of (resonant) leptogenesis. We will develop a series of implementations of quantum transport equations with different levels of sophistication. At the most general level our formalism includes also the particle-antiparticle coherences. From there we move by a well defined reduction process to equations fully describing the flavour coherence separately in the particle and antiparticle sectors. These equations are valid for arbitrary helicities and neutrino masses. We also show how these equations can be further reduced to a helicity-symmetric density matrix equation and eventually to the Boltzmann limit with a soundly motivated form for the CP-violating source in the leptogenesis application. This hierarchy of implementations serves both to compare our results against the existing literature, and provides a "library of methods" from where one can choose the one that is best suitable for the given application at hand.
In the leptogenesis scenario [19] an initial lepton asymmetry n L gets converted to the baryon asymmetry via the B + L-violating sphaleron processes [25] present in the standard model (SM) [26][27][28][29]. Several versions of the leptogenesis mechanism exist [20,21]. In standard thermal leptogenesis n L is generated by CP-violating out-of-equilibrium decays of heavy Majorana neutrinos and the same basic mechanism carries over to resonant leptogenesis. In the latter the neutrino masses are nearly degenerate however, which leads to an enhanced efficiency. The enhancement is maximal when the timescale ∼ 1/∆m of the flavour oscillations is comparable to the timescale ∼ 1/Γ of the change of the neutrino abundances. In this case the flavour mixing has been found to be the dominant source of the lepton asymmetry [20,21]. Thermal leptogenesis can be well treated using standard kinetic equations, but the resonant case needs a more refined treatment to correctly account for the flavour mixing.
Our formalism is a generalisation of the coherent quasiparticle approximation (cQPA) first developed in [72][73][74][75][76][77][78] and further studied in [79]. The cQPA is a two-step approximation where the structure of the Wightman function is solved first in a collisionless approximation in the Wigner representation. This results in a spectral shell structure, including particle and coherence shells, which is then used to solve the full dynamical equation. The present method is not restricted to spectral structures as the equations are solved directly in the two-time representation, which allows taking into account a finite width of the pole propagators. If the width is neglected however, our equations are equivalent to cQPA equations.
Our main results include the quantum transport equations (5.7) and (5.15), the helicity symmetric equations (8.1)-(8.3) and the Boltzmann equation source term (8.10) with the CP-violating parameter (8.11). We provide an explicit implementation for a benchmark model with two Majorana neutrinos and one lepton flavour including decay and inverse decay interactions, but generalising to more complicated fermion sectors and scattering processes would be straightforward. Numerically all approaches are in good agreement in the weakly resonant case, m ∆m Γ, and the helicity-symmetric equation is in good agreement with the full master equation for all parameter values. However, when ∆m Γ the Boltzmann equation results depend strongly on the choice of the CP-violating source, the precise form of which has been debated in the literature [55,62,66,70]. Our result agrees with ref. [62]. We show that the effective width that defines this source corresponds to the off-diagonal damping of flavour coherence in the density matrix equations. We also find that the Boltzmann equations equipped with this source are in best numerical agreement with the full master equation results.
This paper is organised as follows: in section 2 we review the underlying CTP formalism and examine the general structure of solutions. In section 3 we derive the transport equation for the local correlation function S < k (t, t) and the key approximation leading to a closure is introduced in section 3.2. In section 4 we introduce the leptogenesis model, compute the self-energy functions for the decay processes, find the adiabatic background solutions and adapt the transport equations of section 3 to the leptogenesis case. We conclude the section with renormalised master equations for leptogenesis, including an equation for the lepton asymmetry and explicit forms for the source and washout terms. In section 5 we project the neutrino master equation onto different frequency, helicity and flavour quantum states, recasting it as a generalised density matrix equation, which we then further simplify by averaging over the particle-antiparticle oscillations. In section 6 we generalise to the case of an expanding universe and in section 7 we give detailed numerical results for the lepton asymmetry in some benchmark cases. In section 8 we introduce further approximations, first by dropping the helicity dependence and then reducing the master equation to a density matrix equation in the quasidegenerate case and finally into a Boltzmann equation in the decoupling limit, including a semiclassical source term. In section 9 we present detailed comparisons to earlier work in the literature. Further details of the derivation are presented in several appendices. Finally, section 10 contains our conclusions and outlook.

Kadanoff-Baym equations
For completeness and to introduce the notations, we start with a brief review of the CTP formalism. The basic quantity of interest in the CTP quantum transport theory is the contour-time ordered two-point correlation function. For flavoured fermions it is defined by where C is a complex time contour and the expectation value is defined as a trace weighted by the non-equilibrium density operator of the system. We will usually suppress the Dirac (α, β) and flavour (i, j) indices when possible, as they follow the spacetime coordinates u and v of the fermion field ψ. All products involving the two-point and self-energy functions are thus implicitly matrix products in Dirac as well as flavour indices. In this paper we consider  (2.2b) and the retarded and advanced pole propagators 1 Here A(u, v) ≡ 1 2 {ψ(u), ψ(v)} is the spectral function, which satisfies A = i 2 (S > + S < ) = i 2 (S r − S a ). We also denote A = S A below. We also need the Hermitian part of the pole propagators S H ≡ 1 2 (S r + S a ). The fermion self-energy Σ can be divided into real-time components similarly and the various self-energy functions Σ <,> , Σ r,a , Σ A,H satisfy analogous relations.
With the leptogenesis application in mind, we now specialize to spatially homogeneous and isotropic systems where the correlation and self-energy functions F ∈ {S, Σ} have spatial translational invariance: F (u, v) = F (u 0 , v 0 ; u − v). In this case it is convenient to use the two-time representation F k (t 1 , t 2 ) ≡ d 3 r e −ik·r F (t 1 , t 2 ; r), (2.4) which is just the spatial Fourier transform and k is the Euclidean three-momentum vector. The Schwinger-Dyson equation obeyed by the two-point function (2.1) can be cast into four real-time Kadanoff-Baym (KB) equations for the real-time correlation and self-energy functions. In the homogeneous and isotropic case and in the two-time representation (2.4) they are given by where p = r, a and s = <, >. The convolution appearing here is defined by where t in is the initial time of the CTP, and we have taken the limit t f → ∞ for the final time. We will eventually also take the limit t in → −∞ when deriving the kinetic transport equations. Finally, the inverse free propagator in equations (2.5) is given by The time dependent, real mass matrix m(t), can be viewed as a singular contribution to the Hermitian part of the self-energy, but with our immediate application to leptogenesis in mind, it is useful to write it explicitly. Later on we make frequent use of the Hermiticity relations of the propagators and the self-energies. In the two-time representation we can write them as where we defined S <,> ≡ iS <,> γ 0 , S r,a ≡ S r,a γ 0 and S A,H ≡ S A,H γ 0 . The relations (2.8) follow straightforwardly from the definitions (2.2)- (2.4). Note that Hermitian conjugation exchanges the pole propagators in addition to their time arguments. Similar Hermiticity properties hold among the self-energy functions (with S replaced by Σ in equations (2.8)) with the definitions Σ <,> ≡ γ 0 iΣ <,> , Σ r,a ≡ γ 0 Σ r,a and Σ A,H ≡ γ 0 Σ A,H . We will also need the spectral sum rule, which in the two-time representation reads 2A k (t, t) = 1 for any t.
(2.9) Equation (2.9) is a direct consequence of the canonical equal-time anti-commutation relations, and it can be also derived from the definitions (2.3), the analogous relations for the self-energy functions and equations (2.5a) and (2.7). The coupled integro-differential equations (2.5) would be difficult to solve even if the self-energies were some externally given functions. The fact that self-energies are in general functionals of the correlation functions S r,a,<,> , makes (2.5) also non-linear and consequently much more complicated. Our main goal is to find a tractable and efficient approximation scheme for these equations, which still captures the relevant physics.

Formal solutions
Before introducing our approximations, it is useful to study some general properties of the solutions. We start by rewriting (2.5) in an alternative, but equivalent form of Schwinger-Dyson integral equations: where the free propagators S p 0 and S s 0 satisfy and again p = r, a and s = <, >. 2 We suppressed the now obvious coordinate dependencies and used 1 to denote the functional identity matrix (i.e. the delta function in time variables). We can turn (2.10) into formal solutions by iterating and rearranging the infinite series solutions suitably (intermediate steps of the procedure are given in appendix A). The formal pole propagator solutions are given by The Wightman functions S s k can be broken into homogeneous and inhomogeneous parts: (2.14) The inhomogeneous part S s inh,k is a particular solution to the full KB equation (2.5b), while the homogeneous part S s hom,k satisfies the same equation with the right-hand side put to zero. Solutions (2.12)-(2.14) are still completely general, but indeed purely formal because the self-energies in general depend on S k . 3 However, they provide important insight as to how to define and set up approximative solutions and equations. For example, if the self-energy functions are dominated by some known part, such as the equilibrium contribution, they suggest how to split the pole-functions and the inhomogeneous part of the Wightman function to a leading part and a perturbation, following schematically the division Σ k = Σ eq,k + δΣ k . We shall put this observation into good use in section 3 below.

Homogeneous solutions
We can gain further insight to the homogeneous solution (2.13), simplifying it further (still in full generality) by using the KB equations. To this end we first substitute the inverse free propagator (2.7) to equations (2.5) and write them explicitly as Here H k (t) ≡ α · k + γ 0 m(t) is the free Dirac Hamiltonian with α ≡ γ 0 γ, and we used the barred propagators S and self-energies Σ defined below equations (2.8). Using the pole Figure 2. Left panel: schematics of the propagation of the homogeneous correlation from the initial condition set at (t 1 , t 2 ) = (t in , t in ) (red) and the extent of the non-local correlation around a local solution defined at (t 1 , t 2 ) = (t, t) (green, see section 3.2). The strength of the colour illustrates the strength of the correlation. Outside the spheres, with |∆t| 1/γ k , points are no longer correlated (with (t in , t in ) and (t, t) respectively) because of the dissipation. Right panel: shown are the particular time-contours of relevance for the local equation as discussed in section 3.3.

equation (2.15a) and its Hermitian conjugate in equation (2.13) then leads to
To get the second equality we integrated the derivative terms by parts, used equation (2.11b) and its Hermitian conjugate and the definitions (2.3). The only remaining term in the expression is then the boundary term involving the initial time t in . Note that it is essential in this derivation that t in is smaller than both t 1 and t 2 . Finally, we used the sum rule (2.9) at t = t in to identify S s 0,k (t in , t in ) = S s hom,k (t in , t in ). The result (2.16) is exact and it shows that homogeneous solutions describe transients evolving from initial conditions set at a finite initial time t in . Indeed, it is easy to show that the spectral function is the unitary time-evolution operator for two-point functions in the free theory limit: 2A 0 (t 1 , t 2 ) = U (t 1 , t 2 ), for which the spectral sum rule (2.9) imposes the correct normalisation U (t, t) = 1. However, in dissipative systems this time-evolution is no longer unitary. Dissipation shows up as a finite width of the pole functions in the Wigner representation. Going to the Wigner representation, finding quasiparticle poles and computing the corresponding quasiparticle widths γ k , and transforming back to the two-time representation, one could then obtain a more general A k effecting an explicit non-unitary decay of correlations in the relative time [79]: As a result, if the initial time is pushed to the past infinity, t in → −∞, all memory of the initial conditions will vanish, leaving only the inhomogeneous solution. We illustrate this phenomenon schematically in figure 2.
We find the homogeneous solutions of the form (2.16) to be very useful still in another way, which is crucial to our scheme. We will return to this issue in section 3.2 below.

Local quantum kinetic equation
The main source of complexity in equations (2.15) comes from their non-locality and the associated need for a complete accounting of the memory effects. Yet, all physical observables are expressible in terms of the local correlation function S k (t, t), and moreover, we have seen that dissipative processes in general wash out memory effects over time intervals |∆t| ≥ 1/γ k . This suggests to try to find an approximative equation that involves only the local correlation function. Such an equation is easy to set up formally: we first use the chain rule to write 5 Using this with the KB equation (2.15b) and its Hermitian conjugate together with the Hermiticity relations (2.8) yields the equal-time equation This equation is still exact, but of course not closed because it still involves the non-local function S < k (t, t ) with t = t explicitly in the interaction convolution terms Σ r k * S < k and S < k * Σ a k and implicitly within the self-energies. It also depends explicitly on the pole functions S r,a k (t 1 , t 2 ). To make (3.2) self-contained, we need to supply it with enough information of these particular correlations, without going back to the full KB equations.

Perturbation around the adiabatic solution
Let us first address the coupling between the Wightman functions and the pole functions. This issue is intimately connected to finding a good approximation for the inhomogeneous solution for the Wightman functions as is suggested by equations (2.12) and (2.14). We start by formally dividing the correlation functions into some known adiabatic background solution and a perturbation: where α = r, a, <, >. There is some freedom as to how to choose the adiabatic solutions. For example, in a system near thermal equilibrium, the instantaneous thermal equilibrium solutions would be an obvious choice. More formally, the adiabatic Wightman function can be defined as a solution that reduces to the stationary solution of (3.2), when ignoring all 5 In order to get an equation for the local correlation function, it is essential that the derivative acts on both time-variables. Otherwise the limiting procedure introduces other independent functions (first moment of the propagator in the mixed representation) to the left hand side of the equation (3.2).
local time dependence. 6 When S α ad,k fulfil this requirement, inserting the division (3.3) into equation (3.2) leads to The source −i∂ t S < ad,k (t, t) is the leading correction left from the terms involving the adiabatic solution in the full dynamical equation. Note that as the Wightman function can contain also a homogeneous transient, δS < k is not necessarily small even when the adiabatic solution is an equilibrium solution.
The equal-time pole functions, on the other hand, can be taken to be purely adiabatic with no dynamical perturbations: δS p k = 0. In the equal-time limit this condition is actually strictly imposed by the spectral sum rule (2.9) and the relations (2.3) between the pole functions and the spectral function. One can also show that ∂ t S p ad,k (t, t) = 0 exactly at least in a non-interacting theory. This suggests that no perturbations in the pole functions can be included in a truncation to the local limit. It is then remarkable, that in the truncation scheme for the collision terms developed in the next section, the non-local pole function perturbations vanish consistently with their vanishing equal-time counterparts. This means that pole functions become entirely non-dynamical quantities that account for the structure of the phase space only. With this information, equation (3.4) further reduces to This equation no longer depends explicitly on the pole functions. They only have limited influence on (3.5) through the source function and the self-energy functions, as we shall discuss below.
On the choice of the adiabatic solution. Explicit forms for the adiabatic solutions are most conveniently given in the Wigner representation, which is defined as the Fourier transform of the two-time representation (2.4) with respect to the relative time-coordinate: Here t = (t 1 + t 2 )/2 is the average time coordinate and k 0 is the internal energy conjugate to the relative time coordinate r 0 = t 1 − t 2 . Note that we always take the limit t in → −∞ before calculating any Wigner transforms. We now define the adiabatic solutions with instantaneous mass m(t) and self-energies as follows 7 : There still is significant freedom left in these solutions related to the choice of the self-energy functions Σ α ad . For example, one might choose to ignore or include the Hermitian part Σ H ad of the pole self-energy functions, leading to solutions either with vacuum or quasiparticle dispersion relations. Moreover, if one neglects the finite width, setting Σ A ad = 0, the corresponding solutions become spectral (either vacuum or quasiparticle). This is the choice made in the derivation of cQPA-formalism [72][73][74][75][76][77][78][79], as well as in the usual Boltzmann theory [82]. If one includes a finite width however, the adiabatic part of the solution spreads out in phase space, with a consequent change in the source term in the equation (3.5) for the perturbation δS < k . One can even include corrections from the perturbations δS < k in the self-energies, without changing the basic structure of the equation for the perturbation δS < k itself. The point is that the validity of all these approximations is controlled by the coupling constant expansion.

Local approximation
The division into an adiabatic background and the perturbation simplified the original equaltime equation considerably. However, the problem of closure still remains: our equation describes the evolution of the perturbation only along the diagonal t 1 = t 2 in figure 2, but the collision integrals depend on δS < k (t 1 , t 2 ) everywhere in the two-time plane. In a system with dissipation, these memory effects are suppressed however, and there is hope that a strictly local description can be found. We use the evolution of the homogeneous perturbation (2.16) as our guiding principle to reach this goal.
Indeed, imagine first that we somehow have found the correct solution δS s k (t, t) along the diagonal. For t 1,2 not too far from the time t, the true non-local solution δS s k (t 1 , t 2 ) should be correlated with δS s k (t, t), and even reasonably well approximated by a homogeneous solution similar to equation (2.16), with the initial time t in replaced by the local time t. Even when we do not know the local solution beforehand, we can parametrise the non-local solution with the local one. Specifically, we make the local ansatz for any t. (3.8) Note in particular that the spectral function is not a dynamical quantity in (3.8); following the discussion of the previous section, the generalised time evolution operator 2A k (t 1 , t 2 ) is a non-dynamical adiabatic solution that can be computed to the desired accuracy independently from the local correlation function δS s k (t, t). We stress that the ansatz (3.8) will only be used in the convolution terms describing the interactions. Indeed, if taken to hold universally for all t, it would generally be too restrictive and in contradiction with the local equation of motion (3.5).
cQPA and Boltzmann theory limit. It should be stressed that the ansatz (3.8) is an exact relation in free theory and for spectral quasiparticles, where the full solution is homogeneous and the free spectral function is the unitary time evolution operator. Indeed, we can derive the cQPA-correlation function [77,79], and eventually the Boltzmann theory limit directly from (3.8). First choosing t = (t 1 + t 2 )/2 and then Wigner-transforming (3.8), one finds (here we write the results explicitly for S < but analogous results hold for S > ) . This form is still valid for any adiabatic solution for the spectral function. Working to lowest order in gradients and using the free adiabatic spectral function, equation (3.9) reduces to the spectral form where the helicity and energy projection matrices are defined by withĥ k ≡ α ·k γ 5 , H ki = α · k + γ 0 m i and ω ki = |k| 2 + m i 2 . The helicity and energy indices both take values h, s = ±. It is easy to show that the projectors P kh P ± ki γ 0 P ± kj form a complete basis of matrices consistent with homogeneity and isotropy (this will be elaborated further in section 5.1). In the spectral limit the adiabatic solutions and the perturbations can be combined on the common shell functions. Expanding the corresponding full S < k (t, t) in this basis (for the precise definition of P ss khij see equation (5.2) below), we can rewrite equation (3.11) as where ω kij ≡ (ω ki + ω kj )/2 and ∆ω kij ≡ ω ki − ω kj . This is the flavoured cQPA-propagator, up to normalisation, derived in [77] and it carries information of all coherence structures consistent with homogeneity and isotropy in the spectral limit.
If one ignores all coherence information, equation (3.13) reduces to a generalised KBansatz which corresponds to the Boltzmann theory limit with distribution functions that are diagonal in flavour and helicity. The extra sgn(k 0 ) factor could be absorbed to normalisation, but the present normalisation will be more convenient later. Moreover, if one imposes the thermal equilibrium Kubo-Martin-Schwinger (KMS) condition for the cQPA Wightman functions, S > (k, t) = e βk 0 S < (k, t), which now is equivalent to S < ij (k, t) = 2A , the distribution function further reduces to the thermal Fermi-Dirac distribution: f m,± khii = ±f FD (±ω ki ) (cf. [79]).

Local transport equation
For us, the most important utility of the local approximation (3.8) is that it allows closure in equation (3.5), reducing all interaction convolutions containing the non-local function δS < k (t, t ) to simple matrix products involving only the local function δS < k (t, t). For example, where we introduced the effective self-energy Σ r eff,k . We remind that the spectral function is adiabatic A ≡ A ad in these expressions, consistent with our approximation scheme. While the effective self-energy is still a convolution, it can be computed at any time during the solution based only on the local solution itself, or independently of it, depending on the approximation one uses for the adiabatic functions, as discussed in section 3.1.
We now use the local approximation (3.8) to obtain closure in equation (3.5). This amounts to using (3.15) and its Hermitian conjugate with (3.16) in equation (3.5), resulting in the local equation of motion Equation (3.17) is our final quantum kinetic equation (QKE) for non-equilibrium evolution of mixing fermions. The non-local memory integrals have been truncated by the local approximation, so it is an ordinary (matrix) differential equation for the local non-equilibrium correlation function δS < k (t, t). Equation (3.17) still describes both flavour and particleantiparticle coherence effects of the mixing fermions. It also takes into account quantum statistical effects of the thermal medium (within the weak coupling expansion), and it can accommodate thermal corrections to the dispersion relations via the effective self-energy and the adiabatic source term. We have shown that (3.17) encompasses the coherent cQPA-formalism and consequently the usual Boltzmann theory including also semiclassical corrections [79], but it is a more general formulation in that it is not restricted to the spectral limit. We will apply this equation in the leptogenesis setting to describe the evolution of the right-handed Majorana neutrinos in the next section.
On the accuracy of the local ansatz. Despite its wide range of applicability, the ansatz (3.8) should eventually break down if the system develops significant temporal correlations (a memory) over large time intervals. When would this happen and how large would the corrections be? Ultimately one would like to compare the results obtained using the local equation (3.17) with a numerical solution of the full non-local two-time equations (2.15), but we can get a good idea of the size of the memory effects by studying their origin in the Wigner representation.
We first note that the Wigner transform (3.6) encodes all dependence on the relative time r 0 = t 1 − t 2 at constant average time slices t = 1 2 (t 1 + t 2 ) in the frequency components S k (k 0 , t) (see the right panel of figure 2 for illustration). In the weak coupling limit this information gets concentrated on narrow shells in frequency space and eventually to spectral solutions when widths are neglected. The non-local information relevant for equations (3.5) and (3.17) is contained in the convolution integrals that can be expressed as follows: where G = δS < in (3.5) and G = A in (3.17). 8 Clearly, for a fixed t the non-local information of G k (t , t) contained in the two-time convolution (along the contour t 2 = t) is fully encoded in the gradients in the Wigner representation, since G(k, t) only contains information along the contour 1 2 (t 1 + t 2 ) = t. This correspondence is schematically illustrated by the blue arrows in figure 2.
Finally then, the validity of the local approximation (3.8) boils down to the smallness of the gradient corrections and assuming that δS(k, t) has a similar phase space structure as the adiabatic solution. In the leptogenesis application gradient corrections are controlled by the Hubble expansion and hence they are small since H/T 1. The phase space structures of the adiabatic solution and the perturbation δS(k, t) should be similar because the latter is created by the former. Also, both solutions become spectral when the width is zero, so this approximation becomes good also in the weak coupling limit.

Leptogenesis
We now apply our methods to study lepton asymmetry production in the early universe. Leptogenesis has different variants, including the original thermal leptogenesis [19], resonant leptogenesis [83,84] and the freeze-in, or Akhmedov-Rubakov-Smirnov (ARS) leptogenesis [85]. For more discussion see e.g. [20,21]. Our methods apply, with minor modifications, to all these variants, but we will focus to the resonant scenario in the minimal model with two heavy Majorana neutrinos coupled to a single light SM lepton doublet and a Higgs doublet. Generalisation to more neutrino flavours, or more SM lepton flavours necessary e.g. for low scale leptogenesis [86,87] and the ARS-mechanism, would be straightforward. We will also only include the decay and inverse decay interactions, neglecting the 2 → 2 scattering processes. This limits the range of validity of our predictions, but our goal is not the maximal phenomenological reach, but the accuracy of the quantum transport formulation and detailed comparisons between different approximations. Again, generalisation to more complex interactions would be straightforward. 8 There is one subtlety if we take G = δS < as in equation (3.5). In this case one has to account for the gradient operator in the argument of Σ r out , when it acts on the rapidly oscillating coherence solutions in δS < . This ensures that the coherence shell contributions get computed on correct frequency shells in the cQPA-formulation (see e.g. [79]). One of the nice features of the ansatz (3.8) is that it fully automatises this resummation, also when evaluating higher loop self-energy functions [77]. Indeed, the issue clearly does not arise when G = A, because A is an adiabatic function.
In principle all particle species involved in the problem could be treated on the same footing in the CTP-context, resulting in a network of local transport equations. However, a number of simplifying approximations can be made for the SM fields. For example we can neglect the decay widths of the lepton and Higgs fields and assume that they are in kinetic equilibrium due to the SM gauge interactions. To first order, we can also neglect the chemical potential of the Higgs field, which then decouples from the dynamics. The lepton chemical potential µ is essential of course, but we can assume µ /T to be small, which allows us to neglect the backreaction of µ to the dynamics of the Majorana neutrinos. We can then solve the coupled neutrino and lepton equations consecutively instead of simultaneously. For neutrinos we will use the local transport equation (3.17) with full phase space structure and flavour coherence information. Since we are only accounting for one-loop self-energies, we include only the indirect or ε-type self-energy contribution to the CP-violation. Including the direct ε -type contribution would require a two-loop self-energy calculation, but we refrain from doing this, because the indirect contribution is dominant in the resonant regime.
In summary, the objective of this section is to derive an explicit local quantum transport equation for the mixing Majorana neutrinos, coupled with an equation for the particleantiparticle asymmetry of the SM lepton doublet. We will compute explicitly all self-energy functions and the adiabatic background solutions for the neutrinos, as well as the source and washout terms for the lepton asymmetry equation. We will work consistently to the leading order in the coupling constant expansion and discuss the renormalisation procedure necessary for the loop calculations.

The minimal model
Our model contains two singlet Majorana neutrino fields N i , coupled to an active lepton SU(2)-doublet and the complex scalar Higgs doublet φ via chiral Yukawa interactions: We work in the mass basis for the Majorana neutrinos, where m i are the lepton number violating real Majorana masses. The lepton and Higgs fields are assumed to be massless as leptogenesis must take place in a high temperature in the electroweak symmetric phase. Finally, the SU(2)-conjugate Higgs doublet is defined by φ ≡ εφ * where ε is the antisymmetric 2 × 2 matrix with ε 12 = 1. The CP-violating phases necessary for leptogenesis are contained in the complex coupling constants y i . The CTP propagators of the neutrino, lepton and Higgs fields are given by where C is the unitary and antisymmetric charge conjugation matrix. An important consequence of the condition (4.3) is that the neutrino propagator is constrained by Note that here the transpose acts on all of the flavour, Dirac and CTP indices of the propagator.

Self-energy functions
We calculate the self-energies in the 2PIEA formalism, where the self-energies are defined as functional derivatives of the non-trivial part Γ 2 of the effective action with respect to the propagators. The non-trivial part Γ 2 contains contributions of vacuum diagrams with two or more loops. The two-loop contribution (see figure 3a) arising from the Yukawa interaction in the Lagrangian (4.1) is given by Here c w = 2 is a multiplicity factor from the trace over the SU(2)-doublet indices, and C denotes integration over the CTP. To calculate the Majorana neutrino self-energy we need to take the functional derivative of Γ 2 with respect to the neutrino propagator S. As the propagator is constrained by the Majorana condition (4.4), its functional derivative is defined by where the flavour and Dirac indices have been indicated explicitly. For an unconstrained Dirac particle the second term on the right-hand side of equation (4.6) would not be present. The one-loop Majorana neutrino self-energy Σ ij (figure 3b), calculated using equations (4.5) and (4.6), is then given in direct space by and (2.3): for example the Wightman functions are The only difference to the fermionic case is that the bosonic spectral function is a commutator of the fields, rather than an anti-commutator, and we use the standard sign convention ∆ < = ∆ +− for bosons. The lepton self-energy (figure 3c) is calculated similarly, and the result for the SU (2) From these it is straightforward to calculate the real-time self-energies by inserting the CTP indices (which follow the spacetime arguments) and using the relations of the different propagators. Later we need the following neutrino self-energy functions (given now in the Wigner representation): where we defined the statistical propagators S F ≡ 1 2 (S > − S < ) and ∆ F ≡ 1 2 (∆ > + ∆ < ).

Tree level propagators
In order to calculate the self-energies further we now introduce the tree level equilibrium approximations for the lepton and Higgs propagators S and ∆. We assume that the diagonal elements of the SU(2)-symmetric lepton and Higgs correlators are given, in the Wigner representation, by where PV denotes the Cauchy principal value and µ is the lepton chemical potential. The ± sign in equations (4.11c) corresponds to the ≶ sign (+ for < and − for >). The Fermi-Dirac and Bose-Einstein phase space distribution functions are where β = 1/T and we assume a common temperature T for both and φ. Note that the time-dependence of the lepton Wightman functions S <,> (k, t) in (4.11c) comes solely from the chemical potential µ = µ (t).
Similarly to the Majorana neutrino correlation function, it is convenient to split the lepton Wightman functions into the equilibrium and non-equilibrium parts, where S <,> ,eq (k) is given by equation (4.11c) with µ = 0 and the non-equilibrium parts satisfy δS > = −δS < . It then suffices to consider e.g. δS < only, for which we define the non-equilibrium lepton distribution Here we also assumed that the lepton chemical potential µ remains small during leptogenesis.
We can now split the Majorana neutrino self-energies (4.9) and (4.10) using equation (4.13): where Σ eq is the thermal equilibrium part with vanishing µ , and δΣ is proportional to δS and hence linear in µ in the approximation (4.14).
Using equations (4.9)-(4.12), (4.14) and (4.15), together with Σ A ≡ i 2 (Σ > + Σ < ), we can calculate all the needed Majorana neutrino self-energies. We parametrise them as where the various S µ functions are calculated explicitly in appendix B. Due to homogeneity and isotropy there are only two independent functions S 0 (k) ≡ a(k 0 , |k|) and does not contribute to our results in the leading order approximation considered in this work.

Adiabatic neutrino solutions
Now that we have specified the self-energy, we can calculate the adiabatic solutions (3.7) for the neutrino which are needed in the kinetic equation (3.17). Working to the lowest order approximation we only take the equilibrium part of the neutrino self-energy (4.15) into account: We are implicitly assuming that all quantities have been renormalised, so the pole selfenergy Σ p eq appearing here is actually the on-shell renormalised pole self-energy Σ p eq given by equations (4.44) and (4.51). Furthermore, using the KMS relation Σ > eq (k) = e βk 0 Σ < eq (k) and the exact identity A = S r Σ A S a , which holds in direct space as a convolution and as a simple product to zeroth order in gradients in the Wigner representation, we can write the Wightman functions (4.17b) as In perturbative expansions it is more convenient to use the form (4.18) than equation (4.17b). For example, a naive coupling constant expansion of the pole propagators in equation (4.17b) would appear to give a O(y 2 i ) result, whereas the right-hand side of equation (4.18) obviously gives the correct O(1) result. In other words, a consistent expansion of the Wightman function (4.17b) requires resumming the pole propagators, which the identity A = S r Σ A S a implements automatically in equations (4.18).
Given the adiabatic solution, we can now compute the source term −∂ t S < ad,k (t, t) in equation (3.17). To this end, we employ the inverse Wigner transform to calculate the two-time representation function S < ad,k (t 1 , t 2 ) from the Wigner representation (4.18a). Utilising Cauchy's integral theorem to perform the k 0 -integral, we can write the two-time representation as where the k 0 -sum is taken over the poles of S p ad (k, t), and s p = −1, 1 for p = r, a, respectively. We also used the short-hands r 0 ≡ t 1 − t 2 and t ≡ (t 1 + t 2 )/2. The corresponding result for S > ad,k follows by replacing . Note that we neglected any branch cuts in equation (4.20) and kept only the residue contributions of the single particle poles, but otherwise the result is general. Also, the poles of f FD (k 0 ) at the imaginary axis do not contribute due to the KMS relation.
Leading order approximation. Equation (4.20) is written with a general adiabatic pole function, but it will be eventually sufficient to compute it to lowest order in gradients and in the coupling constants. In this case we can use the free tree-level solution with vacuum dispersion relations for the adiabatic pole propagators S p ad (k, t) in equation (4.17a), given by The small imaginary part with ε > 0 ensures the correct boundary conditions for p = r, a.
Inserting this solution into equation (4.20), calculating the residues and taking the limit ε → 0 + yields the leading order result Here / k si ≡ sω ki γ 0 − γ · k and ω i ≡ |k| 2 + m i 2 , and the time-dependent masses m i (t) are evaluated at the average time t = (t 1 + t 2 )/2. 10 Note that the lowest order adiabatic pole and Wightman functions are diagonal in the mass basis in agreement with [63]. We emphasize that in our approach the off-diagonal corrections to the adiabatic propagators are not required to leading order, as we do not need to calculate the homogeneous solution non-dynamically from equation (2.16), as was done in e.g. [62]. Instead, corrections from the full adiabatic pole propagators are already taken into account in the source term of our dynamical equation (3.17) (see also equation (3.4)). More discussion about these corrections is given at the end of section 4.5.
Setting t 1 = t 2 ≡ t in equation (4.22) we then get the equal-time adiabatic function S < ad,k (t, t) needed in the kinetic equation (3.17). The result may also be cast into a compact form where the energy projection matrix P s ki was defined in equation (3.12). To calculate the time-derivative in order to get the source term −∂ t S < ad,kij (t, t) for equation (3.17) is now a simple matter, once the time-dependent masses m i (t) are specified. Equation (4.23) is also needed to calculate the initial values for the (local) non-equilibrium Wightman function δS < ≡ S < − S < ad once the initial value for the full function S < has been specified.

Effective neutrino self-energy
We now calculate the effective self-energy (3.16), which appears in the local kinetic equation (3.17) for Majorana neutrinos. We use the general result for a convolution in terms of the Wigner representation correlation functions: where F out (k, t) and G in (k, t) were defined below equation (3.9). The first line of equation (4.24) is the inverse Wigner transform of the Moyal product (often denoted by ♦). When the gradients are small we may further approximate F out (k, t) F (k, t) and G in (k, t) G(k, t). We will again assume an adiabatic, equilibrium self-energy function and expand to zeroth order in gradients. The general result is then similar to equation (4.20): On the last line we again kept only the residue contribution from the poles of the propagator and neglected contributions from possible branch cuts. To calculate the leading tree level approximation we will again use the free pole propagators like in equations (4.22) and (4.23). The calculation proceeds as before and one finds the lowest order coupling constant result Corrections resulting from resummed pole propagators could be included by using the residue formula (4.26) with complex poles, but as will be argued below, their effect would be formally of higher order (∼ O(y 4 i )) in the coupling constants.
Beyond the leading approximation. To improve on the leading approximation in coupling constants used above, one should solve the complex k 0 -poles from the determinant of S p ad (k, t), keeping the self-energy corrections and use them to calculate the residues in equations (4.20) and (4.26). This is in principle straightforward, but laborious because the block-wise inversion of equation (4.17a) results in complicated formulae due to the flavour mixing. We will not implement these corrections in this paper, because we use the weak coupling approximation where these corrections should be small. We give details of the inversion procedure in appendix C for completeness however, and note that a similar analysis was presented in [62], with a non-relativistic approximation for the self-energy. It should also be noted that these higher order corrections to (4.23) and (4.27) remain parametrically small also in the maximally resonant regime where ∆m 21 /m ∼ |y| 2 . We have verified this explicitly in a generic power counting expansion of the formula (4.20) where ∆m 21 /m and |y| 2 are expanded simultaneously, under the assumption that the difference of the full k 0 -poles to the corresponding free theory poles is O(y 2 ). To get the correct perturbative result it is also crucial to start from the resummed form of equation (4.18a).
The self-energy corrections to the pole propagators can be divided into dispersive corrections due to the Hermitian self-energy and dissipative ones due to the antihermitian self-energy part. The main effect of the latter was already discussed qualitatively in section 2.2. While the actual formulae are very complex, it is easy to see that the main qualitative effect of evaluating (4.20) at a complex pole is to introduce complex damping terms ∼ e −γ ii |t 1 −t 2 | into equation (4.22). Such factors represent the dissipation in the relative time coordinate, but eventually do not affect the source term in the local equation (3.17). Indeed, the only effect on the source, and likewise on the self-energy convolutions (4.26) appearing in (3.17), from a finite width then amounts to shifts ∼ y 2 i to the energy shells where these terms are evaluated. In the weak coupling limit such corrections should be small, which we have verified numerically.
The dispersive corrections could be more interesting. Including Hermitian self-energy corrections would lead to new real-frequency poles for the adiabatic functions. Combining the vacuum Hamiltonian in equation (3.17) with the effective Hermitian self-energy evaluated at these poles, would give rise to an effective matter Hamiltonian for the quasistates. The effective Hamiltonian would in general be no longer diagonal in the mass basis and the energy difference between the matter eigenstates would be a function of time, similar to the case of mixing light neutrinos in the early universe [3,4,7]. Such a dynamical modification of the energy level splitting could be relevant for resonant leptogenesis, although the analysis of ref. [65] (performed in a simplified setup) does not suggest that the effect is quantitatively significant.

Lepton transport equation
The equation for the lepton asymmetry can be derived from the KB equations (2.5) for the lepton propagator. However, as the lepton is massless and we use the tree-level spectral approximation for its propagator, the lowest order adiabatic source term in equation (3.4) vanishes. To derive the leading source for the lepton asymmetry, it is convenient to use a different but equivalent formulation of the KB equations (see e.g. equations 17 and 18 in [73]). The appropriate form of the equation where the finite width and dispersive contributions have been neglected is The corresponding local equation for the local correlation function of the lepton can be derived in the same way that we derived equation (3.2). The result is where α · k is the free Hamiltonian for the massless lepton. We also remind here that equations (4.28) and (4.29) are formulated for the diagonal element S of the SU(2)-symmetric lepton correlator (4.2b). This equation could be solved as such, coupled with the local equation for the mixing Majorana states. However, in practice we can make several further approximations, eventually converging to a simple scalar equation for the lepton asymmetry.
Lepton asymmetry. The lepton asymmetry we are interested in this work can be related to the chiral current density of the left-handed leptons, which is defined by where an implicit summation over the SU(2)-index is assumed. Since we consider a spatially homogeneous and isotropic system, the current depends only on the time x 0 = t and it can be further related to the local two-time Wightman function. Using definitions (2.2), (2.4) and (4.30), we then get Because no asymmetry can exist in thermal equilibrium [88], we can define the lepton asymmetry density n − n as the zeroth component of the non-equilibrium part of the current: We can further relate the asymmetry to the chemical potential of the tree-level lepton propagator. Using equations (4.11c), (4.12), (4.13) and (4.19) we calculate the trace on the right-hand side of equation (4.32). The result, written in terms of lepton and anti-lepton phase space distributions, is where c w = 2 is the weak isospin multiplicity factor. There is no additional spin multiplicity factor because of the chiral projection, i.e. only one spin state couples to the Majorana neutrino and develops an asymmetry in the massless limit. A standard calculation of the integral gives the relation between the asymmetry and the chemical potential: where in the last step we used the linear approximation for µ like in equation (4.14).
Lepton asymmetry equation. We can get the equation of motion for the lepton asymmetry (4.32) by using the split (4.13) in equation (4.29), taking the trace over spinor indices and integrating over momentum. The trace of the commutator term vanishes due to the Dirac structure of the tree level propagator, so that we are left with Substituting the two-time representation of the lepton self-energy (4.8) to the right-hand side then yields This equation can actually be expressed using the already calculated Majorana neutrino self-energy (4.9). Indeed, because of the trace and the equal time arguments it is just a matter of combining the terms differently in the right-hand side of equation (4.36) to rewrite the integral in terms of the Majorana neutrino correlation function and the chirally projected Majorana self-energy Σ L , which results in: One should note that the trace is now taken over both the Majorana neutrino flavours and the Dirac indices and we defined the barred chiral Majorana neutrino self-energy as Also note that the lepton doublet multiplicity c w is now included in the neutrino self-energy.
Results similar to (4.37) are already known in the literature [62], but the novelty of our approach is in the use of the local approximation of section 3 to evaluate the convolution integrals. This is now straightforward because equation (4.37) is written explicitly in terms of the Majorana neutrino propagator. We first use equations (3.3) and (4.15) to split the neutrino Wightman functions and self-energies into the equilibrium and non-equilibrium parts. Then we invoke the local approximation (3.8) to compute convolution integrals with the perturbations δS <,> , along with the constraint δS > k (t, t) = −δS < k (t, t), which is imposed by the sum rule. Finally, we write the convolution integrals involving S <,> ad in the Wigner representation using equation (4.24) and expand consistently to the leading order in coupling constants and gradients to obtain: The first term on the right-hand side of equation (4.38), proportional to Σ A eq , does not contain the lepton chemical potential µ . It is therefore the source term for the lepton asymmetry. The remaining terms, proportional to δΣ α , are linear in µ ∝ n − n (in the approximation (4.14)) and so they contribute to the washout.

General leptogenesis equations
To summarise, we use the local equation (3.17) with equilibrium self-energies to solve the evolution of the Majorana neutrinos and equation (4.38) to subsequently calculate the lepton asymmetry. Our coupled equations for leptogenesis therefore read where the CP-violating source term S CP and the washout terms δW and W ad of the lepton equation are given by The washout terms are proportional to the lepton asymmetry n − n via equation (4.34). The adiabatic source term −∂ t S < ad,k (t, t) of the neutrino equation (4.39) is calculated from equation (4.23), and the effective self-energy Σ r eq,eff,k (t, t) is given by equation (4.27). We expect W ad to be the dominant washout term because it is proportional to the adiabatic functions, as opposed to δW which depends on the non-equilibrium perturbation δS < only.
Equations (4.39)-(4.43) are fully general apart from our using the local ansatz (3.8) to compute the collision terms and the simplifications made in the reduction of the SM sector. The effective self-energies in the lepton source and washout terms (4.41)-(4.43) are all calculated expanding consistently to the leading order in gradients and in the coupling constant expansion (more precisely: they are first order in |y i | 2 , ∂ t m and ∂ t µ combined). This is the most compact form of the equations relevant for the leptogenesis problem. They correspond to an initial value problem for a set of coupled first order equations which is straightforward to solve numerically by discretising the momentum variable. We shall recast these equations into a set of coupled Boltzmann-like equations for the generalised phase space functions in section 5, after a short discussion of the issue of renormalisation.

Vacuum on-shell renormalisation
So far we have implicitly assumed that we are working with finite, renormalised quantities. The renormalisation procedure is slightly intricate due to the flavour mixing. For completeness, we perform the one-loop vacuum renormalisation in our model, following the on-shell method of ref. [89]. This is sufficient for our purposes since we do not consider gauge interactions [90,91]. In this section we denote renormalised quantities by a hat. The renormalised pole self-energies Σ p eq,ij (k, t) can be written in terms of the unrenormalised functions and the counterterms as follows: The complex conjugation is here understood element-wise and not in the matrix sense. The mass counterterms δm i are diagonal in the vacuum mass basis, but the wave function renormalisation factors δZ L,R ij are in general flavour matrices [89]. Because our neutrinos are Majorana fields, the counterterms for different chiralities are related by Because of the Hermiticity of the counterterm Lagrangian, only the dispersive part of the self-energy Σ p eq = Σ H eq + is p Σ A eq contributes to renormalisation. The absorptive parts are finite as such and can be understood as being computed in terms of renormalised parameters throughout. Also thermal corrections are purely finite and may be split from the vacuum parts according to equations (B.1a) and (B.2) given in appendix B. Renormalisation is then associated only with the vacuum part of the Hermitian self-energy function Σ H(vac) eq . The on-shell renormalisation conditions which guarantee that there is no mixing in the external legs, that m i are the renormalised masses and that the residue of the diagonal propagator is unity, are given by [89,92,93] Note that there is no summation over repeated indices here. The dimensionally regularised vacuum self-energy is given by Σ is the spacetime dimension, γ E is the Euler-Mascheroni constant and µ is the renormalisation scale parameter. Using these results we find the following counterterms for the Majorana neutrinos: where i = j in the first equation. The corresponding renormalised vacuum self-energy is The full renormalised pole self-energy (4.44) can now be written as Note that the renormalisation procedure associated with the processes relevant for leptogenesis is not affected by the expansion of the universe. Indeed, as long as curvature effects are not relevant for the physical processes involved, renormalisation can be carried out in the local orthonormal coordinate system, which is locally a Minkowski space. We shall introduce the extension of our equations to the expanding Friedman-Robertson-Walker spacetime in section 6.

Non-equilibrium distribution functions
The local correlation function δS < k (t, t) is a matrix in both Dirac and flavour indices and its components have a complicated time dependence involving many different scales. These scales reflect the complicated phase space structure of the underlying Wigner function δS < k (k 0 , t), and they ultimately arise from the different physical phenomena the correlation function describes. In particular, the vast difference between the particle-antiparticle oscillation time ∼ 1/ω k and the flavour oscillation time ∼ 1/∆m makes equation (4.39) challenging for studying resonant leptogenesis as such. To overcome this problem we will parametrise δS < k (t, t) in terms of phase space distribution functions δf ss khij , and derive their coupled equations of motion. The benefit of this parametrisation, first introduced in [77], is that each phase space function δf ss khij describes separate, clearly defined physics with characteristic time-dependence. This allows us to isolate the physics that we are interested in and to write simplified and yet accurate versions of equation (4.39), that are amenable to efficient numerical solution.

Projection matrix parametrisation
Since we consider a spatially homogeneous and isotropic system, we can construct δS < k (t, t) using only 8 of the 16 basis elements of the full Dirac algebra. The basis matrices of this subalgebra commute with the momentum representation of the Dirac helicity operator, wherek ≡ k/|k| is the unit three-momentum vector. As mentioned above, we will use the parametrisation introduced in [77], and which we already used in the spectral solution (3.13): Here P kh and P s ki are the helicity and energy projection matrices defined in equation (3.12) and are labelled by helicity h = ±1, neutrino flavour i, j and the energy sign indices s, s = ±1. These matrices obviously satisfy the completeness relations P k+ + P k− = P + ki + P − ki = 1 and it is easy to show that they also satisfy the idempotence and orthogonality relations P kh P kh = δ hh P kh and P s ki P s ki = δ ss P s ki . It can also be shown that for any given k, i, j, the four different matrices P s ki γ 0 P s kj (with s, s = ±1) span the same set as {1, γ 0 , γ ·k, α ·k} which commute with the helicity operatorĥ k . Thus the eight matrices P kh P s ki γ 0 P s kj (with h, s, s = ±1) in equation (5.2) can be used as a basis for the entire homogeneous and isotropic Dirac subalgebra.
The normalisation factors N ss khij , which in part define the perturbations δf ss khij , can be chosen freely in (5.2). The choice which gives the most symmetric relations between the phase space distribution functions and the correlation function as well as between the different distribution function components, and leads to simplest source terms in the evolution equations is 11 With this choice the phase space distributions are also correctly normalised in the thermal limit. Note that despite the fact that the definition (5.3) involves the helicity projector, N ss khij does not depend on helicity. It is also symmetric in the energy and flavour indices. We can invert the parametrisation (5.2) to express the phase space distribution functions in terms of the matrix form of the correlation function. Using (5.3) this relation reads simply That is, with the normalisation (5.3), P s s khji is a "correctly normalised" projection operator in our basis.
The basis spanned by P s s khji can be used to define generalised distribution functions for any local correlation function. Below we need the adiabatic distribution functions, which can now be defined analogously to (5.4): f ss ad,khij ≡ tr P s s khji S < ad,kij (t, t) . This shows that the parametrisation (5.2) with the normalisation (5.3) naturally matches to the Fermi-Dirac distribution in the free theory (up to a sign for negative frequency states).

Generalised density matrix equation
Using the parametrisation (5.2) with the normalisation (5.3) for δS < k (t, t) in equation (4.39) we can derive an equation for the non-equilibrium distribution functions δf ss k,hij (t). The calculation is complicated by the fact that also the normalisation factors and the energy 11 In the massless limit (5.3) becomes singular. This is a technical problem similar to the one encountered with massless spinors, and it can be avoided by using a different normalisation. Alternatively one can use (5.3) with finite masses, and take the limit mi → 0 at the end of the calculation when needed. projection operators depend on time due to the time dependent masses m i (t), but the final master equation is structurally very simple: Hereṁ i ≡ ∂ t m i and we combined f k = f ad,k + δf k in the second line, with the adiabatic distribution functions f ad,k defined in equation (5.5), and the collision term is given by C sηs khilj ≡ i tr P s s khji Σ r eq,eff,kil (t, t) P ηs khlj . (5.8) Equation (5.7) describes all particle-antiparticle and flavour-coherence effects in the local limit and includes helicity. This universality is reflected in the large number of indices, which may appear overwhelming at first. However, all terms in (5.7) have simple interpretations. For example the first term on the right-hand side corresponds to the Hamiltonian commutator term in equation (4.39). It falls into this simple form because H ki P s ki = P s ki H ki = sω ki P s ki . The Hamiltonian term causes oscillations in the off-diagonal components and its simple form is pivotal in distinguishing the relevant time scales. For s = −s the oscillations are very fast. These oscillations are essential e.g. for vacuum particle production, but they can be a problem if one is interested only in flavour oscillations, which correspond to s = s but i = j, and usually evolve much more slowly. In the next section we derive an effective equation for flavour oscillations by averaging over the fast oscillations.
The second line in equation (5.7) arises from the time dependence of the basis matrices P ss khij in the parametrisation (5.2). The precise form of these terms depends on the normalisation, and the choice (5.3) turns out to give the most compact form. These terms are again essential for vacuum particle production, but they can be neglected in the leptogenesis application. The source terms ∂ t f ss ad,khij result from the projection of the adiabatic matrix source in (4.39). Lastly, the collision terms C sηs khilj can in general be separated into dispersive and absorptive parts, just by using Σ r = Σ H − iΣ A , and consequently we define The dispersive term C H,sηs khilj can be broadly characterised as a generalised matter Hamiltonian. It is of course a very general structure, but e.g. in the case of ordinary light neutrino mixing it can be shown [94] to exactly reproduce the neutrino effective potential in matter. It could be interesting for resonant leptogenesis as well, because it would make the energy splitting ∆ω k between the neutrinos a dynamical quantity. We will not consider the dispersive corrections numerically in this paper, but also this topic will be pursued elsewhere. Finally, the absorptive part C A,sηs khilj of the collision term contains both flavour-diagonal and offdiagonal damping rates, and importantly for leptogenesis, cross couplings between the diagonal and off-diagonal distribution functions. These cross coupling functions together with the CP-violating flavour oscillation are the mechanism which generates the lepton asymmetry.

Lepton source and washout terms
It is easy to write also the CP-violating source term and the washout terms of the lepton equation using the parametrisation (5.2) and normalisation (5.3) for the Majorana neutrino correlator. For the source term (4.41) and the washout term (4.42) we get, respectively, For the washout term (4.43), we again expand the adiabatic Wightman functions to leading order to get the O(y 2 i ) result (5.12) Expanded forms of equations (5.10)-(5.12) are given in appendix F, where we use the leading order approximation (4.27) for the effective self-energies and perform the traces.

Mass shell equations
As stated above, the most important benefit of the parametrisation (5.2) is that it allows to separate the physics corresponding to different time scales. In particular we can distinguish the mass shell functions corresponding to s = s (but including the flavour coherences i = j) from the fast oscillating coherence shell functions for which s = s . 12 For a graphical illustration see figure 2 in [77]. We denote these functions by Indeed, from equation (5.7) we see that for δf c,± khij the Hamiltonian term is proportional to ∓i(ω ki + ω kj ) = ∓2iω kij , corresponding to very fast particle-antiparticle oscillations (zitterbewegung). For the mass shell solutions δf m,± khij the Hamiltonian term is proportional to ∓i(ω ki − ω kj ) = ∓i∆ω kij , corresponding to flavour oscillations for i = j, at a frequency which is suppressed for large |k| or a small mass difference |m i − m j |. This is the case of interest for resonant leptogenesis.
If particle-antiparticle oscillations are much faster than the flavour oscillations, we can drop the coherence functions δf c k in equations (5.7), since their effect on the flavour oscillations averages out. 13 This results in a much simpler coarse-grained master equation 12 This naming scheme follows the earlier cQPA notation [77] and the one we already used in (3.13), although in our current treatment the phase space structures are not a priori restricted to a spectral form. 13 We can formally justify this as follows. First, generically δf c± are some functions that vary only in the flavour scale. Next, take the convolution of (5.7) with some appropriate normalised weight function W , e.g. the Weierstrass transform with W (t, t ) ∼ exp(−(t − t ) 2 /2σ 2 ), where for the mass shell functions: Removing the coherence shell solutions greatly facilitates the numerical solution, in particular for quasi-degenerate Majorana neutrinos. In addition, a number of relations hold between different components of δf k , which further simplify the numerical task; we will give these relations in appendix D. In the following sections we will solve equations (5.15) numerically using the leading order expansion (4.27) for the effective self-energies in the collision terms functions. Detailed expressions for C sss khilj are given in appendix E.
Hierarchical limit. Let us briefly comment on the validity of our master equations in the hierarchical limit of leptogenesis [20,42], where ∆m/m ∼ O(1). The general master equation (5.7) is of course applicable also in this limit. However, the condition ∆ω kij 2ω kij might hold only to a limited degree (see footnote 13), so that neglecting the coherence shell functions δf c k might not be warranted, possibly reducing the accuracy of the mass-shell equations (5.15). However, using equation (5.7) to model this case would be numerically very challenging, because due to the large hierarchy H m, both oscillation time scales are very fast compared to the heavy neutrino decoupling time. Luckily, due to the very same reason, one can in this case work in the decoupling limit (see section 8 below) and derive an effective Boltzmann equation for the lepton asymmetry. To our knowledge the contribution from the coherence shell functions δf c k has never been included, however. This could be done by starting from the master equation (5.7) and it would indeed be interesting to study the size of these corrections.

Expansion of the universe
So far all our equations have been formulated in the flat Minkowski spacetime, but we eventually need to work in the expanding Friedmann-Lemaître-Robertson-Walker (FLRW) background. In this section we show how the generalisation to an expanding universe can we can choose 1/∆ω k σ 1/2ω k . This has no effect on the mass-shell contributions, because they do not vary significantly over the time σ. However, the terms involving the coherence solutions behave as where D ± k represents whatever term that is multiplying the coherence solution. Because ωσ 1, these terms are extremely suppressed and completely drop from the averaged mass-shell equations. be done by a simple reinterpretation of all variables in the comoving frame. We finish this section by rewriting our master equations explicitly in the expanding universe in terms of the scaled inverse temperature.

Lagrangian in curved spacetime
First we need to generalise the Minkowski Lagrangian (4.1) to curved spacetime: Here g ≡ det(g µν ) is the determinant of the metric and ∇ µ is the covariant derivative given by the spin connection. We also added the non-minimal coupling of the Higgs field φ to the scalar curvature R and the factor √ −g originates from the volume form of the curved spacetime action integral. We then assume a spatially flat FLRW metric where a = a(η) is the dimensionless scale factor and the conformal time η is defined by dt = a(η) dη. With the metric (6.2) the contracted covariant derivative in equation (6.1) becomes where ∇ is the usual flat spatial derivative vector, ∂ 0 = ∂/∂η and a ≡ ∂a/∂η. Using equations (6.2) and (6.3) with √ −g = a 4 and scaling all fermion fields ψ according to ψ → a −3/2 ψ and the Higgs field as φ → a −1 φ, the Lagrangian (6.1) is transformed to where the scalar curvature is given by R = 6a /a 3 . From the Lagrangian (6.4) we see that the only effects of the expanding universe compared to the Minkowski theory (4.1) are that the time variable is replaced by the conformal time η, spatial coordinates become the comoving ones, neutrino masses are multiplied by the scale factor a and the Higgs field gets a geometric mass term. This mass term vanishes for a conformal coupling ξ = 1/6 or when the universe is radiation-dominated: a(η) ∝ η, which is the case to a high accuracy in leptogenesis. We shall thus continue to assume that the Higgs field is massless.
The comoving frame. Based on the above, all expressions and equations in the earlier sections given in the Minkowski background remain valid also in the expanding flat spacetime when we make the replacements where k com ≡ ka is the comoving momentum, T com ≡ T a is the comoving temperature of the relativistic SM heat bath and m i are the physical constant masses. Note that the phase space distribution functions f are dimensionless scalars and have the same values in both physical and comoving variables. We assume that the universe is dominated by relativistic SM particles, which are kept in equilibrium by fast gauge interactions, and that the universe expands adiabatically. In the absence of entropy production the physical temperature then scales as T ∝ a −1 , so that T com remains constant. The comoving momentum k com is also constant, as the physical momentum redshifts as k ∝ a −1 . We can also write the Hubble parameter as where m Pl is the Planck mass and g * is the effective number of relativistic degrees of freedom, which is g * = 110.25 when including all SM fields and two massless Majorana neutrinos. Note that the scale factor can now be written as a(η) = H(T com )η. The entropy density is given by s = 2π 2 g * T 3 /45. Overall, the leptogenesis equations retain the same form they had when formulated with a generic time dependent mass term in the Minkowski background, when using the replacements (6.5) and scaling the equations by appropriate powers of the scale factor a.

Final master equation in expanding space time
For the numerical implementation it is convenient to formulate the equations using dimensionless variables. The most relevant temperature scale for leptogenesis is around T = m 1 , where m 1 is the mass of the lightest Majorana neutrino. We thereby use the variables where H 1 ≡ H(m 1 ) is the Hubble parameter (6.6) evaluated at T = m 1 . The z-parameter is directly proportional to the scale factor so it serves as the time evolution parameter. Due to the constraints (D.9) there are only four independent mass shell functions for the two mixing Majorana neutrinos. We choose them as follows: δf kh ≡ δf m,+ kh11 , δf m,+ kh22 , Re δf m,+ kh12 , Im δf m,+ kh12 . (6.9) We can now formulate the k-dependent neutrino equation (5.15) as a vector equation for the components (6.9). Including also the lepton asymmetry equation (4.40), our final master equations written using the dimensionless variables are dδf kh dz = ∆ω k −C kh δf kh − df ad,kh dz , (6.10) where Y L ≡ (n − n )/s, and we used (4.34) to relate the lepton chemical potential to the asymmetry. The dimensionless tree level oscillation coefficient ∆ω k and the collision term coefficientC kh of the neutrino equation, as well as the CP-violating lepton source termS CP and the washout term coefficientsW (with W = δW, W ad ) are given by (6.13) Equations (6.10) and (6.11) are formally analogous to the momentum dependent Boltzmann equations, which we present in equations (G.5) and (G.6) of appendix G. The difference is that the quantities ∆ω k and C kh in equations (6.10) and (6.12) are 4 × 4 matrices consisting of the coefficients for different components of the column vector δf kh . To avoid confusion with earlier notation, we denote these components (when needed) with bracketed indices: for example δf kh [3] = Re δf m,+ kh12 , and C kh [12] is the coefficient of δf kh [2] in the equation of δf kh [1] . The matrix ∆ω k corresponds to the tree level flavour oscillation term and its only non-vanishing elements are ∆ω k [34] = −∆ω k [43] = ω k1 − ω k2 . The collision term coefficient matrix is given by where C + ilj ≡ C +++ khilj were defined in equation (5.8). Explicit expressions for the collision terms C +++ khilj are given in appendix E, computed using self-energies in the leading order expansion (4.27) and explicitly written in appendix B. Finally, explicit expressions for the lepton source and washout terms can be found in appendix F.

Numerical results
In this section we solve numerically the system of equations (6.10) and (6.11) for the Majorana neutrino distribution functions δf kh and the (normalised) lepton asymmetry density Y L in the case of two Majorana neutrinos and one lepton flavour. We start by considering possible initial conditions. Then, before going to the discussion of the final lepton asymmetry and its dependence on model parameters, we establish the time scales relevant for the problem and study the momentum dependent neutrino distribution functions.

Initial conditions
While we took t in → −∞ for the CTP to calculate the interactions, we can of course choose any finite time t = t 0 (or z = z 0 ) as the starting point of our calculation with arbitrary initial conditions for the correlation functions. In particular, we will consider both vacuum and thermal initial conditions for the Majorana neutrinos, while we assume that the lepton and Higgs distributions stay in local equilibrium in the common SM plasma temperature T . In both cases we assume a vanishing initial lepton asymmetry: Y L (t 0 ) = 0.
Vacuum initial conditions. Here we assume that at t = t 0 the Majorana neutrinos are decoupled and thus effectively in zero temperature (e.g. if only the SM particles get reheated after inflation). The full local neutrino correlator is then given by which is calculated from (4.23) by taking the limit T → 0. For the non-equilibrium part we then get, using equations (3.3) and (7.1), Note that this perturbation is not necessarily small as it measures the deviation of the full correlator S < (now initially in vacuum with zero temperature for the Majorana species) from the adiabatic equilibrium correlator S < ad (with the high temperature T of the SM particle species). For the mass shell functions the initial condition (7.2) with (7.1) becomes where we used equations (5.4) and (5.6).
Thermal initial conditions. Here we assume that also the Majorana neutrinos are in local thermal equilibrium with the SM particles at t = t 0 , corresponding to some high initial temperature T 0 m i . The full local neutrino correlator is then given by which trivially implies that δS < k (t 0 , t 0 ) = 0 and that all non-equilibrium distribution functions vanish initially: δf m,± khij (t 0 ) = 0. Note also that in this case the Majorana neutrinos deviate from equilibrium only due to the dynamical source term −∂ t f ad .

Physical scales and parameters
The minimal leptogenesis mechanism is mainly controlled by four time scales, corresponding to the expansion rate of the universe H, the Majorana neutrino decay rate Γ, and the oscillation frequencies ∆ω kij and 2ω kij . We already discussed the role of the flavour and zitterbewegung oscillations in section 5.4, finding that the latter may be relevant in the hierarchical limit, but can be safely ignored when ∆ω kij 2ω kij . We work in this limit here, but the question of the relative sizes of the three slow time scales H −1 , Γ −1 and ∆ω −1 kij still remains. The resonant leptogenesis mechanism turns out to be most efficient when all these scales are roughly comparable.
Indeed, if H ∆ω kij , there is no time for flavour oscillations to develop, while for H ∆ω kij the source terms ∼ H become suppressed. Also, in the very strong washout limit, Γ H the asymmetry is suppressed by thermalisation due to interactions, whereas for Γ ∆ω kij , the magnitude of CP-violation is suppressed, since the source in the asymmetry equation is proportional to the couplings y i . Additionally, if Γ ∆ω kij the flavour oscillations become over-damped [69], analogously to light mixing neutrino systems [3,4,7]. This leaves us to seek parameters in the range Γ ∼ ∆ω kij ∼ H for maximal resonant enhancement.
Benchmark parameters. Our simple leptogenesis model has 5 physical input parameters: the magnitudes of the two Yukawa couplings y 1 and y 2 , their relative CP-violating phase θ 12 , the lighter Majorana neutrino mass m 1 and the relative mass-squared difference ∆m 2 21 /m 2 1 ≡ (m 2 2 − m 2 1 )/m 2 1 . We define the Yukawa phase as We use the following set of benchmark values as a baseline: In addition one has to define the number of effective relativistic degrees of freedom g * and the initial temperature T 0 where we start the calculation. For these we use g * = 110 and z 0 = 10 −2 . For the benchmark parameters this corresponds to the initial temperature T 0 = m 1 /z 0 = 10 15 GeV and the Hubble expansion rate H 1 ≡ H(m 1 ) ≈ 1.428 × 10 8 GeV. Some more comments are still in order. First, we have chosen a high mass scale 10 13 GeV for the lightest Majorana neutrino, typical for traditional thermal leptogenesis [20]. However, it can be seen that all terms in equations (6.12) and (6.13) scale either as ∼ Γ/H 1 or ∼ ∆m/H 1 , so that the dynamics does not depend on the mass scale m 1 as long as Γ/H 1 is kept constant (and if ∆m ∼ Γ). More precisely, the dynamics can be effectively characterised by three parameters: the washout strength parameters and the number of flavour oscillations in a Hubble time: N 12 ≡ |∆m 12 |/(2πH 1 ). Using Γ (0) i = |y i | 2 m 1 /(8π) (see appendix G.1), we find that in the benchmark case K i ∼ O(10) (corresponding to strong washout) along with N 12 ≈ 1.5.
These estimates agree qualitatively with what is previously known from the semiclassical Boltzmann approach, where the CP-asymmetry parameter ε CP i ∝ sin(2θ 12 ) is resonantly enhanced for ∆m ∼ Γ. The CP-violating angle θ 12 was chosen to be maximal in this sense in the benchmark case, but it can be used to adjust the value of the final asymmetry downwards, as it affects the results mainly as an overall scaling factor with only a small impact on the dynamics.

Neutrino distribution functions
Solving the master equations (6.10) accurately requires of the order of hundred discrete momentum variables. Thousands of collision integrals C +++ khilj are then needed at each timestep, each of which contains a one-dimensional integral. It is clear that these cannot be feasibly computed during the evaluation. Fortunately, to the order we are working, they can be computed and fitted before solving (6.10). Moreover, the source terms in the neutrino ad,kh11 (bottom). On the right: the collision term coefficient functionsC kh [11] , which is the same for both helicities (top), and the helicity-odd combination (C k,+1 [14] −C k,−1 [14] )/2 (bottom). equation (6.10) are localised in momenta and temperature, which facilitates the fitting process.
In the left panels of figure 4 we show heat maps of the flavour diagonal adiabatic distribution function (5.6) (top) multiplied by a phase space factor, κ 2 f m,+ ad,kh11 , and the corresponding source term (bottom) −κ 2 df m,+ ad,kh11 / dz as a functions of z = m 1 /T and κ = |k|/T . Both functions are indeed localised around |k| T in momentum and the source term is also localised around z 1, while the distribution function f FD (ω k1 ) exhibits exponential fall off to zero in the same region. Plots for f m,+ ad,kh22 are qualitatively similar, with the fall off region moving according to the value of m 2 .
In the right panels of figure 4 we show heat maps ofC kh [11] = 2 Re(C +++ kh111 ) z/H 1 and (C k,+1 [14] −C k,−1 [14] )/2 = Im(C +++ k,+1,121 − C +++ k,−1,121 ) z/H 1 , given by equations (6.12) and (6.14). The former is the damping rate for the flavour-diagonal distribution function δf kh11 , and the latter is the helicity-odd part of the function which couples the flavour diagonal distribution to the off-diagonal distribution Im(δf kh12 ). The diagonal damping rate is the same for both helicities in our leading order approximation. Note that the colour bar scales in the figure are logarithmic for both positive and negative directions separately, except for a small region around zero (up to one tick in both directions), where linear scaling is used. Figure 4 was created using the benchmark resonant leptogenesis parameters given in subsection 7.2. Note that we have dropped the dispersive self-energy Σ H everywhere, as was discussed earlier. All other components of C +++ khilj are qualitatively similar to the ones shown here. In particular, the helicity even parts of the off-diagonal coupling functions 2 Re(C +++ kh121 ) and 2 Re(C +++ kh211 ), which are important for leptogenesis, have similar forms to the diagonal damping rates, and in the quasi-degenerate case m 1 m 2 also their scales are similar. Also, when Σ H is dropped, the diagonal damping rate C kh[ii] is actually just the momentum and temperature dependent decay rate of the Majorana neutrinos. This is because the leading order Majorana neutrino self-energy Σ A used in this work corresponds only to the decay and inverse decay processes.
In figure 5 we show similar heat maps of some non-equilibrium distribution functions δf m,+ khij , obtained from a numerical solution of equation (6.10), using the same benchmark parameters as in figure 4. In the top row we show the helicity-even parts of the flavour-diagonal and off-diagonal functions: κ 2 (δf m,+ k,+1,11 + δf m,+ k,−1,11 )/2 (left panel) and κ 2 Im(δf m,+ k,+1,12 + δf m,+ k,−1,12 )/2 (right panel). In the bottom row we show the corresponding helicity-odd parts, κ 2 (δf m,+ k,+1,11 − δf m,+ k,−1,11 )/2 (left panel) and κ 2 Im(δf m,+ k,+1,12 − δf m,+ k,−1,12 )/2 (right panel). We will later see that the helicity-even part of Im(δf m,+ kh12 ) gives the main contribution to the CP-violating lepton source term. The results show that the off-diagonal components are localised around z 1 and complete approximately one period of flavour oscillation before being exponentially suppressed (the blue colour indicates negative and the red colour positive values). The localisation around |k| T is due to the phase space factor κ 2 together with the exponential suppression at high momenta. The negative values of the h-even diagonal distribution (top left panel), extending to very small z result from the vacuum initial conditions used here. The other independent components of δf m,+ khij are again qualitatively similar to the ones shown here.

Results on lepton asymmetry
Having now established the scales of the problem and that the numerical solution of the neutrino correlation function is under control, we turn to study the lepton asymmetry evolution. Because of the smallness of the lepton chemical potential µ , we can solve the equations (6.10) and (6.11) sequentially. The neutrino equation is solved first and its solutions δf kh are used to calculate the lepton source termS CP and the washout term coefficients δW andW ad according to equations (6.13), (F.1a), (F.2a) and (F.3). These are then used to solve the lepton equation (6.11). Note that equation (6.10) does not couple helicities or momenta so we can solve each mode δf kh separately. In practice we reformulate the neutrino equations in terms of the helicity-even and helicity-odd combinations (δf k,+ ± δf k,− )/2, which are more convenient for the calculation of S CP and δW , as described in appendix F.
Benchmark case. In section 7.3 we presented some intermediate results for the phase space functions δf kh using the vacuum initial condition (7.3) and the benchmark parameter values (7.7). In figure 6 we show the source termS CP and the washout term coefficients W ad , δW of the lepton equation (6.11) as functions of z for the same parameters and initial conditions. We find that dividing the dimensionless momentum κ = |k|/T to 100 bins in logarithmic scale between 10 −2 and 10 4 already ensures that the results are not sensitive to the cutoff or the discretisation. Indeed, one can see from figures 4 and 5 that the largest contributions come from the range κ ∈ [0.1, 10]. In figure 6 we also compare the results from our full quantum kinetic equations (QKEs) (6.10) and (6.11) to those following from the traditional semiclassical Boltzmann equations (BEs) given in (G.5) and (G.6) and the corresponding momentum integrated rate equations (REs) (G. 19) and (G.20). The latter two sets of equations are well known, but we provided them explicitly for completeness. Also, we wrote all equations using a similar notation, which greatly facilitates the comparisons. Both BEs and REs require an externally provided CP-violating parameter ε CP i . At this point we are using ε CP i corresponding to the "mixed" regulator (G.24) [83,84] (see also e.g. [95]). Other regulators will be discussed in detail below.
It turns out that both the BEs and the REs significantly overestimate the source term S CP in the relativistic region z 1. The QKE-source starts to grow and changes the sign later, but all sources start to converge for z 2. On the other hand, the washout terms W ad and δW have only minor differences. As expected,W ad is by far the dominant of the two. It also appears to be identical in the full QKE and in the BE approach, and indeed it is: this term originates from the flavour-diagonal equilibrium part of the Majorana neutrino distributions, which we have calculated to the zeroth order in gradients in the QKEs. The function δW , which is proportional to the non-equilibrium perturbation δf , is approximately two orders of magnitude smaller and could be neglected with practically no effect on the final asymmetry. Also, there is no corresponding function δW in the leading order rate equations (G. 19) and (G.20) because in the Boltzmann approach this part results from the Pauli blocking and stimulated emission factors, which are dropped from the rate equations.
In figure 7 we show the lepton asymmetries Y L as a function of z in the benchmark case. Left panel corresponds to the vacuum initial conditions used above. The asymmetry behaves as expected from the source termS CP shown in the left panel of figure 6. While the Y L (z)-evolution obtained using the BE or the RE deviate strongly from that found using the full QKE, the final asymmetries differ by less than a factor of 2. In the right panel we show the case with thermal initial conditions for the Majorana neutrinos. The early evolution of the asymmetry is of course very different from the vacuum case, but the final asymmetries are identical with both initial conditions. This behaviour is due to the strong washout assumed in the benchmark case (K 1 ≈ 10.0 and K 2 ≈ 27.9), which efficiently erases the early evolution of the lepton asymmetry. The final asymmetry then mostly depends on the source at the very end of the integration range, where theS CP computed using different methods were found to agree better.
In figure 8 we show for comparison similar plots in the weak washout case, with |y 1 | = 0.01, corresponding to K 1 ≈ 0.279. The left panel again corresponds to the vacuum and the right panel to thermal initial conditions. Now the early evolution deviates less in different approaches. Also, as expected, the final asymmetries are no longer the same for vacuum and thermal initial conditions. The differences between the BE and RE predictions and our QKE results for the final asymmetry also remain significant. These results show that the lepton asymmetry evolution and its asymptotic value depend in an essential way on the model parameters. Also, it is clear that to obtain accurate results, one should use the full QKEs instead of the less accurate BE-or RE-approaches.
Varying mass difference. In the top panels of figure 9 we show the asymptotic lepton asymmetry as a function of (m 2 2 − m 2 1 )/m 2 1 ≈ 2∆m 21 /m 1 , keeping other benchmark parameters fixed. We again compare the full QKE results with the BE and RE predictions, now for four different CP-asymmetry parameters ε CP i that have been discussed in the literature [61,62,66,70,95,96]. The different choices, which we denote as 'mixed', 'difference', 'sum' and 'effective' vary in how the resonance near m 1 = m 2 is regulated. Explicit forms for the ε CP i -parameter and the regulators are given in appendix G.3. All approaches give qualitatively similar ∆m 21 -dependence with a single maximum at ∆m 21 ∼ Γ. However, a more close look reveals significant quantitative differences between the full QKE results and the BE and RE approximations, as well as between using different regulators in the latter two approaches.
The location of the resonance peak is approximatively given by ∆m 2 2 | for the difference and sum regulators in the Boltzmann approaches. We have shown these locations in the top-right panel of figure 9, denoting them by ∆ max diff. and ∆ max sum . The BE and RE results using different regulators fall either above or below the correct QKE-result shown by the solid blue line, varying by an almost order of magnitude for ∆m 21 Γ. The effective sum regulator given in [70] (see equation (8.13) below) is designed to work in the strong washout case; it is thus not surprising that it works best in our benchmark case. On the other hand, for ∆m 21 > Γ, where we enter the rapid flavour oscillation regime, all results converge. This is expected, since the regulators become irrelevant in the CP-parameter (G.23) and the diagonal elements decouple from the off-diagonals in the QKEs in this limit (we will show this explicitly below). We also observe that the BEs always give a slightly lower final asymmetry than do the REs, as was also observed in [97]. The difference between the BE and RE results is smaller, however, than the difference arising from using different regulators and eventually the correct QKEs. That is, treating the quantum physics part of the problem correctly is more important than the momentum dependence of the phase space distributions.
The situation gets even more interesting when we begin to vary the couplings. In the bottom panels of figure 9 we show the results for hierarchical (left panel) and for almost degenerate (right panel) Yukawa couplings. In the left panel the washout is weaker, K 1 ≈ 0.279 and K 2 ≈ 27.9, whereas in the right panel the washout is strong, K 1 ≈ 22.6. The different CP-asymmetry regulators now lead to even more dispersion in the BE and the RE results when ∆m 21 Γ. In the hierarchical case using the mixed regulator leads to two peaks at ∆m 21 Γ 1 and ∆m 21 Γ 2 , corresponding to the different regulators used for the two Majorana neutrinos in this case. The full QKE result, again shown by the solid blue line, shows no such structure. Also the effective regulator is somewhat less accurate here. In the right panel, with almost degenerate Yukawas, the mixed and sum regulators are in better agreement with the QKEs, but the difference regulator has an extra spurious enhancement because the regulator vanishes and the CP-asymmetry is unbounded in the double limit m 2 → m 1 and Γ 2 → Γ 1 .
The main take-home message from this section is that the Boltzmann equation and the rate equation approaches are inaccurate and strongly sensitive to the choice of the regulator in the resonant and quasidegenerate region ∆m 21 Γ i . The Boltzmann approach reproduces the full QKE results accurately in all cases only in the region ∆m 21 > Γ i when the regulator is already negligible and one is approaching the hierarchical mass limit. Also, the most accurate regulator over varying coupling strengths is the sum regulator of [62]. We shall show below how the sum regulator indeed consistently emerges when we reduce the QKEs to BEs in the helicity-symmetric decoupling limit.
Flavour oscillation. Let us look more closely at the role of the flavour oscillations in resonant leptogenesis. In the left panel of figure 10 we show the evolution of the helicityeven part of the off-diagonal Majorana neutrino phase space function Im(δf m,+ kh12 ) with |k| = 3T and ∆m 21 = 0.01m 1 . This is a case with rapid flavour oscillations corresponding to N 12 ≈ 110. The modes of this δf -component in the range |k|/T ∼ 0.1-10 give the dominant contribution to the integrated lepton source term. In the right panel we plot the contributionS CP,k of the same mode to the lepton source term, normalised according toS CP ≡ d 3 k/(2π) 3S CP,k /T 3 . Both the mode and its contribution to the lepton source display a strong oscillation pattern with a quickly dying amplitude. This decay of the oscillations is precisely the reason for the emergence of the semiclassical limit (shown as the green dash-dotted line) from the QKEs. We will make this explicit in equation (8.8) below. We show also in the right panel the full integrated source term from the QKEs (dotted line) and from the Boltzmann approach with the sum regulator (dash-dot-dotted line). Even in the QKE-result all oscillations are smoothed out in the integrated source, due to the phase differences between different modes.  Figure 10. Flavour oscillation of Im(δf m,+ k,h12 )(z) (helicity-even component) andS CP,k (z) for a single k-mode, with ∆m 21 = 0.01m 1 . Other parameters have the benchmark values, and vacuum initial conditions were used. On the right we show the single k-mode as well as the integrated lepton asymmetry source from the full QKEs and the BEs.

Helicity-symmetric approximation
In this section we will derive a series of approximations to the QKEs (5.15) (or equivalently (6.10)), eventually reducing them to the semiclassical Boltzmann limit. This process also leads to a simplified source term in the asymmetry equation (4.40), which eventually reproduces the CP-asymmetry parameter with the sum regulator.
Looking more closely at equations (4.41) and (F.1a), one can see that the lepton asymmetry is sourced mainly by the helicity-even combination of the imaginary part of the off-diagonal function δf kh12 in the non-relativistic or mildly relativistic case. Based on this observation, we drop the tracking of the helicity asymmetry in the equations for δf khij . We can then write a simpler set of equations for the h-even part, which we denote simply by δf kij henceforth, and a simpler form for the lepton source term S CP including only Im(δf k12 ). We will also work with vacuum dispersion relations, setting Σ H ≡ 0. The resulting helicity-symmetric equations are and where Γ k12 ≡ (Γ k11 + Γ k22 )/2 and Γ kil ≡ 2 Re(C +++ khilj )| Σ H →0 with C +++ khilj given by equation (E.1). Note that all Γ kil are now real so the diagonal functions δf kii couple directly only to Re(δf k12 ). The diagonal damping rate admits the factorisation Γ kii = (m i /ω ki )Γ (0) i X ki where m i /ω ki is the time dilation factor, Γ (0) i = |y i | 2 m i /(8π) is the tree-level vacuum decay width of the Majorana neutrino and X ki is the thermal quantum statistical correction factor, which obeys X ki → 1 when |k|/T → ∞ (see equation (G.7)). The damping rate Γ k12 in the off-diagonal equation, given by the average of the diagonal rates, agrees with the flavour coherence damping rate found in the density matrix formalism for mixing neutrinos [3,4,7]. This is an expected result, as the two phenomena are of course closely related. Indeed, if we further assume that Γ k21 Γ k12 , we can write equations (8.1)-(8.3) in a simple density matrix form: where H k ≡ diag(ω k1 , ω k2 ). A similar equation has been used earlier to study light neutrino mixing in the early universe [3,[6][7][8] and in the resonant leptogenesis context it was first derived in [61]. 14 The helicity-symmetric equations (8.1)-(8.4) provide an excellent approximation to the full QKEs with helicity-even initial conditions, as can be seen in figure 11 where we compare the two for our benchmark parameters. This is so because the neutrino source terms are helicity-symmetric and the helicity-asymmetry is only generated by loop effects. Also, it is evident from (F.1a) that the contribution to the source from the helicity-odd δf m k -function is thermally suppressed compared to the one coming from helicity-even part. This shows that resonant leptogenesis is dominated by the helicity-independent flavour mixing, although this conclusion partly relies on the fact that the asymmetry is mostly generated in the non-relativistic regime T m 1 . Helicity would play a more important role if the Majorana neutrinos were relativistic when the asymmetry is generated, like in some low scale leptogenesis scenarios [86]. Indeed, equations similar to (8.5), but keeping the helicity degree of freedom, were recently used to study leptogenesis [69,86,98]. Our full QKEs (5.7) as well as the mass-shell equations (5.15) are more general than (8.5) and the helicity dependent equations in [86,98,99]. Note in particular, that for more than two Majorana neutrino flavours the collision terms in (5.15) cannot in general be reduced to the canonical form for a density matrix equation.
Decoupling limit. Equations (8.1)-(8.4) still incorporate all essential flavour mixing consistent with full resummation of the interaction terms. Now we simplify these equations further in the case where the flavour oscillations are fast (∆m 21 Γ). We use the same reasoning as in section 5.4 (see footnote 13) to argue for dropping the flavour off-diagonal terms in the diagonal equations (8.1) and (8.2). We call this approximation the decoupling limit. The diagonal equations, written with the expansion of the universe, are then identical to the semiclassical Boltzmann equation (G.5). The washout term can also be approximated with the Boltzmann version or even with only the W ad contribution (G.10). Because the flavour-diagonal functions now decouple from the off-diagonal ones, they can be solved independently and their solutions can be treated as external sources to the off-diagonal function δf k12 .
Assuming that initially δf k12 (t 0 ) = 0 (a non-zero initial value could be easily added as a special solution to the homogeneous equation), the off-diagonal differential equation (8.3)  Figure 11. Comparison of the full results (solid line) to the helicity-symmetric approximation (dotted), decoupling limit (dash-dotted) and Boltzmann results with the sum regulator (dashed). We show the results for the lepton asymmetry source term (left) and the lepton asymmetry itself (right) as functions of z. The benchmark parameters with vacuum initial conditions were used.
can be integrated to give Substituting this into equation (8.4) we then get a closed formula for S CP which now defines the source term in the semiclassical Boltzmann equations; indeed comparing to the Boltzmann formula (G.8), this improved form shows that the combination Im(y * 1 y 2 ) Im(δf k12 ) acts as an effective dynamical CP-asymmetry parameter.
In the quasidegenerate case m 1 m 2 (i.e. at this point we assume the weakly resonant regime Γ ∆m 21 m 1 ) we can further approximate Γ kij Re(y * i y j )Γ kjj /|y j | 2 and Γ k11 /|y 1 | 2 Γ k22 /|y 2 | 2 . Then we can write the lepton source term into an even more suggestive form: This result shows that the lepton asymmetry is cumulatively sourced by the non-equilibrium perturbations in the diagonal mass-shell functions δf ii (u), such that past contributions are suppressed by the flavour coherence damping rate Γ 12 . We show the approximation (8.7) for S CP (with the diagonal solutions δf ii calculated from the ordinary Boltzmann equations) and the resulting lepton asymmetry by the red dash-dotted lines in figure 11. Even though in the figure we used the benchmark parameters where ∆m 21 Γ, equation (8.7) is still a relatively good approximation to the full QKEs, in particular at early times z 1. However, when z 1 it deviates from the full result, and at late times z 1 the asymmetry coincides perfectly with the Boltzmann result (see below) instead, shown by the dashed line in figure 11.

Boltzmann limit
To get an even closer comparison to the existing literature, we now make stronger approximations to evaluate the time integrals in equation (8.7) analytically. Indeed, if one assumes that the source functions are roughly constant, one can take them outside the u-integral. 15 If one further assumes that Γ k12 and ∆ω k12 are constants (we take X ki → 1 in Γ kii , assuming Γ kii (m i /ω ki )Γ (0) i , and neglect the Hubble expansion for the mass and momentum), one can perform the integrals to get A similar expression was also found in [57], but without the exponential damping factor in the oscillating term. This is important, because now we see that taking the limit t 0 → −∞, the oscillating part is damped to zero. 16 We already saw this effect in figure 10: while the source was there computed using the full QKE, the strongly damped rapid oscillation we observed corresponds to the second term in equation (8.8).
Substituting (8.8) back to equation (8.7) and dropping the damped oscillating term, one finds This form already greatly resembles the Boltzmann result (G.8) with the constant CPasymmetry parameter (G.23). It is now also clear that the CP-asymmetry is regulated by the coherence damping rate Γ 12 in the degenerate limit m 2 → m 1 . We also note that dropping the damped oscillating term is the main reason why the BE approach tends to initially overestimate the asymmetry. This effect is clearly visible in figures 6, 7 and 11. We can go even further and extract the CP-asymmetry parameter by using again the quasidegeneracy m 2 m 1 , whereby ∆ω k12 ∆m 2 12 /(2ω ki ), and taking into account the time dilation factor in the damping rate Γ kii = (m i /ω ki )Γ (0) i . The result (8.9) can then be written as

10)
15 This is consistent if the damping time is much shorter than the time of variation of the diagonal δfii-functions, which is given by the Hubble time. That is, an approximation of the type of equation (8.9) is valid in the strong damping limit: Γ12 H. 16 Authors of ref. [57] argued that the second term in (8.8) vanishes due to averaging out over oscillations, but this is not the correct explanation; note that (8.8) is valid even in the limit ∆ω12 → 0. Also, a similar expression and including the damping, albeit with a different damping factor (Γ11 − Γ22)/2 (giving rise to the 'difference' regulator), was found in [61].
Equation (8.10) is still different from the Boltzmann result (G.8) in its dependence on the flavour-diagonal functions δf ii . But using the quasidegeneracy argument once more to approximate (m 1 /ω k1 )X k1 (m 2 /ω k2 )X k2 , it can finally be written exactly in the same form as (G.8).
Strong washout limit. In articles [64,66,70] the strong washout limit was considered, where one assumes that the diagonal rate parameters are large separately. In this limit one can find an approximative late-time solution by putting all derivative terms to zero on the left-hand side of equations (8.1)-(8.3). One can then solve all distribution functions algebraically, and eventually the CP-violating source term (8.4) takes the form No other approximations have yet been made at this point. In the quasidegenerate limit this source then has the same form as (8.9) except that the off-diagonal backreaction to the diagonal functions is taken into account by a modification of the regulator term Γ 12 . 17 This approximation can be again reduced to the Boltzmann equation with yet another effective CP-asymmetry parameter, similar to equation (8.11), but with Here we used the fact that in the single lepton flavour and quasidegeneracy limit Γ k12 Γ k21 cos 2 θ 12 Γ k11 Γ k22 , where the Yukawa phase θ 12 was defined in equation (7.6). The result (8.13) agrees with the effective sum regulator defined in [70]. We have shown the results using this effective sum regulator in figure 9 as the black dashed line. This is indeed the best approximation in the strong washout limit. On the other hand, this approximation does not work well outside the strong washout case. We show in figure 12 the case where both Yukawa couplings are small, corresponding to the washout strength parameters K 1 ≈ 0.28 and K 2 ≈ 1.1. In this case the regulator (8.13) is worse than the simpler sum-regulator we found in the decoupling limit. These examples show that using different approximations one can accommodate the most relevant quantum effects in different parametric regions. However, no approximation remains quantitatively accurate throughout the parameter space.  Summary. We have used a series of controlled approximations to reduce the Majorana neutrino QKEs and the lepton asymmetry source term to the Boltzmann limit. The decoupling limit CP-asymmetry parameter (8.11) corresponds to the sum regulator (G. 26) in agreement with [62][63][64]. The previous work used some extra assumptions about the model parameters however, and their validity in the doubly degenerate limit have been questioned [96]. Our derivation is very different and does not rely on these assumptions. We also find a clear origin and interpretation for the regulator, corresponding to the flavour coherence damping in the neutrino equation. Our numerical work confirmed that the sum-regulator is indeed the most accurate in the Boltzmann approach in the sense that it provided qualitatively best results in all regimes, as can be seen in figures 9 and 12, while the effective sum regulator (8.13) is the most accurate one in the late time limit in the strong washout case. As a byproduct, we demonstrated quantitatively how the quantum oscillations are suppressed and the semiclassical limit arises in the hierarchical strongly damped limit ∆m 21 Γ 12 H. Beyond this case, the BE approach does not provide highly accurate results and QKEs are needed. However, we found that the helicity-symmetric QKEs give a very accurate final asymmetry throughout the resonant leptogenesis parameter range.

Comparison to earlier work
Leptogenesis has been studied extensively before , and many of the results presented here have been found in some form previously. In this section we will provide a more in-depth comparison of our results and other studies based on first-principles CTP methods, which we believe are most closely similar to ours.
In several studies of leptogenesis based on the CTP method (e.g. [46,62,65,68,100]) the non-equilibrium part of the Majorana propagator is identified with a homogeneous transient, while the inhomogeneous part is taken to be the thermal equilibrium solution.
In this approach the only non-equilibrium source is in the initial conditions, as there is no dynamical source for the asymmetry generation. Such approaches have been used to model the lepton asymmetry generation during the initial approach of the Majorana neutrinos to equilibrium [46,100]. Our formalism contains this effect, which shows up as the initial negative Y L dip in the left-hand side panels in figures 7 and 8. However, as these figures show, a moderate washout can erase this asymmetry. This restricts the use of pure transient methods to the weak washout case.
A very careful analysis using the transient approach was given in [62]. In particular off-diagonal pole propagators, which are crucial for the evolution of the non-equilibrium initial state in the two-time approach, were calculated in detail. These results are similar to our (4.17a). However, in our approach where the local correlator is evolved dynamically, it is sufficient to compute the pole propagators to a leading order approximation. In [62] it was also found that the sum regulator (G.26) is a reasonable choice for the CP-asymmetry parameter (G.23) in the Boltzmann approach. The Hubble expansion was not included in [62], but the issue was later addressed in [63,64]. Our derivation is more general, and does not impose restrictions on Yukawa couplings [62,63] or on the size of the deviation from equilibrium [64]. Our derivation also reveals the physical origin of the sum regulator as corresponding to the flavour coherence damping rate.
The approach in [61] (see also [66,69]) is similar to ours, but it relies heavily on the Wigner space representation. The neutrino correlator S < k is also expanded in a different basis, used earlier in the EWBG context [12,13] and in the cQPA approach [72][73][74][75][76][77][78]. In this basis the division of S < k into components with a characteristic time dependence is obscured, making it difficult to separate the particle-antiparticle and flavour coherence effects. Ref. [61] then used several approximations in integrating the particle-antiparticle coherences, in reduction to the spectral shell limit and a restriction to the quasidegenerate case, that we do not need to make. On the balance, the final QKE of [61] agrees with our equation (8.5) in the quasidegenerate limit (and includes the small backreaction to the neutrino equation from the lepton chemical potential, which we omitted.) Overall our derivation is more general and displays a hierarchy of QKEs, which separate the different physical scales. We also provide detailed numerical examples and comparisons, with the Hubble expansion included.
In [67] resonant leptogenesis was studied in the interaction picture method [101], using the double momentum representation, and it was found that the lepton asymmetry source term contains two distinct contributions from mixing and oscillations. 18 These results suggest that the usual Boltzmann approach, which contains only the mixing contribution, potentially captures only half of the actual late time lepton asymmetry. The findings of [67] were supported by an analysis [68] performed in a simplified non-equilibrium setup in the weak-washout regime. Both articles found that the mixing contribution is mainly due to the flavour diagonal functions and the oscillation contribution mainly due to the off-diagonal functions in the Majorana neutrino correlator. Ref. [68] found also another contribution resulting from the interference of the mixing and oscillation terms, which tends to cancel the other contributions in some cases. These results have been further discussed in [55].
Our results do not support these findings. Indeed, our result for the lepton asymmetry converges to the usual Boltzmann result in the weakly resonant case (see figure 9) and a difference by a factor of two should be clearly visible. It is then curious to note that, similarly to refs. [67,68], our lepton asymmetry source (F.1a) also contains two distinct contributions from the flavour diagonal and off-diagonal phase space functions. However, in our case the flavour-diagonal contribution is helicity suppressed as explained in section 8. This helicity dependence may thus be the source of the discrepancy, especially given that refs. [67,68] are based on scalar toy models where the effective Majorana states do not have Dirac structure. Also, while the results of [67] were essentially reproduced with true Majorana neutrinos in [96], that analysis was based on a semiclassical approach, which again may not necessarily implement the correct helicity structure.
Finally, flavour coherent equations similar to ones used to study light neutrino mixing in the early universe in [2][3][4][5][6][7][8], were used to study leptogenesis in [86,98,99]. We showed how these equations arise from our more general formalism. We also showed that a further reduction to even simpler, but very accurate helicity-even equations (8.5) is possible in the resonant leptogenesis case. Our formalism is also self-contained giving explicit expressions for all self-energies involved.

Conclusions and outlook
We have developed a general and comprehensive formalism for problems involving quantum coherence effects in spatially homogeneous and isotropic systems with mixing fermions. Our methods are applicable to various problems ranging from vacuum particle production to neutrino physics and resonant leptogenesis, which we used as an example and a platform for the more detailed developments of the method. In particular we concentrated on a benchmark model with two Majorana neutrinos and one lepton flavour including decay and inverse decay interactions.
Our method uses the CTP formalism [30,31] and the 2PIEA methods [34,35]. An essential part in our derivation of tractable quantum kinetic equations including coherence was finding a closed equation for the local correlation function S < k (t, t). This required a method to evaluate collision integrals which contain the full correlation function S < k (t 1 , t 2 ). We explained how this can be done when the system has dissipation. Two key elements in the process were the identification of proper adiabatic background solutions and the ansatz (3.8) to parametrise the perturbation δS < k (t 1 , t 2 ) in terms of the local correlation function in the collision terms. However, no assumptions were made to restrict the spin or the flavour structure of the local correlator, which makes it well suited for studying dynamical mixing. Another essential element was the use of the projector basis (5.2), which provides a clean separation of physics related to different time scales, e.g. the particle-antiparticle oscillations and the flavour oscillations. Our formalism can be seen as a generalisation of the cQPA method developed in [72][73][74][75][76][77][78] and further studied in [79].
Our main results include the quantum kinetic equation (5.7) which contains complete coherence information including particle-antiparticle mixing, and its coarse-grained version (5.15), which completely incorporates the flavour mixing in both particle and antiparticle sectors separately. Note that these QKEs cannot in general be written as a traditional density matrix equation, if there are more than two flavours, due to the complicated flavour structure of the collision term (5.8). For the two-flavour case we derived even more simplified but very accurate helicity-symmetric QKEs (8.1)-(8.4) and further wrote them into a density matrix form (8.5) in the nearly degenerate limit m 1 m 2 . Eventually we reduced our QKEs to the diagonal Boltzmann limit, deriving the CP-asymmetry parameter ε CP i of leptogenesis with the 'sum'-regulator (8.11) directly from the quantum transport formalism. We also pointed out that the sum-regulator physically corresponds to the coherence damping rate of the Majorana neutrinos in the underlying QKEs.
The question of the correct CP-asymmetry regulator in the semiclassical approach has been under some discussion recently [55,62,66,70]. We performed careful numerical comparisons of our QKEs and the Boltzmann equations endowed with different choices for ε CP i . We found that the sum-regulator derived here in the decoupling limit, and first found in [62], agrees best qualitatively with the QKEs throughout the parameter range of interest for resonant leptogenesis. We also implemented the integrated Boltzmann rate equations and compared them with the momentum dependent BEs and the QKEs. It turns out that the difference between the REs and BEs was always much smaller than the difference stemming from the choice of different regulators. While the sum-regulator was generically the best choice for the BE and RE approach, their results can still be wrong by a factor ∼ 2-4 for ∆m 21 Γ. In the strong washout case the modified sum regulator (8.13) of [66,70] gives even more accurate late-time results. For high accuracy results throughout the parameter range, QKEs are needed however. We found that the helicity-symmetric QKEs (8.2)-(8.4) give a very accurate approximation to the full flavour QKEs (5.15) for all parameters. Finally, all of these approaches are in good agreement in the weakly resonant case, m i ∆m 21 Γ, as expected. Leptogenesis has been studied extensively using the CTP approach  and many of the results shown here have been found previously. We believe that our treatment stands out in displaying the most complete set of quantum kinetic equations, with a clean separation of different physics and by giving a comprehensive account of the approximations made in deriving them. Based on our results one can easily compute the effective self-energy functions to different levels of approximation in the coupling constant expansion, and including also the coherent propagators in the internal lines. In some accounts we did less than what has been done before. Definitely the phenomenological reach of our results is compromised by our not including the 2 → 2 scattering processes or multiple lepton flavours in our equations. We will leave these to a future work.
In this article we mainly concentrated on resonant leptogenesis, but our methods are applicable as such, or easily modifiable to other versions of the leptogenesis mechanism [20,21]. We already briefly discussed thermal leptogenesis in the hierarchical limit in section 5.4. In this case it would be interesting to compute corrections to the lepton asymmetry source arising from the mixed particle-antiparticle flavour correlation functions (δf c,± kh12 in the notation of section 5.4), starting from the full QKEs (5.7) and working in the decoupling limit. Our equations can also easily accommodate dispersive corrections to the neutrino QKEs. Such corrections would generalise the vacuum Hamiltonian to include the matter effects, similar to the well known case with light neutrinos. This would replace the mass difference ∆m 21 with an effective dynamical quantity and potentially change the quantitative predictions in resonant leptogenesis.

Note added
The Mathematica code package that was used to compute all numerical results in this paper is publicly available at https://doi.org/10.5281/zenodo.5025929.
Here we present the derivation of the results (2.13) and (2.14), starting from the Schwinger-Dyson equations (2.10). We will shorten the notation explained below equations (2.11) even further, by leaving out the momentum indices and convolution signs: e.g. S p 0,k * Σ p k * S p k → S p 0 Σ p S p . It is also important to note that the inverse free propagator S −1 0 contains a derivative operator, and thus changing the direction of its operation generates additional surface terms.
First, the pole equation (2.10a) for p = r, a can be formally iterated as This is consistent with the Hermiticity properties (2.8) of the pole propagators, which relate equation (A.1) for p = r to equation (A.3) for p = a and vice versa. Equation (2.10b) for S < (similarly for S > ) can now be iterated and rearranged as follows: The second and third equalities explicitly show the first and second iteration of the first equation (A.4). In the fourth step this iteration is assumed to continue indefinitely. To get the final two lines we then used equations (A.2) and (A.3). Note that the final line can also be written in the form where the vacuum part was given (using dimensional regularisation) in equation (4.47) and the temperature dependent part is The temperature dependent integral (B.3) has been worked out in [102]. We further parametrise the integrals (B.1b), (B.1c) and (B.3) with and a α k 0 , |k| = Tã α k 0 /T, |k|/T , whereã α ,b α are dimensionless functions, and α = A, H, <, > and S = S eq , δS. After some calculation we get the following results (with κ 0 ≡ k 0 /T , κ ≡ |k|/T ): The integrals given above forã A eq ,b A eq , δã ≶ and δb ≶ may be further calculated in closed form using logarithm and dilogarithm functions (like in [42]). The results given here hold only for |κ 0 | > κ > 0 (timelike four-momentum). The corresponding results for κ > |κ 0 | > 0 (spacelike four-momentum) are attained by replacing the integral operator above as follows: Also, in the results (B.7) it was assumed that the energy parameter κ 0 is real. If these results are continued to complex values of κ 0 (as required when considering finite widths) then the implicit sign functions in the absolute values must be applied to the real parts only (e.g. sgn(κ 0 )κ 0 = |κ 0 | should be replaced by sgn(Re κ 0 )κ 0 = √ κ 0 2 ). Finally, the a α and b α functions have the following k 0 -symmetry properties:

C Adiabatic pole propagator inversion
In this paper we have used the leading order spectral approximation for the adiabatic solutions, which indeed is a good approximation in the weak coupling limit. In this appendix we show how to obtain more complete solutions for the adiabatic pole propagators S r,a ad (k, t). We start by writing equation (4.17a) in an equivalent form and with explicit flavour indices: Generally, the inverse of a 2 × 2 block matrix A ij , satisfying 2 l=1 A il S lj = δ ij 1, is given by S ij = (A ji − A jk A −1 lk A li ) −1 with k = i and l = j. For given i and j the indices k and l are fixed, because the block dimension is only 2. This block-wise inversion formula can be generalised to larger block matrices, which would indeed be necessary if there were more flavours. We now solve the flavour components of the pole propagators from (C.1) as follows (equivalent formulae were given in [83]): The inverses here are taken with respect to Dirac indices only and we suppressed the (k, t)arguments of the self-energies. These solutions can also be obtained formally by expanding the inverse of / k − m − Σ p eq as a geometric series in powers of the off-diagonal self-energy Σ p eq,ij (with i = j) and performing a resummation of the series.
The solutions (C.2) and (C.3) are still general. We now use the leading order equilibrium self-energy Σ eq given by equations (4.16a) and (4.16b). We neglect here the vacuum part of Σ H eq which acquires additional Dirac and flavour structures from the vacuum on-shell renormalisation and would make the resulting formulae below more complicated. However, near the poles the vacuum part should only give small corrections. Hence, we use here the pole self-energies Σ p eq,ij (k) = c w y i y * j P L + y * i y j P R γ µ S p eq,µ (k).

(C.4)
Using this form, the Dirac matrix inverses in equations (C.2) and (C.3) can be written explicitly as with the definitions Note that since D p ij = D p ji , the denominator of all flavour components S p ad,ij is the same, D p 12 . Using the formulae presented in this appendix, one can straightforwardly implement more accurate approximations for the adiabatic correlation functions (4.20) and the effective self-energies (4.26).

D Independent components of the neutrino Wightman function
The Hermiticity property (2.8a), the sum rule (2.9) and the Majorana condition (4.4) can be used to derive relations among the components of the neutrino Wightman functions. First, the Majorana condition (4.4) and equation (2.4) imply that S > k (t 1 , t 2 ) = −CS < −k (t 2 , t 1 ) T C −1 in the two-time representation. We can then write the aforementioned three relations as where we also used 2A = S > +S < in the sum rule. These identities hold for the full Majorana Wightman functions S <,> kij (t, t). Assuming that they hold independently for the adiabatic functions S <,> ad,kij (t, t) implies that the non-equilibrium local correlator δS < kij (t, t) satisfies the constraints δS < kij (t, t) = δS < kji (t, t) † , (D.4) δS < kij (t, t) = −γ 0 C δS < −k,ji (t, t) T C −1 γ 0 . In the case of two neutrino flavours (i, j = 1, 2) these relations imply that for a given helicity h the sixteen different (s, s , i, j) components of δf kh can be reduced to only six independent distributions with 10 degrees of freedom, given explicitly by the following equations: Note that these relations are satisfied for perturbations when the sum rule (D.2) is satisfied by the full adiabatic solution. This is true e.g. for a full thermal solution and for the free particle solution (4.23). One can still use approximative forms for the adiabatic solutions in various expressions for the sources and collision terms.

E Neutrino collision term traces
Here we give explicit results for the collision term trace functions C sss khilj which we use in the mass shell equation (5.15) of the Majorana neutrinos. We actually need only the s = +1 component C +++ khilj because of the relations (D.9) and because we only solve the positive energy solutions. Using the definition (5.8) and the leading order expansion (4.27) for the effective self-energy, we get first C +++ khilj = i tr P ++ khji Σ r eq,il (k +l )P ++ khlj . Here the self-energy is given in the Wigner representation and evaluated at the four-momentum (k +l ) µ ≡ (+ω kl , k).
Next we observe that in the case with two flavours (i, l, j = 1, 2) P ss khlj P ss khji = P ss khli , which implies that C +++ khilj = i tr P ++ khli Σ r eq,il (k +l ) . Notably, the collision function is in this specific case independent of the last flavour index j. Using equations (4.16) and (4.51) for the equilibrium Majorana self-energy functions with Σ r = γ 0 Σ r and Σ r = Σ H − iΣ A , and performing the Dirac matrix trace yields the result C +++ khilj = c w N m kil 2ω ki ω kl Re(y * i y l ) m i k +l + m l k +i µ − ih Im(y * i y l ) m i k ⊥ +l + m l k ⊥ +i µ × S A eq (k +l ) + iS H(T ) eq (k +l ) µ . (E.1) The self-energy four-vector functions S µ are calculated in appendix B. We also defined N m kij ≡ N ss khij (see equation (5.3)) and (k ⊥ ) µ ≡ (|k|, k 0k ).

F Lepton CP-source and washout terms
Explicit forms for the CP-source and washout terms were given in equations (5.10)-(5.12). We now substitute the self-energy functions defined in (4.16a) and (4.16c) to these expressions and use the symmetry properties (B.10b) and (B.10c). We also use the lowest order result (4.27) for the effective self-energy, adapted for Σ A L ≡ P R Σ A and the constraints (D.9) for the non-equilibrium distribution functions. We keep the coherence shell functions here for completeness, but split the results to separate mass and coherence shell parts, S CP = S m CP + S c CP and δW = δW m + δW c . After performing the traces, we get Here µ is the lepton chemical potential, β = 1/T , and we further defined (k ⊥ ) µ ≡ (|k|, k 0k ).
The self-energy functions S A eq,µ and δS A µ = (δS > µ + δS < µ )/2 are calculated in section B. We also defined here N m kij ≡ N ss khij and N c kij ≡ N s,−s khij using equation (5.3). Finally, the adiabatic washout term W ad can be simplified to The results (F.1a) and (F.2a) show explicitly how the helicity sums of the non-equilibrium distribution functions δf kh enter the leading order CP-source and washout terms. In particular, from equation (F.1a) we see that CP-violation is only sourced by the helicityodd combinations Re(δf k,+1 − δf k,−1 ) of the real parts and helicity-even combinations Im(δf k,+1 + δf k,−1 ) of the imaginary parts.

G Semiclassical Boltzmann equations
For comparison with our main quantum kinetic equations, we implement and solve numerically the semiclassical Boltzmann equations (see e.g. [84,88,97,103,104]) and the momentum integrated rate equations in our model (4.1). We include only the decay and inverse decay contributions supplemented by a RIS-correction term required to cure the problem of spurious equilibrium source [88,97,104]. We write the equations first following ref. [103], to establish the notation and to facilitate comparison with the literature. We then present the equations in a more compact form which is directly comparable to our main equations and also more suitable for a numerical implementation.

G.1 Momentum-dependent equations
As in the main text, we assume that SM-fields are in kinetic equilibrium and that Higgs field chemical potential vanishes, µ φ = 0. However, we make no assumption about the form of the phase space distribution functions f i of the Majorana neutrinos i = 1, 2, and consider full quantum statistics with Pauli blocking and stimulated emission factors. In the expanding universe, the Boltzmann equations for f i and the lepton asymmetry n L can then be written as Here f x (p x , t) is the phase space distribution function, ω x ≡ |p x | 2 + m 2 x is the on-shell energy and dπ x ≡ d 3 p x /(2π) 3 /(2ω x ) is the phase space integration element for particle species x = i, , φ. The RIS-subtraction has been performed as in [88,97,104], which ensures that also the term associated with the CP-asymmetry parameter ε CP i in equation (G.2) vanishes in thermal equilibrium. For the number densities and distribution functions we are using the following notations: where f FD and f BE are the Fermi-Dirac and Bose-Einstein distribution functions (4.12). The factor c w = 2 is the SM doublet multiplicity and g i = 2 is the number of spin (or helicity) states of the Majorana neutrino N i . Finally, the tree level matrix element for the neutrino decay process, summed over the Majorana neutrino spins and the SM doublet multiplicity is: where m = m φ = 0 was assumed in the end. The corresponding total decay width is then Γ (0) i ≡ |M i | 2 /(16πg i m i ) = |y i | 2 m i /(8π). In the following we use these tree level results for the leading approximation.
We now give equations (G.1) and (G.2) in a more compact dimensionless form similar to (6.10) and (6.11), in terms of the variable z = m 1 /T . We also omit the lepton backreaction term (proportional to ε CP i f − ) in equation (G.1), as we did in (6.10). Working consistently to linear order in perturbations we then find where Y L = n L /s, s is the entropy density, δf ki ≡ f ki − f eq ki and f eq ki = f FD (ω ki ). The dimensionless neutrino collision term coefficient, CP-violating lepton source term and the lepton washout term coefficients are here given bỹ with H 1 ≡ H(m 1 ) and X ki = T |k| log e (ω ki +|k|)/T − 1 e ω ki /T − e |k|/T , (G.11) For the total z-derivative of f eq ki we use df eq ki / dz = f eq ki (f eq ki − 1)m 2 i /(m 1 ω ki ), assuming that |k|/T is independent of z. Note also that n LWad H 1 /z is exactly equal to the W ad defined by (F.3), when we use n L ≈ c w T 2 µ /6.

G.2 Rate equations
Simplified rate equations for the number densities can be derived from the Boltzmann equations (G.1) and (G.2) with two additional assumptions: kinetic approximation for the Majorana neutrino distributions and Maxwell-Boltzmann statistics for all particle species. To this end, we use f i = (n i /n eq i )f eq i , replace f FD and f BE by f MB (k 0 ) ≡ e −k 0 /T and remove all extra terms originating from the Pauli blocking and stimulated emission factors. Integrating equations (G.1) and (G.2) over momenta then gives: Here K n are the modified Bessel functions of the second kind. Here we kept also the lepton backreaction term which is the last term on the right-hand side of equation (G.14). Again, we can write equations (G.14) and (G.15) in a compact dimensionless form: The lepton source term in these equations is given byS CP = i ε CP i D i δY i and for the washout term δW = 0 andW ad = − i D i Y eq i /(2Y eq ). We used these equations also to check that the lepton backreaction was indeed numerically negligible in all examples (with vanishing initial lepton asymmetry) that we studied.

G.3 CP-asymmetry parameter
The Boltzmann equations (G.1), (G.2), (G.5), (G.6), (G.14), (G.15), (G. 19) and (G.20) given above do not depend on the exact form of the CP-asymmetry parameter ε CP i . There we used the generic definition (in the unflavoured case) where Γ(N i → φ) + Γ(N i → φ) = Γ i is the total (vacuum) decay width of the Majorana neutrino to the lepton and Higgs doublets. The CP-asymmetry parameter vanishes at tree-level and is calculated from the interference of tree level and higher order amplitudes. This is highly nontrivial, as the calculation breaks down when using ordinary perturbation theory in the degenerate limit m 2 → m 1 .
We only consider the self-energy contribution (also called indirect or ε-type CP-violation), neglecting the vertex correction (i.e. the direct ε -type CP-violation), and use the generic form ε CP i,x = Im (y * 1 y 2 ) 2 |y 1 | 2 |y 2 | 2 (m 2 2 − m 2 1 )m i Γ This result is specific to the case of two Majorana neutrinos and one lepton flavour. Here R ij,x is the regulator which removes the singularity that would occur when m 2 → m 1 . We consider four alternative regulators [62,70,83,84,105] (the subscript x): the 'mixed' regulator, the 'difference' regulator, the 'sum' regulator and the 'effective' sum regulator given by The relative phase θ ij of the Yukawa couplings was defined in equation (7.6). Note that ε CP i,diff with the difference regulator is still singular if |y 1 | = |y 2 |. This means that it is unbound and leads to unphysical results when approaching the doubly degenerate limit |y 2 | → |y 1 | and m 2 → m 1 . For more discussion about the validity of different regulators, see [96].