Kadanoff-Baym approach to the thermal resonant leptogenesis

Using the non-equilibrium Green function method (Kadanoff-Baym equations) in the expanding universe, we investigate evolution of the lepton number asymmetry when the right-handed (RH) neutrinos have almost degenerate masses $|M_i^2-M_j^2| \ll M_i^2$. The resonantly enhanced $CP$-violating parameter $\varepsilon_i$ associated with the decay of the RH neutrino $N_i$ is obtained. It is proportional to an enhancement factor $(M_i^2-M_j^2) M_i \Gamma_j/ ((M_i^2-M_j^2)^2 +R_{ij}^2)$ with the regulator $R_{ij}=M_i \Gamma_i+M_j \Gamma_j$. The result is consistent with the previous result obtained by Garny et al., in a constant background with an out-of-equilibrium initial state. We discuss the origin of such a regulator, and why it is not like $R_{ij}=M_i \Gamma_i-M_j \Gamma_j$.


Introduction
Origin of the Baryon asymmetry in the universe is one of the unsolved issues in particle physics. Although the standard model (SM) satisfies the Sakharov's three conditions [1], sufficient number of baryon asymmetry cannot be produced due to the smallness of the CP -asymmetry in the CKM matrices and the modest electroweek phase transition. On the other hand, the neutrino oscillation which implies tiny neutrino masses demands that some extension of the SM is necessary. Introducing right-handed (RH) neutrinos N i with large Majorana masses M i gives a natural solution to explain the smallness of the neutrino masses via see-saw mechanism, but it also naturally explain the Baryon number asymmetry in the universe through the leptogenesis [2]. (See e.g., a very nice recent review [3].) In this scenario, RH neutrinos are produced thermally by the reheating after inflation. As temperature decreases with the expansion of the universe down to the Majorana mass scale, RH neutrinos become out of thermal equilibrium and their CP -asymmetric decay into the SM leptons and the Higgs produce lepton number asymmetry in the universe. The lepton number asymmetry produced is then converted into the baryon number asymmetry through the nonperturbative B + L -violating process of sphalerons in the SM [4].
If the Majorana masses of the RH neutrinos have a hierarchical structure, the lightest Majorana mass must satisfy the Davidson-Ibarra(DI) bound [5], M > ∼ 10 9 GeV in order to produce sufficient lepton number asymmetry. When at least two of the RH neutrinos are degenerate in their masses, the DI bound can be evaded. In this case, quantum oscillation of almost degenerate RH neutrinos resonantly enhance the CP -violating decay and hence lepton number asymmetry can be produced sufficiently even for RH neutrino masses as light as TeV scale. This scenario is known as the resonant leptogenesis [8] [9] [10]. Such light RH neutrinos might induce detectable non-unitarity of the mixing matrix of active neutrinos [11] [12] and have attracted much attention.
TeV scale leptogenesis has attracted enormous attention in light of the LHC experiment [13]- [40]. The scale can be made even smaller if the leptotenesis occurs through CP -violating oscillations between RH neutrinos far away from the thermal equilibrium. The mechanism plays an important role in the model of νMSM [41]- [45].
Furthermore, light RH neutrinos do not give large radiative corrections to the Higgs boson mass and are safe in view of the naturalness [46]. Related to the naturalness of the electroweak weak against higher physical scales, one of the authors proposed a classically conformal U (1) B−L extension of the SM [47] [48]. In this model, B − L gauge symmetry is spontaneously broken via the Coleman-Weinberg mechanism which triggers the electroweak gauge symmetry. In [49], we further showed that if the Higgs potential is flat at the Planck scale, the model naturally predicts the Higgs boson mass at around 126 GeV and TeV scale B − L breaking (or the leptogenesis). This motivated us to investigate the TeV scale leptogenesis in the U (1) B−L model [28].
In the resonant case, the CP -asymmetry in the decay of N i mainly comes from an interference of the tree and the self-energy one-loop diagrams (See Figure 1: Tree and one-loop diagrams of the RH neutrino decay. In the resonant case, an interference of the tree and the self-energy diagram [6] [7] gives a dominant contribution to the CP -violating parameter. where h is the neutrino Yukawa coupling and Γ i ≃ (h † h) ii M i /8π is the decay width of N i . The resonant enhancement of the CP -violating parameter was discussed in [50]. Systematic considerations were performed by Pilaftsis [9][51] [52], and he found that the regulator in the denominator is given by R ij = M i Γ j . If the mass difference is larger than the decay width, we have |M 2 i − M 2 j | ≫ R ij , and ε i is suppressed by Γ i /M ∼ (h † h) ii . However, in the degenerate case, |M i − M j | ∼ Γ and ε can be enhanced to O((h † h) 0 ) ∼ 1. Hence the determination of the regulator R ij is essential for a precise prediction of the lepton number asymmetry in the resonant leptogenesis. The authors [53] calculated the resummed propagator of the RH neutrinos and obtained a different regulator R ij = |M i Γ i − M j Γ j |. By using their result, the enhancement factor becomes much larger. The origin of the difference of the regulators is discussed in [54] [55]. Since the lower scale of the leptogenesis is strongly sensitive to the form of the regulator, it is very important to systematically evaluate the precise form of the regulator.
Conventionally, leptogenesis is often calculated based on the classical Boltzmann equation which describes the time evolution of the phase space distribution function of on-shell particles [85]. In the Boltzmann equation, the interactions between particles are taken into account through the collision terms that comprise the S-matrix elements calculated separately in the framework of quantum field theory. The authors [57] applied the non-equilibrium Green's function method with the Kadanoff-Baym (KB) equations developed in studies of the transport phenomena [87] [88] and derived the full-quantum evolution equation for the lepton number in the hierarchical mass case. Using this method, one can systematically take into account quantum interference, finite temperature and finite density effects. 1 The method was intensively used in the leptogenesis in various papers [58]- [69]. In the resonant leptogenesis, since the quantum interference effect is crucial to the evaluation of the CP -violating parameter, we can expect importance of such a full-quantum mechanical formulation based on the KB equations. In [70], the authors used the method to obtain an oscillating CP -violating parameter in the flat space-time. Then applying it to the Boltzmann equation in the expanding universe, they calculated the lepton number asymmetry. In the strong washout regime, the oscillation is averaged out and the lepton number asymmetry is expressed with an effective CP -violating parameter. Then the maximal value agrees with the case of R ij = M i Γ j [71]. The authors of [72] also found an oscillatory behavior by a different calculation, and discussed an implication to the leptogenesis in the expanding universe. The quantum oscillations in the flavored leptogenesis was also performed in [73][74] [75].
Recently Garny et al. [56] systematically investigated generation of the lepton asymmetry in the resonant leptogenesis using the formulas developed in [58] [59]. In the investigation, they considered a non-equilibrium initial condition in a time-independent background and calculated generation of the lepton number asymmetry. Starting from the vacuum initial state for the RH neutrinos, they read off the CP -violating parameter from the generated lepton number asymmetry. The effective regulator they derived is R ij = M i Γ i + M j Γ j , which differs from the previous results, R ij = M i Γ j by [9] or R ij = |M i Γ i − M j Γ j | by [53].
The purpose of the present paper is to perform systematic investigations of the thermal resonant leptogenesis based on the KB equations. We scrutinize various properties of the Green functions of the RH neutrinos, and directly extract the CP -violating parameter ε i from the evolution equation for the lepton number in the expanding universe, with an emphasis on the quantum flavor oscillations.
The paper is organized as follows. In section 2.1 and 2.2, we first summarize the basic properties of various Green functions and the Kadanoff-Baym (KB) equations that these Green functions must satisfy. Then we derive the evolution equation of the lepton number in the expanding universe in section 2.3. The evolution equation is written in terms of the propagators of the RH neutrinos, the SM leptons and the Higgs. In section 2.4 we explain how the KB equation is reduced to the ordinary Boltzmann equation. The most important ingredient necessary to solve the evolution equation for the lepton number is the Wightman functions of the RH neutrinos. The flavor diagonal component is directly related to the distribution function, but more important for the lepton asymmetry is its off-diagonal component.
In section 3, we investigate how the flavor oscillation affects the off-diagonal component of the propagators. In the section, we focus on the resonant oscillations in the thermal equilibrium. In section 3.1, 3.2 and 3.3, we study the properties of the retarded and advanced propagators in which information of the spectrum is encoded. Then we study the Wightman functions with information of the distribution functions.
In section 4, we scrutinize the behavior of Green functions out of equilibrium. In the expanding universe, Green functions are approximated in the leading order approximation by the thermal values at the local temperature. But in order to calculate the lepton asymmetry, deviations from the thermal values are important. We show in section 4.3 that the deviations of the flavor offdiagonal Wightman functions from the thermal values behave quite differently from behaviors of other Green functions.
In section 5, we apply the calculated deviations of the Wightman functions of the RH neutrinos into the evolution equation derived in section 2, and obtain the quantum Boltzmann equation for the lepton number asymmetry. The deviations are classified into 3 types. One of them generate the lepton number asymmetry while the other two wash out the generated asymmetry. In section 5.4, we read off the CP -violating parameter ε and show that the regulator is given by In section 6, we give a physical interpretation why the regulator R ij = M i Γ i +M j Γ j appears instead of R ij = M i Γ i −M j Γ j . In particular, we show that if we neglect what we call the off-shell contributions the regulator is erroneously given by In section 7, we summarize our results.
In appendix A and B, we give a brief introduction to the closed time path (CTP) formalism and the KB equations. In appendix C, we introduce the 2PI formalism and then in appendix D we derive the self-energies for the RH neutrinos and the SM leptons based on the 2PI formalism. In appendix E and F, some useful identities in calculating convolutions are given. From appendix G to J, we give details of the calculations of various Green functions. In appendix K, we give anther derivation of the off-diagonal component of the Wightman functions out of equilibrium. The calculation explains why the regulator R ij = M i Γ i + M j Γ j naturally appears. Appendices L and M are calculations of some equations in the paper.

Evolution equations of lepton numbers
A systematic method to investigate the evolution of lepton asymmetry is the Kadanoff-Baym (KB) equations. The advantage of the KB equation to the Boltzmann equation is that it gives a quantum evolution equation of various correlation functions which does not distinguish on-shell and off-shell states. Accordingly it can take into account quantum coherence of the system and memory effects. Also the doubling problem in the scattering processes with on-shell internal lines can be systematically resolved (see [68] and references therein). The KB equation can be reduced to the classical Boltzmann equation only in special cases where the memory effects can be neglected.
Time-evolution of a quantum system is determined by the Hamiltonian of the system and the initial wave function at the initial time t = t i . Such timeevolution is described by the wave function at later times, or instead, a set of all n-point Green functions. Of course, it is practically impossible to study the evolution equations containing all the n-point functions and we need to select an important set of observables. In the classical approach based on the Boltzmann equation, one-particle distribution function on the phase space is selected. In the quantum Boltzmann approach, two-point Green functions are selected.
In this section, we summarize notations of various Green functions and their basic properties in the thermal equilibrium. We also summarize the nonequilibrium evolution equation (KB equation) for the Green functions. More details are given in Appendices. After brief reviews in section 2.2 and 2.3, we derive the evolution equation of the lepton number in section 2.4 and 2.5.

Our model
The model we consider is an extension of the SM with RH neutrinos ν R,i . i is the flavor index, i = 1, 2, 3. We set N i = ν R,i + ν c R,i . The Lagrangian is given by where α, β = 1, 2, 3 and a, b = 1, 2 are flavor indices of the SM leptons ℓ α a and isospin SU (2) L indices respectively. M i is the Majorana mass of N i and h iα is the Yukawa coupling of N i , ℓ α a and the Higgs φ a doublet. P R(L) are chiral projections on right(left)-handed fermions. In the present paper, we consider the case of almost degenerate Majorana masses at TeV scale. Then the Yukawa couplings become very small h iα ≪ 1 so as to generate tiny neutrino masses through the see-saw mechanism. Hence the decay width Γ i ≃ (h † h) ii M i /8π is much smaller than the mass M i .

Green functions and KMS relations
Various Green functions are introduced in field theories. Consider a fermion field ψ. The statistical propagator G F and the spectral density G ρ are defined as The statistical propagator G F contains information of the particle density of the state on which operators are evaluated. On the other hand, the spectral density G ρ gives information of the particle's mass and decay width. Because of the anti-commutator, γ 0 G ρ (x 0 , y 0 ) becomes proportional to the spatial delta function δ 3 (x − y) at the equal time x 0 = y 0 : where 1 is an identity matrix in the flavor and the spinor indices.
Other useful Green functions are the Wightman functions and the retarded and advanced Green functions are given by The spectral function can be written as In the present paper, we assume homogeneity along the spatial directions so that we can always use the Fourier transform in the 3-dimensional space. If the state is described by the thermal equilibrium state, we can further Fourier transform in the time direction. 2 In the thermal equilibrium at temperature T , the Green functions G(x, y) are anti-periodic in the time direction with an imaginary period iβ = i/T and their Fourier transforms satisfy the KMS (Kubo Martin Schwinger) relation Here f (q 0 ) is the Fermi-Dirac distribution function f (q 0 ) = 1/(e q0/T + 1). In presence of the chemical potential µ, q 0 is replaced by q 0 − µ. Since the relation relates the fluctuation described by the Wightman function to the dissipation described by the retarded Green function, it is also called the fluctuationdissipation relation. By this relation, the spectrum of the system determines all the Green functions. When the system becomes out of equilibrium, the KMS relation is violated. The violation plays an important role in the leptogenesis. As a final remark in this section, let us recall that the explicit forms of the Wightman functions of free charged fermions (bosons) are given by where ω q is the energy of the on-shell particle, andĝ ± = (±ω q γ 0 − q · γ + m), η = +1 for Dirac fermions with their mass m andĝ ± = 1, η = −1 for bosons. f q andf q are distribution functions of on-shell particles and anti-particles with their spatial momentum q respectively. They are not necessarily the equilibrium distribution functions.

Kadanoff-Baym equations
If the system is out of equilibrium and the state is time-dependent, we cannot use the ordinary perturbative method based on the Feynmann propagators. A general formalism is given by the closed-time-path (CTP) formalism in which perturbative vertices are inserted on the closed-time-path C = C + + C − . See Appendix A for brief review and Figure 4 there. One of the self-consistent approximation of the Schwinger-Dyson equations in the CTP formalism is called Kadanoff-Baym (KB) equation. Derivation of the KB equations is given in Appendix B and C. The equations for the retarded and advanced Green functions are is the free kinetic operator whose derivatives act on a field at x. Π R/A is the self-energy and defined in eq.(B.21). They have the same properties as G R/A , e.g., Π R (x, y) = 0 for x 0 < y 0 is satisfied. Note that the integration range in (2.11) is constrained between x 0 and y 0 : because of the step functions in Π R/A and G R/A . Therefore G R/A (x, y) is determined by the local information between x 0 and y 0 . Namely G R/A does not depend on the information of the system in the past: there is no memory effect for G R/A . Other Green functions G * ( * = F, ρ, ≶) satisfy (2.13) In the second equality, we have used the properties of R/A functions. By using eq.(2.11), this equation can be solved formally in terms of the self-energy function and the R/A Green functions as (2.14) In the last line * -operation denotes the convolution operation. Let us see the memory effect of G * . Generally speaking, the integrals in (2.13) over z are performed from the past at the initial time t int to x 0 or y 0 . This makes Green functions dependent on the state of the system in the past before x 0 or y 0 . This is indeed the case for G F and G ≶ , but for the spectral density G ρ , there is no memory effect. It can be seen by using Π ρ = Π R − Π A . Then the integral of (2.13) can be rewritten as Or it can be directly seen from the relation In the thermal equilibrium, since the system is translationally invariant, (2.14) can be Fourier transformed and (2.16) These equations (2.11), (2.13) are not closed within the two-point Green functions because the self-energy Π contains n(> 2)-point functions. Hence, in order to solve them explicitly, we need to make an approximation to express n(> 2)-point functions in terms of the two-point functions. 2PI effective action method is one of the simplest and self-consistent methods. (See Appendix C for brief explanation.) By using it, the self-energies Π in the above equations (2.11) (2.13) are represented as a sum of 1PI diagrams made of full propagators, and consequently these equations can be interpreted as simultaneous equations for various propagators in the system. These self-consistent equations among the propagators are especially called the Kadanoff-Baym equations.

Evolution of lepton number in the expanding universe
Now we investigate the KB equations of lepton numbers in the expanding universe. We first define Green functions, G, S and ∆ for the RH neutrinos, the SM lepton doublet and the Higgs doublet respectively: The classical inverse propagators are given by In this paper, we consider the spatially flat space-time with the scale factor a(t): (2.23) We useμ,ν, . . . as the space-time indices and µ, ν, . . . as the local Lorentz indices. γ matrices are written as γμ(t) = γ µ eμ µ where the vier-bein field eμ µ satisfies eμ µ eν ν gμν = η µν . In the following we mainly use t-independent γ µ = (γ 0 , γ) instead of t-dependent γμ(t). The delta-function becomes δ g (x − y) = δ 4 (x − y)/a 3 (x 0 ). In the background, the covariant derivative (2.22) becomes Since the spin connection is given by Ωμ = aH[γ µ , γ 0 ]/4, the covariant derivative for spinors in (2.20), (2.21) is given by Here the Hubble parameter is defined by H(t) =ȧ/a. In the radiation dominant universe, it is given by Lepton number density n L is a matrix with flavor indices α, β and isospin indices a, b. It is given by theμ = 0 component of the lepton number current Here tr{· · · } is the trace of the spinors. Because of the spatial homogeneity, divergence of the current j L is equal to On the other hand, it can be rewritten as 3 In the second equality, we have used the definition of S −1 0 (x, z) in (2.21). 3 Here, we have defined the derivative operator ← − / ∇y as By using the KB equation of (2.13) for the SM lepton Green function S> < , we have where Σ is the self-energy of the SM lepton. The second equality is obtained by using the relations (B.9) and (B.10).
Acting iS −1 0 from the right, a similar equation can be derived: (2.31) By using these equation, (2.28) becomes This is the evolution equation for the lepton numbers in the expanding universe. The right hand side (r.h.s.) is written as an integral of the full propagator S of the SM lepton and its self-energy Σ. Since the self-energy Σ contains various diagrams, some systematic simplification of Σ is necessary for practical calculations. A well-known approach is to use the 2-particles-irreducible (2PI) formalism briefly reviewed in Appendix B. In the 2PI formalism, the self-energy diagrams are obtained by taking a variation of 2PI diagrams made of full propagators with respect to the full propagator.
In the leading approximation, the self-energy Σ is obtained from the simplest 2PI diagram of Figure 2. Note that each propagator represents a full propagator, and the self-energy of the SM lepton is obtained by cutting the propagator ℓ. The next simplest 2PI diagram is given by Figure 5 in Appendix D, but in most of the present analysis, we consider only the contribution from Figure 2. It gives a good approximation if the RH neutrinos have almost degenerate masses. The contribution to the lepton self-energy Σ from Figure 2 is written in terms of the full propagators: Recall that (i, j) are flavor indices of the RH neutrinos. Then, summing the 2) with Yukawa interactions. Each line represents a full propagator of the SM lepton, Higgs and the RH neutrino. By taking a functional derivative with respect to each propagator, we can obtain the self-energy for the corresponding particle.
lepton flavor α, β and SU (2) L isospin a, b indices, we have Here we used the fact that the electroweak symmetry is restored at the temperature T TeV we are in mind and hence the propagators are written in SU (2) symmetric forms: S αβ ab = S αβ δ ab , ∆ ab = ∆δ ab . g w = 2 is the number of d.o.f. of SU (2) L doublets. Since the third and the fourth terms are complex conjugate to the second and the first terms, we can simplify the above equation as where we have defined τ = z 0 and This is the equation we evaluate in the following investigations. As we mentioned above, the r.h.s. contains only the contribution from the simplest 2PI diagram of Figure 2. This corresponds to taking the processes of decay and inverse-decay of the RH neutrinos. The effects of scattering can be taken into account 4 by considering the next simplest diagram of Figure 5. A systematical study of the KB equation including the scattering effects is given in [68][69].

Boltzmann equation for the lepton number
The evolution equation (2.35) of the lepton number is determined by the behavior of full propagators of the RH neutrinos G, the SM leptons S and the Higgs ∆. In sections 3 and 4, we investigate detailed properties of the propagator G of the RH neutrinos. In this section, we will see how an ordinary Boltzmann-type equation can be derived from eq.(2.35) by using the quasi-particle approximation for the SM particles described by S and ∆. The quasi-particle approximation is an approximation to express the Green functions in terms of distribution functions of quasi-particles with a mass m and a width Γ. Hence the propagators in this approximation are obtained from the free Wightman function of eq.(2.10) by introducing the decay width Γ. For a moment, we neglect the time-dependence of the background. For the SM leptons, we have where ω p = m 2 ℓ + |p| 2 /a 2 and / p ± ≡ ±ω p γ 0 − p · γ/a. Here we assumed the flavor independence of the lepton propagators, S αβ ∝ δ αβ , for simplicity. 5 Similarly the Wightman functions of the Higgs boson becomes The thermal mass and width are given by m ℓ,φ ∼ gT , Γ ℓ,φ ∼ g 2 T where g is the SM gauge coupling g. The effects of the thermal plasma play very important roles, and are systematically investigated in [68] [69]. For example, the thermal mass of the Higgs becomes larger than the RH neutrino masses at very high temperature. In the present paper, we focus on the largeness of Γ ℓ,φ as an important thermal effect and do not consider other effects.
In these expressions we defined (−1) ǫ = ±1 for ǫ = ± respectively. The distribution functions are assumed to be in the kinematical equilibrium and given by the Fermi-Dirac or the Bose-Einstein distributions at temperature T with a chemical potential: For anti-particles, the signs of the chemical potentials are reversed and their distributions are given by In the second equalities of eq. (2.37) and (2.38), we have defined which satisfy Now we come back to the time-dependence of the background. Since the scale factor a(t) is time-dependent, temperature T , thermal mass and width are dependent on the time t and we need to specify at which time these quantities in the quasi-particle approximation of eq. (2.37) and (2.38) are defined. If the temperature of the universe is sufficiently low (e.g., ∼ 10 TeV), the decay width is much larger than the Hubble expansion rate: and the propagators damp quickly at |x 0 − y 0 | ≫ 1/Γ l,φ . For such short period, time-dependence of the physical quantities such as the scale factor in the propagators (2.37) and (2.38) are suppressed by H/Γ ℓ,φ , and we can approximate these quantities as being constant in the integration of τ in (2.35). Then the physical quantities can be evaluated at time t = X xy = (x 0 + y 0 )/2 as we see in (4.4). By Fourier transforming in the spatial direction and using the above approx-imation, (2.35) becomes where t = x 0 and d 3 q = d 3 q/a 3 (t). Using the quasi-particle approximations (2.37) and (2.38), π> < (τ, t; q) are given by . In order to solve the evolution equation (2.44), we need detailed information of the Wightman function G ij ≷ (x, y) of the RH neutrinos. In the following sections, we obtain behaviors of the Wightman Green functions, especially deviations from the thermal equilibrium values in the expanding universe.
Here we briefly comment on the basic structures of the r.h.s. First, contributions from the flavor diagonal part i = j are evaluated by the quasi-particle approximation. G ii ≶ is proportional to (f ǫ Ni ) or (1 − f ǫ Ni ) respectively where f Ni is the distribution function of the RH neutrino N i . Therefore, combined with the distribution functions from π ≷ , flavor diagonal term gives the tree-level decay or inverse-decay of Figure 1, and wash out the generated lepton number asymmetry 6 .
On the other hand, by using the formal solution of the Wightman function in (2.14), contributions from the flavor off-diagonal part i = j are interpreted as an interference effect between the tree and 1-loop diagrams as follows. The formal solution is written in terms of the self-energy Π ≷ as (2.48) In the leading order approximation with the 2PI diagram of Figure 2, the selfenergy Π ≷ is written as a functional of the full propagators of the SM lepton and the Higgs as in eq.(D.6). Hence it can be interpreted as an interaction vertex of ℓφ ↔ N at u ∼ v. The RH neutrino propagates from u ∼ v to another interaction vertex at x ∼ y. By inserting this expression into eq.(2.44) and taking the on-shell limit of the RH neutrinos, CP -asymmetric interference between the tree and the one-loop self-energy diagram can be obtained. If the RH neutrinos propagating between these vertices are off-shell, the contribution is interpreted as s-channel scattering processes. Hence flavor off-diagonal terms in the r.h.s. of (2.44) give both of the CP -asymmetric decay of the RH neutrinos and the washout of the lepton numbers via s-channel scattering of leptons and Higgs.
In the resonant leptogenesis where the RH neutrinos have almost degenerate masses, however, it is not legitimate to separate the on-shell and off-shell contributions as above since the RH neutrinos are coherently mixed between different flavors, as has been mentioned in [61] [68]. Therefore we need to scrutinize the behavior of the Wightman functions G ij ≷ (t, τ ; q) in the expanding universe.

Summary of this section
The evolution of the lepton number is given by (2.35) or its Fourier transform (2.44). They are the basic equations we evaluate in the following sections. If we adopt the quasi-particle approximations of (2.37) or (2.38), an ordinary classical Boltzmann equation is derived. But in the resonant leptogenesis, quantum coherence between different flavors of N i plays an essential role and such an approximation is not valid for the RH neutrinos. An evaluation of the r.h.s. of (2.35) by scrutinizing the behavior of the off-diagonal components of the Wightman functions G ij ≷ , which is formally solved as (2.48), in the expanding universe is the main issue in the following sections.

Resonant oscillation of RH neutrinos
In this section, we study how the RH neutrinos with almost degenerate masses behave in the thermal equilibrium. Deviation from the thermal equilibrium is investigated in the next section 4.
We consider two flavors i = 1, 2 whose masses are almost degenerate. The third flavor RH neutrino is assumed to have larger mass. In order to calculate the evolution of the lepton asymmetry in (2.35), we need to know the Wightman functions G ≷ of the RH neutrinos. And, since the KB equation of G ij ≷ is formally solved by the convolution eq.(2.48), it is necessary to investigate the properties of the retarded (advanced) Green functions G ij R/A first. We first study both of the flavor diagonal (i = j) and off-diagonal (i = j) components of G ij R in the equilibrium. Then we will see the behaviors of the Wightman functions G ij ≷ in the thermal equilibrium. Throughout the paper, G d (also Π d for the self-energy) and G ′ (Π ′ ) denote the flavor diagonal i = j and off-diagonal i = j components respectively: (3.1)

Retarded/Advanced propagators
We first define the spatial Fourier transform of G R/A by Similarly, for the self-energy, we define Then using (2.25), the KB equation (3.2) becomes This is the basic equation for G R/A . We then decompose the propagator and the self-energy into flavor diagonal and off-diagonal parts: Using this decomposition, we solve the KB equation (3.5) iteratively.
First we define the differential-integral operator D d x 0 by In terms of the operator, the flavor diagonal component of the KB equation (3.5) becomes Similarly the KB equation of the flavor off-diagonal component is written as We then introduce the kernel G with a retarded (advanced) boundary condition. Using G R/A , we can integrate the equations (3.8), (3.9) as Then we can iteratively solve the above equations by expanding it with respect to the small off-diagonal component of the Yukawa (3.14) Here * denotes a convolution in the time-direction. The second term G d (2) R/A in the flavor diagonal propagator (3.13) is the second order of (h † h) ′ and smaller than G R/A or G ′ R/A . Hence we drop it and write G d(0) as G d for notational simplicity in the following.
We note that the above integrals do not have the memory effect. This is because the convolution is written explicitly as, e.g., (3.15) and the integration region is limited between x 0 and y 0 . Namely, the retarded (advanced) propagators are "local" functions of time during x 0 and y 0 and insensitive to the past (t < x 0 , y 0 ). This is different from the convolution contained in the Wightman functions (2.48) in which the integration range of time is extended to the past.

Diagonal G d R/A in thermal equilibrium
We will first look at the flavor diagonal component of the propagator G d R/A (x 0 , y 0 ; q) in the thermal equilibrium at temperature T . The scale factor a is also fixed at a 0 = a(x 0 ) = a(y 0 ). Because of the translational invariance in the time direction, G R/A (x 0 , y 0 ; q) can be further Fourier transformed: Then the KB equation (3.11) becomes and can be solved The real part of the self-energy gives the mass and wave-function renormalization. In the following we assume that they are already taken into account in the bare Lagrangian and focus only on the imaginary part Π d . The one-loop diagonal self-energy in the thermal equilibrium is expressed as Π d ρ = γ µ Π d ρ,µ . From the imaginary part of the pole of the propagator G d(eq) R (q), we see that the decay width Γ q of the RH neutrino is given by where Ω ǫi ≡ ǫω iq − iΓ iq /2 (3.21) and In the real time representation, it becomes 7 Γ q is multiplied by the Lorentz boost factor as Γ q ≃ (M/ω q ) × Γ where Γ ≡ Γ q=0 is the decay width of the RH neutrino. In the present paper, we consider a situation that two RH neutrinos are almost degenerate in their masses In the following, we sometimes use the averages denoted by quantities without the flavor index i, j , Ω ǫ = Ω ǫi + Ω ǫj 2 , etc.

Off-diagonal G ′ R/A in thermal equilibrium
We then study the behavior of the flavor off-diagonal component G R/A of the retarded (advanced) propagators in the thermal equilibrium. From (3.14), it is given by The q 0 integration can be performed by summing residues of the poles. Eq. (3.20) shows that the retarded propagator G d(eq)ii R has poles at q 0 = Ω ±,i and the advanced propagator G d(eq)jj A has poles at q 0 = Ω * ±,j . The self-energy Π R/A consists of the SM lepton and the Higgs propagator, and hence it has poles at q 0 = ǫ ℓ ω p + ǫ φ ω k ∓ iΓ ℓφ /2 with a large imaginary part. Because of this, the residues of the poles of the self-energy are suppressed by Γ i /Γ ℓφ ≪ 1. Noting the relation we can see that the contribution ǫ = −ǫ ′ is also suppressed by ∆M/M compared to the ǫ = ǫ ′ contribution. Hence, dropping these suppressed contributions, we have .

(3.29)
We also used the approximation Π( The minus signs in the parentheses come from the relative minus sign of the residue in (3.27). Because of this, the off-diagonal Green functions vanish at x 0 = y 0 : This should generally hold by the definition of Note that the flavor off-diagonal components of the retarded (advanced) propagators are enhanced by the factor 1/(Ω i − Ω j ) (or its complex conjugate). Such a large enhancement comes from the large mixing of the RH neutrinos with almost degenerate masses. For the self-energies h (ǫω q ) = 0 as in Appendix D, the following expressions [56] are reproduced: Here we have used the relation which is valid for ω q ≃ ω iq ≃ ω jq and ω q Γ q ≃ M Γ. Hence, the enhancement fac- As shown in section 3.6, the same enhancement factor, that is, the same regulator appears in the off-diagonal Wightman function in the thermal equilibrium.
For the deviations of the off-diagonal Wightman functions out of equilibrium, however, we show in section 4.5 that the enhancement factor is changed to be 1/(Ω i − Ω * j ). This corresponds to the regulator (M i Γ i + M j Γ j ). Finally we note the validity of the expansion with respect to the off-diagonal components of the Yukawa couplings (h † h) ′ . From the expressions (3.32), the iterative expansions (3.13) and (3.14) turn out to be valid when the real part of the off-diagonal components of Yukawa coupling Hence the expansion is understood as an expansion of the ratio (

Wightman functions
The Wightman functions can be solved as (2.14) or (2.48). If we take terms up to the first order of (h † h) ′ , the flavor diagonal component is given by Similarly the flavor off-diagonal component is given by By using (3.14) and (3.34), (3.35) can be also rewritten as which makes it clear that the off-diagonal part of the self-energy causes the flavor mixing of the RH neutrino. 9

Diagonal Wightman G d ≷ in thermal equilibrium
In the thermal equilibrium, the Wightman function can be easily obtained by using the KMS relation.
can be written as Let f (q) be the thermal distribution function for the RH neutrinos. Note that f (q) is a function of q 0 , which is not equal to the on-shell energy ω q . The KMS relation for the self-energy function is Using the solution of the KB equation for the spectral density It is nothing but the KMS relation (2.9) for the Green function.
Performing the q 0 integration, it becomes Here we have dropped the contributions from poles of the distribution function f (q 0 ) since they are suppressed by Γ/T ≪ 1. Furthermore we used the distribution function by dropping the imaginary part of the pole Ω ǫi in f (q) because it is suppressed again by the factor Γ ≪ T . Recall that it satisfies the relation in the thermal equilibrium. The off-diagonal component also satisfies the KMS relation and we (3.42) for x 0 > y 0 . We have used similar approximations by dropping suppressed contributions by Γ/T and Γ/Γ ℓφ . The off-diagonal component of the thermal Wightman functions are enhanced by the same factor 1/(Ω ǫi −Ω ǫj ) as in (3.28). Hence the flavor oscillation of the Wightman function in the thermal equilibrium is enhanced by a factor with the regulator At the temperature T ≫ ∆M we have in mind, f i and f j can be almost identified.
The off-diagonal Wightman functions in the thermal equilibrium vanishes at the equal time x 0 = y 0 : Later this property becomes very important to evaluate the deviation of the off-diagonal component of the Wightman function when the system is out of thermal equilibrium.

Summary of this section
In this section, we calculated various propagators of the RH neutrinos in the thermal equilibrium. We especially focused on the resonant enhancement of the flavor oscillation of N i . Retarded or advanced propagators are composed of two propagating modes, i and j flavors. The flavor diagonal components are given by (3.20) or (3.23). Since their masses are almost degenerate, the flavor offdiagonal component is largely enhanced due to their oscillation as in (3.28) or (3.29). The enhancement factor is proportional to 1/(

Propagators out of equilibrium
Now we study effects of the expanding universe into account. First we summarize various time-scales in the system. An important time scale is given by the Hubble expansion rate H of the universe. Other scales are the decay widths of the SM particles Γ φ , Γ ℓ and of the RH neutrino Γ i . Another important time scale in the resonant leptogenesis is given by the mass difference ∆M of the RH neutrinos because it gives the frequency of the flavor oscillation.
In type I sea-saw model studied in the present paper, the decay width Γ i of the RH neutrino is approximately given by Γ i ∼ (h † h) ii M i /8π. The ratio of Γ i to the Hubble parameter (2.26) at temperature T = M i is rewritten in terms of the "effective neutrino mass"m i as (see e.g. [3]) where v is the scale of the EWSB. Hence if we take the Yukawa coupling so as tom i ∼ 0.1 eV, the ratio becomes K i ∼ 100. This corresponds to the strong washout regime. The Yukawa coupling itself is very small (h ∼ 10 −5 for M ∼ 10 TeV), and we have the inequalities

Deviation of self-energy from the thermal value
Under the condition (4.2), we can expand the scale factor as The other physical quantities such as temperature are correlated with the change of the scale factor, and can be similarly expanded.
In order to calculate the out-of-equilibrium behavior of various Green functions in the expanding universe, we need to evaluate the change of the selfenergies Π(x, y). The self-energy of the RH neutrino is a rapidly decreasing function with the relative time as ∼ e −Γ ℓφ (x 0 −y 0 ) due to the SM gauge interactions. So in the leading order approximation, the self-energy Π(x, y) can be evaluated by the thermal value with the local temperature at the center-of-mass time x 0 ∼ y 0 ∼ X xy . Therefore it is convenient to write the self-energy as The first equation of of (4.4) is the definition of Π(X; s; q). In the second equality, we replaced Π by its thermal value Π (eq) since the SM leptons and Higgs are in the thermal equilibrium and the self-energy of the RH neutrinos is well approximated by its thermal value. Π (eq) (X xy ; s) means the thermal self-energy in the thermal equilibrium evaluated at time X xy . In evaluating the Wightman function G ≷ of the RH neutrinos, we need to know a difference of the self-energy Π(u, v) from the thermal value at a later time t. For example, in (2.48), the difference of the self-energy Π(X uv ; s) at X uv and the thermal value Π (eq) (t; s) at t = X xy controls the behavior of G ij ≷ . In this case, the time difference between X uv and t = X xy is given by the inverse of the decay width Γ i of the RH neutrino N i . Since the derivative expansion of the self-energy around the thermal value is a good approximation: The second term is of order O(H/Γ i ) owing to (4.6). The third term comes from the chemical potential of leptons generated by CP -violating decay of the RH neutrinos. So it is the genuine deviation of the self-energy from the thermal value at the same time X uv . In this section, we mainly focus on the change of the physical quantities, namely the second term because the back reaction of the generated lepton asymmetry to the evolution of the number density of the RH neutrinos is very small. The effect of the chemical potential becomes important in the generation of the lepton asymmetry and is considered in section 5.

Notice for notations
As already used in (4.4), Π(X; s) is the self-energy at the center-of-mass time X with the relative time s. For the thermal value Π (eq) (X; s), X is not necessarily at the center-of-mass time, but, more generally, denotes the reference time when it is evaluated. s is always the relative time. For the thermal value, we also use its Fourier transform (4.8) In order to avoid complications of appearance, we use the same notations Π for Π(X; s) and its Fourier transform Π(X; q). They can be distinguished by their arguments, s or q, if necessary. We always use s for the relative time and q for its conjugate frequency. For the first argument (the reference time), we use X or t. The same notation is used for the thermal Green functions. We hope it does not cause any confusion to the readers.

Retarded propagator out of equilibrium ∆G R
First we study how the retarded (advanced) propagators of the RH neutrinos deviate from the thermal value in the expanding universe. Consider the flavor diagonal component G d R/A first. We write the deviation around the thermal value G d(eq) by ∆G d : Note that ∆G d R/A depends on the reference time t at which the equilibrium value is evaluated. It is calculated in Appendix G and given by (4.10) The first term is the change of the physical parameters such as mass or width in Ω ǫ and Z ǫ . The second term represents a change of the spinor structure due to an expansion of the universe in the propagator during the propagation. The retarded (advanced) propagator does not have the memory effect, and the deviation is essentially determined by the change of the local temperature. By taking a variation of (3.14), the deviation of the off-diagonal components G ′ R/A can be expressed in terms of the deviation of the diagonal components The above formula is used to evaluate the deviation of the Wightman functions of the RH neutrinos in the latter section 4.5. Since the above relation (4.11) is sufficient for latter calculations of ∆G ′ ≷ , we do not calculate an explicit form of ∆G ′ R here. We note that, since the retarded (advanced) propagators do not have the memory effect, its deviation is essentially determined by the change of the local temperature. Also note that the enhancement factor is proportional to 1/(Ω i − Ω j ) as the Green functions in the thermal equilibrium since there is no chance to mix G R and G A .

Diagonal Wightman out of equilibrium
The deviation of the flavor diagonal Wightman function ∆G d > < (x 0 , y 0 ) can be calculated by taking a variation of (3.34): There are three terms. The first two terms are interpreted as the change of the spectrum in the expanding universe contained in G R/A . On the other hand, the third term reflects the memory effect. The third term is explicitly written 10 as This shows that the Wightman function is sensitive to the change of the background before x 0 and y 0 unlike the retarded or advanced Green functions. Writing the self-energy in terms of the center of mass coordinate X uv = (u+v)/2 and the relative coordinate s uv = u − v, its deviation from the thermal self-energy at time t = x 0 is written as (4.14) Note that |s uv | 1/Γ ℓφ due to the rapid damping of SM leptons and Higgs propagators. In the second equality the KMS relation for the thermal selfenergy (3.38) is used. As explained in eq.(4.4), the self-energy function out of equilibrium can be approximated by the equilibrium self-energy Π (eq) of (3.38) at the local temperature. Note that the distribution function f (q 0 ) = 1/(e q0/T + 1) is time-dependent through the time-dependence of the temperature T = T (X).
The calculation of the deviation of the diagonal Wightman function ∆G d ≷ is performed in Appendix I. For x 0 > y 0 , it is given by Each term of (4.15) is classified into three types of terms. The first term of ∆G d ≷ in the square bracket reflects the change of the spectrum in the propagators G R and related by the KMS relation (3.39). It reflects a change of the local temperature during the period x 0 and y 0 .
The term proportional to (X xy − t) comes from a difference between the distribution function f q (t) at the reference time t and f q (X xy ) = f q (t) + (X xy − t)d t f q at time X xy . The time-dependence of f q comes from both of the local temperature and the physical frequency ω q as shown in the definition of the derivative operator d t . The term with s xy is similar. If x 0 = y 0 , the distribution function at X xy is affected by the information at the past.
The most important part is the term proportional to 1/Γ i , which reflects the memory effect of the Wightman function. Since the Wightman function is written as a convolution G d ≷ (X xy ; s xy ) = −(G R * Π ≷ * G A )(X xy ; s xy ), they depend on the information in the past at X uv where X xy − X uv ∼ 1/Γ i (see (4.13)). In the expanding universe, the temperature is higher in the past and the number density of leptons and Higgs are larger than the present density. Accordingly the number density of the RH neutrinos is also larger by an amount of Hence the term with 1/Γ i is directly related to the memory effect of G d ≷ . In applying ∆G ≷ to the evolution equation of the lepton asymmetry, it always appears as a product with the propagators of the SM particles (leptons and Higgs) as in eq. (2.44). Since these propagators damp quickly with the decay widths Γ ℓ,φ , we can drop all the terms in (4.15) except the term containing 1/Γ i . Furthermore, during the period 1/Γ ℓφ , RH neutrinos are almost stable: Γ i ≪ Γ ℓφ . Hence we can replace the frequency Ω i by its real part ω i .
Let us write this simplified form of ∆G as ∆G: The definition of Z i ǫ is given in (3.22).
As a final remark in this section, we mention that the above simplified form is directly obtained from the classical Boltzmann equation as follows. The Boltzmann equation for the RH neutrino distribution function is given by All external momenta are on-shell. Leptons and Higgs are assumed to be in the thermal equilibrium.
is the square of the tree-level decay amplitude of a RH neutrino into a lepton and a Higgs. The spin in the initial state is averaged and the isospin sum in the final state is performed. By using the relation ( Here, we have used the definition of the decay width (3.19) with (D.9). 11 The solution of (4.20) is given by and (

becomes ∆G
(4.22) In this expression, we have assumed that the reference time t is very close to X xy , and the conditions |X xy − t|, |s xy | 1/Γ ℓφ are satisfied. Such conditions appear when we use the Wightman functions in evaluating the evolution equation of the lepton number. We also took the leading order terms with respect to Γ/Γ ℓφ ∼ Γ/T . (4.22) is of order (H/Γ). 12 We have also identified Ω i ≃ Ω j in e −iΩǫsxy since the mass difference ∆M and the widths Γ i are much smaller than the typical scale of 1/|s xy | = Γ ℓφ .
Here is an important comment. As discussed in (3.6), the off-diagonal components of the Wightman function in the thermal equilibrium (3.43) is enhanced by a large factor 1/(Ω i − Ω j ) because of the resonant oscillation between flavors. But in the limit x 0 → y 0 it vanishes as in (3.45). Both of these properties are related to the behavior of G ′ R/A through the KMS relation and the fact that G ′ ≷ is separated into the retarded and advanced propagators as in (3.42).
The deviation ∆G ′ ij ≷ does not satisfy either properties. First, the enhancement factor is replaced by 1/( The replacement of the enhancement factor by 1/(Ω i −Ω * j ) reflects the mixing between the retarded and advanced propagators. Such mixing is naturally generated because the off-diagonal component of the Wightman function is solved as in (3.36) to contain both types of Green functions. Since the retarded and advanced propagators have poles at q 0 = Ω ǫi and q 0 = Ω * ǫj respectively, the appearance of the term 1/(Ω ǫi − Ω * ǫj ) by q 0 integration can be naturally understood. In the equilibrium case, since the retarded and advanced propagators are decoupled by the KMS relation, such mixings of poles at q 0 = Ω ǫi and at ≷ (x 0 , y 0 ) in the evolution equation of the lepton number, the arguments x 0 , y 0 are restricted to the region s xy = x 0 − y 0 < 1/Γ ℓφ ∼ 1/T as mentioned above. During such short period, the decay of N i is neglected and we can safely replace Ω ǫi in e −iΩǫsxy by its real part ω ǫ . We write the simplified version of ∆G The second term in the square bracket with the real part of the self-energy can be dropped by imposing Π h = 0 by the mass renormalisation. If we include the effect of the temperature dependent mass, Π h is not always zero.

Summary of this section
In this section we studied the deviation of various Green functions from the thermal equilibrium. The deviation of the retarded Green function ∆G d R is mainly caused by the local change of the physical quantities. It is also true for the diagonal component of the Wightman function (4.18). It is because the diagonal component in the time-dependent background is determined by the distribution function at the local temperature.
In contrast, the off-diagonal components behave differently. The off-diagonal component of the retarded (advanced) Green functions ∆G ′ ij R is largely enhanced by the factor 1/(Ω i − Ω j ) due to the flavor oscillation.
The behavior of the off-diagonal components of the Wightman functions ∆G ′ ≷ in (4.22) is different. First it is not written as a change of the local equilibrium Green function G ′ (eq) ≷ in (3.44): which is controlled by the KMS relation. Or in other words, the process of taking a variation ∆ and the flavor oscillation do not commute as in (4.25). We come back to this property in section 6.
The evolution equation of the lepton asymmetry is given by the KB equation (2.44). The r.h.s. is written as a functional of the Wightman functions of RH neutrinos, SM leptons and SM Higgs. Since the SM leptons and Higgs are almost in the thermal equilibrium, their distribution functions are approximated by the thermal values at the local temperature. But the RH neutrinos decay much slower, and furthermore the RH neutrinos with almost degenerate masses coherently oscillate between different flavors during their decay. Hence the Wightman functions G ij ≷ of the RH neutrinos must be treated in a full quantum mechanical way by using the KB equation, not by the classical Boltzmann equation. Once G ij ≷ are obtained, they can be inserted into the r.h.s. of (2.44) to obtain the Boltzmann equation for the lepton asymmetry.
Here we summarize the basic ingredients of the evolution equation for the lepton asymmetry. The evolution equation is given by where P L π> < P R = π> < as defined in (2.45). After Fourier transformation of the r.h.s., the frequencies q 0 , p 0 , k 0 of the Green functions, G ≷ (q 0 ) and S> < (p 0 ), ∆> < (k 0 ) in π> < , satisfy the relation q 0 = p 0 + k 0 . Furthermore, in the thermal equilibrium, the Wightman functions are related to the retarded (advanced) propagators through the KMS relation (2.9), (3.38) and (3.39). Then, by using the relation two terms in the square bracket cancel each other. Hence there is no generation of lepton asymmetry in the thermal equilibrium.

Lepton asymmetry out of equilibrium
In the expanding universe, there are three sources for changing the lepton asymmetry, and the r.h.s. of (5.1) can be classified into the following three terms: Here we rewrite the sum over i, j into the flavor diagonal part K = d and the offdiagonal part K = ′ . Namely K = d corresponds to a summation of i = j = 1 and i = j = 2 while K = ′ corresponds to a summation of i = 1, j = 2 and i = 2, j = 1.
The first term C K ∆f comes from the deviation of the Wightman functions of the RH neutrinos (i.e., the distortion of the distribution function ∆f ) from the thermal value and is given by This generates the lepton asymmetry in the expanding universe. The second term comes from the deviation of π ≷ : which is caused by the deviation of the distribution functions of the SM leptons and Higgs. C K W is written as This gives washout effect of the lepton asymmetry. The third term comes from the back reaction of the generated lepton asymmetry to G ′ ≷ , namely to the distribution function of the RH neutrinos. It is written as Here ∆ µ G is defined as the back reaction of the generated chemical potential of the lepton and Higgs to the RH Wightman function.

Effect of ∆G ≷ on the lepton asymmetry: C ∆f
The deviation of the Wightman function from the equilibrium value generates the lepton asymmetry out of equilibrium. First let us look at the contribution of the flavor diagonal (K = d) part of C K ∆f . Inserting (4.18) 13 and (2.46) into (5.5), we have Here we took all the ǫ's, ǫ in (4.18) and ǫ ℓ , ǫ φ in (2.46), the same ǫ = ǫ ℓ = ǫ φ because the temperature considered is not so high that a process like φ → ℓ + N does not occur. Hence the flavor diagonal component does not generate the asymmetry. In the last equality, we used the relation for the thermal distribution function. Next we calculate the off-diagonal term C ′ ∆f with K = ′ . Inserting (4.24) 13 We note again that (Xxy − t) and sxy of the arguments of G ≷ (x 0 , y 0 ) are smaller than 1/Γ ℓφ due to π > < (τ, t) ∼ e −(t−τ )Γ ℓφ /2 . Hence the use of ∆G is justified. into (5.5), we have Here, using the definition of π> < in (2.46), we have defined π ρ = i(π > − π < ) = (π R −π A ), π h = (π R +π A )/2 and their Fourier transform in the time direction, to separate the self-energies Π ′ (eq) ρ/h in (4.24) into the Yukawa coupling (h † h) ′ and the equilibrium values of π ρ/h (see (D.8) and (D.14)). If we use the vacuum values 14 for the self-energy calculated in Appendix D, i.e., π ρ (ǫω q ) = −g w iǫ/ q ǫ /(16π) and π h (ǫω q ) = 0, the second term in the square bracket is dropped and (5.10) is simplified as The factor δ|M| 2 can be interpreted as the CP -asymmetric part of the decay amplitudes, which gives the CP -asymmetry of the decay rates Γ Ni→ℓφ −Γ Ni→ℓφ . The term (5.11) produces the lepton asymmetry through the CP -asymmetric decay of the RH neutrinos that are out of the thermal equilibrium. The distortion of the distribution function is given in (4.17). An important point in (5.11) is that the enhancement factor of the CP -asymmetry is given by , and the regulator R ij relevant to the CP -asymmetric decay of the RH neutrinos is given, not by

Washout effect on the lepton asymmetry: C W
The term C K W washes out the generated lepton asymmetry. In order to calculate ∆π, we first perform the Fourier transform of π ≷ (τ, t; q) defined in (2.46): where ∆D This gives a washout effect on the generated lepton asymmetry and it is physically interpreted as the inverse decay of the RH neutrinos. Next let us see the flavor off-diagonal component, K = ′ . Because of the property (3.45), it vanishes in the leading order approximation: Hence only the diagonal component plays a role of washing out the generated lepton asymmetry.

Backreaction of the generated lepton asymmetry: C BR
Finally let us see the back reaction of the generated lepton number asymmetry (i.e., the nonzero chemical potential of the SM leptons) to the Wightman functions of the RH neutrinos.
By using (D.6) and the flavor symmetry S αβ = δ αβ S, the deviation of the self-energy in the presence of the chemical potential is written as ∆π> < is the CP -conjugate of ∆π> < and obtained by changing the sign of the chemical potential of the SM leptons and the Higgs. ∆ µ G d Similarly the off-diagonal contribution becomes Details of the calculations are given in Appendix L. In the above calculations, we took the weak coupling limit discussed in Appendix D. This term represents the effect of back reaction of the generated lepton asymmetry on the Wightman functions of the RH neutrinos. Such a term appears because we first solved the propagators of the RH neutrinos in the background of the SM leptons and the Higgs. The relative sign of the back reaction to the washout effect C d W in (5.15) is opposite so that the back reaction tends to reduce the washout of the generation of lepton asymmetry. If we solve the KB equations for the lepton asymmetry and the Wightman functions of the RH neutrinos simultaneously, the generated lepton asymmetry (namely the effect of the chemical potential) makes the RH neutrinos further away from the equilibrium. It is the reason why the back reaction reduces the washout.

CP -violating parameter
The CP -violating parameter can be read off from (5.11). δ|M| 2 of (5.12) gives the CP -asymmetry of the decay rates Γ Ni→ℓφ − Γ Ni→ℓφ . Since the tree decay amplitude is given by |M| 2 tree = g w (h † h) ii (q · p), the CP -violating parameter ε i is given by Hence the regulator discussed in the introduction is given by The result is consistent with the result obtained in [56]. In the paper [56], the CP -violating parameter is obtained indirectly from the generated lepton asymmetry in a static background with an out-of-equilibrium initial condition. In our calculation, we directly obtained the same result in the expanding universe. It shows that the result obtained by Garny et al. is universal and can be applied to the thermal resonant leptogenesis.

Summary of this section
By using ∆G ′ ij ≷ calculated in the previous section 4 in the r.h.s. of (5.1), we obtained the evolution equation (5.3) with three terms. C ′ ∆f generates the lepton asymmetry and corresponds to the CP -asymmetric decay of the RH neutrinos. C W gives the washout effects on the generated lepton numbers. C BR is the effect of the back reactions of the generated lepton asymmetry on the distribution functions of the RH neutrinos. From C ′ ∆f , we extracted the CP -asymmetric parameter ε i given in (5.20). The enhancement factor due to the degenerate masses is regularized with an regulator R ij = M i Γ i + M j Γ j , which reflects the enhancement factor of ∆G ′ ≷ .

Physical interpretation of the regulators
In this section, we give a physical interpretation of the appearance of the regu-  These equations mean that the information of the distribution function of the RH neutrinos in the Wightman functions G ij ≷ are encoded in the self-energies Π kl ≷ in the past, and transferred from the past to the present t = x 0 , y 0 . The self-energies Π kl ≷ encode the information of the distributions functions of the SM leptons and the Higgs in the past (see Fig.3). In the flavor diagonal case of (3.34), all flavors of the RH neutrinos are the same in the leading order approximation. On the other hand, in the flavor off-diagonal case of (3.35), the flavor oscillation plays an important role.
Here we note that, as shown in (3.28) and (3.29), G ′ ij R/A is a coherent sum of two terms, each of which corresponds to a propagation of the i-th (or j-th) flavor RH neutrino. We divide it as follows: A precise definition is given in (M.80).
6.1 On-shell and off-shell separation of G ′ (eq) ≷ Now let's investigate G ′ ij > < . By looking at the first term of (3.35), it contains G djj A which describes the propagation of the j-th RH neutrino. The propagator G ′ ij R in the first term contains both of the propagations of i-th and j-th flavor neutrinos. If the j-th neutrino propagates in G ′ ij R , only a single (j-th) neutrino propagates from the past, when the decay/inverse-decay represented by Π ≷ takes place, to the present at t = x 0 , y 0 . We call this type of contributions the "on-shell" contributions. 15 These contributions are all taken into account in the classical Boltzmann equation.
On the contrary, if the i-th neutrino propagates in G ′ ij R , two different flavors propagate from the past to the present. This type of contributions are essentially "off-shell". In the classical Boltzmann equation, we first calculate the S-matrix elements of various processes and the external lines are taken to be on-shell. Hence this type of "off-shell" contributions are not taken into account by ordinary methods. 16 If we neglect the off-shell terms and take only the on-shell terms, G Note that the sum of the on-shell contributions do not vanish even at x 0 = y 0 and f i ≃ f j : It is different from the property of the full contributions given in (3.43).
6.2 On-shell and off-shell separation of ∆G ′ ≷ We next investigate ∆G ′ ≷ . We show that neglecting the off-shell contribution in ∆G ′ ≷ , we get an enhancement factor for the CP -violating parameter with a regulator In Appendix M.7, we separate ∆G ′ ≷ into on-shell and off-shell contributions: 17 The full ∆G ′ ≷ is given in (J.34). The on-shell contribution is given by (for This on-shell contribution has two important properties. First, it satisfies on-shell is given in (6.3). The on-shell contribution of (M.107) is simply obtained by replacing f (eq) by its variation ∆f in (6.3). This replacing means that the process of the flavor oscillations and the process of taking a variation from the thermal values are commutative if we neglect the off-shell contributions. For full quantum calculations, (3.44) cannot be obtained by such a replacement from (3.43). This is because the flavor oscillations and the deviation from the thermal values are coherently mixed and these processes are not commutable. Namely, dropping the off-shell contributions corresponds to neglecting the interference between the flavor oscillations and the deviation of the distribution functions from the thermal equilibrium.
Second, compared with the full result (J.34), the enhancement factor 1/(Ω i − Ω * j ) is replaced by 1/(Ω i − Ω j ). It is related to the above non-commutativity of taking ∆ and flavor oscillation effects.
By inserting the on-shell formula (6.3) and (M.107) into (5.5), and supposing π ρ (ǫω q ) = −g w iǫ/ q ǫ /(16π), π h (ǫω q ) = 0, we have an on-shell approximation C ′ ∆f on-shell of C ′ ∆f : initial condition in the flat space-time, they found two different behaviors of the generated lepton number. One is the ordinary term common in the conventional Boltzmann equation. The other term is specific to the quantum treatment by the quantum KB approach. The latter oscillates in time and reduces the eventual lepton number. "Off-shell" contribution here corresponds to the latter effect. However, note that in the present case the CP -violating parameter, and hence the resulting lepton number does not oscillate. the oscillatory behavior is averaged out because the deviation from the equilibrium is caused by the expansion of the universe, and its expansion rate H is much smaller than the oscillation scale ∆M ≃ Γ. This averaging also occurs in the analysis by [71] in the strong washout regime.
Hence the regulator |M i Γ i +M j Γ j | is replace by |M i Γ i −M j Γ j | if we take only the on-shell terms. It is not valid in general, especially in the resonant leptogenesis. If the masses are hierarchical, it becomes identical with the correct value in (5.11).

Summary of this section
As emphasized above, if we neglect the off-shell contributions that are not included in the ordinary Boltzmann type analysis, we get a result (6.6) which is different from the correct one given in (4.22). The only difference is the enhancement factor, and if the mass difference is much larger than the width they coincide. But the difference is enlarged when the masses are almost degenerate. This reflects the fact that the flavor oscillation becomes important only for degenerate masses. Another important point is that the property of the noncommutativity Based on this observation, we give another derivation of the properties of ∆G ′ ≷ in Appendix K by directly solving the KB equations. If we assume the vanishing condition (K.64) of G ′ ≷ which is equivalent to (3.45), we show that the enhancement factor with a regulator M i Γ i + M j Γ j appears as in (K.65). On the other hand, if we erroneously assume that it does not vanish, it leads to a much larger enhancement factor.

Summary
We investigate the Kadanoff-Baym equations of the thermal resonant leptogenesis in the expanding universe. The lepton asymmetry is generated in the CP -asymmetric decay of the RH neutrinos which are deviated from the thermal equilibrium. If the RH neutrinos have almost degenerate masses, they coherently oscillate during their decay into the SM leptons and the Higgs. In such a situation, the classical Boltzmann equation is not valid because the decays and the inverse-decays of the RH neutrinos cannot be separated into different processes, and the full quantum mechanical approach is necessary. A systematic approach is given by solving the KB equations.
Kadanoff-Baym approach to the resonant leptogenesis was performed in the previous analysis [56]. In the paper, the authors studied the coherent oscillation of the RH neutrinos in a time-independent background with a non-equilibrium initial condition. In the process of approaching the equilibrium, the RH neutrinos coherently oscillate and decay into the SM particles. From the generated amount of the lepton number, they extracted the CP -violating parameter ε i of the i-th RH neutrino (1.1) and obtained the regulator R ij = |M i Γ i + M j Γ j |. Since the resonant enhancement is the key to the resonant leptogenesis, it is essential to obtain the correct form of the regulator.
In the present paper, by extending the analysis in the static background [56] to the thermal resonant leptogenesis in the expanding universe, we obtain an analytical expression of the evolution equation of the lepton number asymmetry. The CP -violating parameter is obtained as in (5.20). The regulator we obtained is consistent with the result [56].
The difference between the regulator R ij = |M i Γ i +M j Γ j | and R ij = |M i Γ i − M j Γ j | obtained by [57] comes from different forms of the enhancement factors of flavor off-diagonal components of the RH neutrino propagators. Since the RH neutrinos have almost degenerate masses, they coherently oscillate very much. We show that the resonant oscillations between different flavors have two different types, one proportional to 1/(Ω i − Ω * j ) and the other proportional to 1/(Ω i − Ω j ). Here Ω i = ω iq − iΓ iq /2 is a position of the pole of the ith RH neutrino. In the thermal equilibrium, the resonant oscillations in the flavor off-diagonal Green functions have the type 1/(Ω i − Ω j ) ( or its complex conjugate) as shown in (3.28) or in (3.44). Since 1/(Ω i − Ω j ) is rewritten by (3.33), the enhancement of flavor oscillation corresponds to the regulator R ij = |M i Γ i − M j Γ j |. However, the deviation of the off-diagonal components of out of equilibrium has a different enhancement factor 1/(Ω i − Ω * j ) as shown in (4.22), which corresponds to the regulator R ij = |M i Γ i +M j Γ j |. Physical interpretation of the change of the regulator is given in Section 6. The off-shell contributions to the off-diagonal component to the Wightman functions are essential which can be incorporated only by the KB equations. The property (3.45) is also important for this change of the regulator. As we show in Appendix K, if we erroneously assume that the off-diagonal component of the Wightman function is non-vanishing and its deviation is given by the change of the local temperature, it leads to much more enhanced oscillation similar to the regulator R ij = |M i Γ i − M j Γ j |.
In the present paper, we focus on the formalism of the thermal resonant leptogenesis and derivations of the CP -violating parameter in the expanding universe. Phenomenological investigations of the amount of the generated lepton asymmetry and the lower bound of the leptogenesis scale consistent with the neutrino oscillation data are left for future investigations.

A CTP formalism
In this appendix, we briefly review the closed time path (CTP) formalism.
In equilibrium field theories, we implicitly assume that the initial and final states asymptotically approach the ground state of the free Hamiltonian. But this does not hold in general, especially in time-dependent backgrounds such as the evolving universe. The final state is generally different from the initial state. The CTP formalism, or the Schwinger-Keldish formalism, is the general formalism to calculate physical quantities for time-dependent wave functions. Suppose that a system is described by a HamiltonianĤ 0 +Ĥ 1 , whereĤ 0 andĤ 1 are free and interaction Hamiltonians, and that the system is in the initial state |ψ i at time t = t i . In the interaction picture, the expectation value of an observableÔ at time t is given by Here the operator in the interaction pictureÔ I (t) is related to the operator in the Heisenberg picture aŝ In equilibrium cases, the final state at time t = t f is assumed to be proportional to the initial state U I (t f , t i )|ψ i = e iθ |ψ i where θ(t f , t i ) is a c-number phase. Then we can factorize O(t) of (A.1) as Similarly, an expectation value of the time-ordering product of two operatorŝ O I 1 (t 1 ) andÔ I 2 (t 2 ) is given by This formula gives an ordinary perturbative expansion of correlation functions in equilibrium field theories. Namely, if we take t i → −∞ and t f → ∞, the interaction verticesĤ I (t) are inserted in −∞ < t < ∞.
In non-equilibrium cases where the final state is no longer proportional to the initial state, the factorization property does not hold and we have In perturbative expansions, the interaction vertices are inserted not only on the path C + from t i to t f , but also on the backward path C − from t f to t i . Figure 4 shows the closed time path (CTP), C = C + + C − . In this formalism, the final state is not specified at all and we can calculate time-dependence of various quantities as in (A.5). The time-ordering T C is defined on the CTP as a path-ordering along C = C + + C − .

B Evolution equations of various propagators
We define various propagators and give a brief derivation of their evolution equations. In the following we consider a real (Majorana) fermion fieldψ and write its conjugate byψ =ψ t C where C = iγ 2 γ 0 . In the CTP formalism, a generating function of time-ordered products of operators is given by For notational simplicity, we use the following abbreviation unless the explicit dependence of the measure on x is necessary. The stationary condition (δ/δΨ(x))Γ[Ψ] = 0 at J = 0 gives the equation of motion of Ψ. By taking the second derivative of the effective action Γ with respect to Ψ, we obtain the Schwinger-Dyson (SD) equation Π is the self-energy and only 1PI diagrams contribute to it. The connected Green function G on C is defined by where Θ C (x 0 − y 0 ) is the step function on C and are the Wightman Green functions.
is an inverse of the free propagator and Π(x, y) is the self-energy of the fermion field ψ.
The statistical propagator G F (x, y) and the spectral function G ρ (x, y) are defined by G F contains information of the distribution function of the specified state while G ρ depends only on the spectrum of the system. In this sense, G F is dynamical while G ρ is kinematical. Especially, G ρ (x, y) becomes proportional to the spatial delta-function δ (3) (x−y) in the equal-time limit. We further define the retarded and advanced Green functions by They are related to G ρ as In terms of G F and G ρ , the Green function (B.6) can be written as where the sign-function on C is defined by By convoluting (B.5) with the full propagator G, we have (B.14) Here δ g C (x − y) is the delta-function on C with the space-time metric g, and satisfies C d 4 z g δ g C (x − y) = 1. By denoting x on C ± as x ± respectively, the delta-function on C can be expressed by a 2 × 2 matrix: where a, b takes + or −. The minus sign on C − comes from the backward integral of the time variable and corresponds to the anti-time-ordering of the Green function G in (B.6). δ g (x − y) is an ordinary delta-function for (x − y). The 2-point function G(x, y) of (B.6) with x, y ∈ C can be similarly decomposed (depending on whether x, y are on C + or C − ) into a 2 × 2 matrix form as where T, T denote time and anti-time orderings respectively, and Θ(x 0 − y 0 ) is the ordinary step-function. By using (B.12) and (B.11), we have We also decompose the self-energy Π(x, y) as the matrix form of the self-energy is obtained as Using these matrix forms of G ab and Π ab , the equation (B.14) becomes The matrix c cd between Π and G comes from the backward integration of the time variable in the original integral in (B.14). By multiplying (U t ) −1 on the left and U −1 on the right, using (B.18) and (B.22) and noting we obtain the following set of the evolution equations: From the equations (B.26), we obtain the evolution equation for the spectral density G ρ = G R − G A : The Wightman Green function (B.28)

C 2PI formalism
In this appendix, we give a brief review of a systematic approach to evaluate the self-energy based on the 2PI formalism (see, e.g., [92] for more details). The generating functional Z[J, R] in the presence of sources J(x) and R(x, y) 18 is given by .
By taking a variation with respect to the source fields J(x) and R(y, x), we have Here ζ, η represent Spinor indices. Ψ is defined as Ψ(x) ≡ Ψ t (x)C and the connected Green function G is given by By taking the Legendre transform of W [J, R] with respect to the sources J, R, we obtain the effective action in the presence of source fields TrGR .

(C.4)
Tr in the last term represents a trace in the Spinor indices and an integration over the closed time path C. Now we decompose the effective action into The first term is the classical action. The second and the third term are '1-loop' type contributions to the effective action. The meaning of the decomposition can be understood by taking a functional derivative 19 with respect to G: Compared with the SD equation (B.5), the last term can be identified as the self-energy Π: In this way, the proper self-energy Π is obtained by differentiating Γ 2 with respect to the full propagator G. Since the proper self-energy Π is calculated as a sum of contributions from 1PI diagrams, Γ 2 becomes a sum of contributions from 2PI diagrams with respect to the full propagator. In another word, the proper self-energy can be systematically obtained by taking a functional derivative of 2PI diagrams (in which all internal lines are full propagators) with respect to the full propagator. The SD equation (C.6) can be interpreted as a self-consistent equation for the full propagators. By rewriting this equation on the forward time line C + , For Dirac or Weyl fermions, we introduce additional source terms The self-energy is similarly obtained as a functional of the full propagators: In this appendix, using the 2PI formalism, we give an expression of the selfenergy function for the RH neutrino Π(x, y) in terms of the lepton and Higgs propagators. The simplest and the most important contribution to the 2PI effective action Γ 2 in the model (2.2) is given by the 2-loop diagram of Figure 2. The second simplest 2PI diagram is given by Figure 5. Note that each internal line represents a full propagator of the SM lepton, Higgs and the RH neutrino. If the RH neutrinos have almost degenerate mass, we need to use the resummed propagators for the RH neutrinos. Once resummed, we can use an ordinary perturbative expansion with respect to the Yukawa coupling h iα . Hence it will not be a bad approximation to use the simplest 2PI diagram to evaluate the self-energy. In terms of the full propagators, the contribution from the diagram Figure  2 becomes Here G, S, ∆ are full propagators of the RH neutrino, the SM lepton doublet and the Higgs doublet respectively. (i, j), (α, β), (a, b, a ′ , b ′ ) represent the flavor indices of the RH neutrino, the flavor indices of the leptons and the SU (2) L indices of the SM doublets respectively. By using the formula (C.8), the self-energy of the SM lepton doublet is given by taking a functional derivative of Γ 2 with respect to the lepton propagator S: Here we have used the Majorana property G ij (x, y) = CG tji (y, x)C −1 of the RH neutrinos. In the second equality, we have used the fact that the lepton and the Higgs propagators are SU (2) L symmetric and proportional to δ ab , S ab = Sδ ab , ∆ ab = ∆δ ab , in the early universe where the SU (2) L symmetry is restored. This is indeed the case in the era of the lepton asymmetry generation through the decay of the RH neutrino. Similarly the self-energy of the RH neutrino is obtained by taking a functional derivative of Γ 2 with respect to G: where P = γ 0 . In the first equality, we have used In the following, we derive the self-energy of the RH neutrino Π (eq) under an assumption that the lepton and the Higgs are in the thermal equilibrium. The approximation is justified in the leading order calculation since the SM leptons and the Higgs particles interact faster than the Hubble expansion rate in the era of the leptogenesis. See (2.43). Hence the deviation from the equilibrium can be neglected in the calculation of Π. In the equilibrium, the lepton and the Higgs propagators become CP -symmetric and satisfy S(x, y) = S(x, y) , ∆(x, y) = ∆(x, y) . (D.7) By using the quasi-particle approximation for the propagators (2.37) and (2.38), the Fourier transform of the self-energy Π ρ = i(Π > − Π < ) = Π R − Π A of the RH neutrino becomes with the definition of D ≷(p,k) in (2.47), and satisfy the relation D and is followed by the relation of In the calculation we have used the integral (in the limit t int → −∞) The contribution from the boundary at τ = −∞ vanishes because of the damping factor ∼ e −Γ ℓφ (t−τ )/2 . In the weak coupling limit of the SM gauge couplings, Γ ℓφ becomes much less than the typical energy transfer (q 0 − ǫ ℓ ω p − ǫ φ ω k ) ∼ T where T is the temperature at which the leptogenesis occurs. In such a limit, the above integral becomes proportional to δ(q 0 −ǫ ℓ ω p −ǫ φ ω k ), and the exact energy conservation is satisfied instead of the Lorentz type in (D.9). Furthermore, in order to simplify the form of the self-energy (D.9), we neglect the medium effects (e.g., the Pauli exclusion of the SM lepton and the induced emission of the Higgs) encoded in D ρ in (D.9) and drop the distribution function f .
Adopting these two simplifications of the weak coupling limit and neglecting the medium effects, the self-energy (D.9) reduces to the vacuum one: Since the main purpose of the present paper is to obtain the effect of quantum oscillations of almost degenerate RH neutrinos, we use this simplified form of the self-energy. The full treatment is investigated by using the integral form (D.9) of the self-energy instead of (D.13).

(D.15)
It satisfies the relation Note that π (eq) ρ (q) is pure imaginary while π (eq) h (q) is real. The real part π h (q) contains a diverging integral which is subtracted by the mass renormalisation. In the body of the paper, we have implicitly assumed that the self-energy π (eq) h (q) is already regularized. The imaginary part π (eq) ρ gives a decay width of the RH neutrino.

E Kramers-Moyal product
The convolution * is defined on bi-local functions f (x, y) and g(x, y): The Wigner representation of bi-local functions is also defined as the Fourier transform of the relative coordinate as Then it is straightforward to show that ( f * g)(X; p) = dp 1 2π where ∂ f X is a Xderivative on the function f . The non-commutative product is called the Kramers-Moyal product. In the leading order approximation of the derivative expansion, the commutator of the * -product is reduced to the Poisson bracket:

F Useful identities
The Green functions are written in terms of convulsions, and by taking a variation, various functions are inserted in the integral. Thus we ofter encounter the following types of integrals: Here G i or Π are assumed to be functions of the relative coordinate only. X uv = (u + v)/2 etc. are the center of mass coordinates. In order to evaluate these integrals, we consider the following identities. For the integral (F.1), By acting (i∂ Q − t) on both sides and setting Q = 0, we obtain ∂ (i) are derivatives acting on G i . Taking higher derivatives with respect Q, we can obtain other relations. For the next type integral (F.2), we start from the following identity: 2 )sxy−i(Q1+Q2+Q)Xxy .

G Calculation of ∆G d R/A
The deviation of the retarded (advanced) propagators out of equilibrium can be calculated by taking a variation of (2.11). In this section, we consider the diagonal component G d R/A . By dropping the higher order term Π ′ G ′ , the diagonal component of the equation (2.11) is written symbolically as Then expand D and G around the thermal equilibrium values (at the reference time t) as D = D (eq) t + ∆D and G = G (eq) t + ∆G. Inserting them into above we have R/A = −δ xy and dropped the higher order term ∆D∆G of deviations from the equilibrium. It can be solved as The deviation ∆D R/A comes from the change of the physical momentum through a(t) and the change of the self-energy Π R/A . Thermal corrections to the mass M are included in Π.
Writing the integral explicitly, ∆G d becomes For notational simplicity, we did not write the reference time t explicitly in ∆D R/A (u, v).
Note that the equilibrium self-energy Π is a function of the relative coordinate u − v with the temperature at time t. Since the self-energy Π is obtained by loop integrals of the SM particles, it decreases rapidly as e −Γ ℓφ (u−v)/2 . Hence we can adopt derivative expansions of the selfenergy Π(u, v) around the thermal value Π The first term is the change of the physical momentum. The second term is the change of the background SM plasma, i.e. the change of the distribution functions and the spectrum of the SM leptons and the Higgs. Let us write the Fourier transform of the coefficient in the square bracket of (G.6) with respect to (u − v) as The retarded (advanced ) Green functions in the (local) equilibrium are given by We also define G . Physical quantities are determined by the local temperature at time X. Then by taking derivatives with respect to X and q 0 , we have Let us now calculate ∆G d R/A (x 0 , y 0 )) for x 0 > y 0 . Inserting (G.6) into (G.5) and using the identities (F.4) and (F.6), we get (G.14) By using (G.13), U 1 becomes In the first equality, the triple pole contributions are canceled out each other.
In the second line, only the lowest order terms in Yukawa couplings are taken.
In the third line, we have performed the q 0 integration and neglected higher order terms in Γ/M ≃ h † h. 20 U 2 can be calculated much easier because the time derivative in U 2 is commutable with the q 0 integral. Then we get Summing up (G.15) and (G.16), we get the time representation of the devi-ation of the retarded(advanced) propagator: The first term comes from U 2 , It is written as the change of the mass and width in Ω ǫ = ǫω ǫq − iΓ q /2 and the physical quantities in Z ǫ . The second term comes from U 1 . It represents a change of the spinor structure in the propagator during the propagation in the expanding universe.

H Calculation of ∆G d ρ
The deviation of the spectral density from the equilibrium value is obtained by taking a variation of the relation G ρ = −G R * Π ρ * G A : In this appendix we explicitly evaluate these terms since it is similar to and instructive for more complicated calculations of ∆G d ≷ in Appendix I. By using (G.14) and (F.7), the first term in the r.h.s. of ∆G d ρ is given by (v, y 0 ) ≡ T 11 + T 12 + T 13 + T 14 , (H.2) (T 12 + T 13 + T 14 ) came from U 2 while T 11 came from U 1 . By using (G.14) and (F.8), the second term is Finally, by using (4.14) and (F.6), the third term becomes Let us look at the above terms separately. First, the terms proportional to (X xy − t) become T a : Next, the terms proportional to s xy can be rewritten as T b + T 44 : (H.6) T b is a total derivative and vanishes, but we keep it for later convenience in calculating ∆G d ≷ in the next section. The other terms and T 44 are combined to become where we used the relation as expected. Here we used (G.14). In the second equality, T b is dropped since it is a total derivative. In the calculation of ∆G d ≷ , T b is modified to contain a function f (q 0 ) outside the derivative ∂ q0 , and contributes to the final result.
As explained in (4.4), the self-energy function can be safely replaced by its equilibrium value at the local temperature: The distribution function f (q 0 ) = 1/(e q0/T (t) + 1) is time-dependent through the temperature T (X and we have V 1 is obtained by inserting the distribution function in the integrand of (H.8) for ∆G ρ . V 2 comes from T b -type term of (H.6) in the previous appendix. In the calculation of ∆G ρ , this term vanishes because it is a total divergence. In the present case, since the distribution function depends on q 0 , the derivative can act on it and the term remains. V 3 and V 4 come from time-dependence of the distribution function f (q 0 ). Let us consider the region x 0 > y 0 , and perform the q 0 integration. First, V 4 and T a -type contribution (H.5) in V 1 turn out to represent the change of local temperature as (I.4) T c -type contribution (H.7) in V 1 is of higher order with respect to Γ i /M and can be neglected. 21 Therfore, we get Let's move on to V 3 . By using the explicit expression (G.13) of ∂ q0 G d(eq) R/A , we get It is nothing but the one obtained by inserting the distribution function into U 2 in the retarded propagator. In (G.15), we neglected higher order contribution in Γ/M . In the present case, although q 0 derivative, coming from the double pole integral, can act on the distribution function f (q), the contribution becomes a higher order contribution because of ∂q 0 f (q 0 ) ∼ f (q 0 )/T ∼ f (q 0 )/M , and can be neglected again. where Note that the number of the factor G S /M in V 33 is less than those in the other two terms. Hence V 33 turns out to be a higher order with respect to Γ/M and can be neglected. The q 0 integration of V 31 and V 32 can be similarly performed as in U 1 , and we get In the first line, half of the term containing 2/Γ 2 q comes from V 32 , and the other half of it arises from the q 0 derivative of G AS /M . The term s xy /Γ q comes from ∂ q0 e −iq0sxy in V 31 . Note that time-derivatives in (I.6) act on the distribution function only through the time dependence of temperature. In the second line, we have used the definition of the width.
Noticing that V 2 is obtained from V 3 by changing its sign and exchanging ∂ q0 and ∂ t . Then we can immediately get the time representation for V 2 . Similarly to (G.13), we get An exchange of ∂ q0 and ∂ t leads to an exchange of ξ and ζ. From the definitions of ξ in (G.12) and ζ in (G.7), the exchange results in in (I.7) in the leading order. Hence we get the expression of V 2 as (I.10) Putting (I.7) and (I.10) together, is obtained. Finally, from (I.5)+(I.11), we get the final result: In the leading order approximation, off-diagonal components of the Wightman function are given by eq. (3.36). Its deviation from the thermal value is obtained by taking a variation: We will calculate these 9 terms in the following sections. Among all these 9 terms, it will turn out that leading order contributions to ∆G ′ ≷ that eventually remain in the limit come from three terms containing variations of the Wightman functions ∆G d ≷ or ∆Π ′ ≷ , namely the 3rd term, the 4th term and the 8th term in (J.1). These terms represent variation of the distribution function ∆f , and are shown to be of order O(H/Γ). After summation of various terms, the other terms become O(H/Γ) × O(Γ/Γ ℓφ ) etc. and negligible in the above limit compared to the leading terms. We first give order estimations of the following quantities. First, since the distribution function f (q 0 ) contains a factor e −q0/T (t) , the derivatives are We have usedṪ ∝ HT and that the typical frequency in the analysis is q 0 ∼ T . On the other hand, the derivatives of Green functions (G.10) and (G.11) are Green functions contain a factor 1/(q 0 −Ω j ) or 1/(q 0 −Ω * j ) where Ω = ω i −iΓ i /2. In q 0 integrations, q 0 is replaced by positions of poles coming from other similar factors, and these factors become 1/(Ω i − Ω j ) for i = j or 1/(Ω i − Ω * j ). Both of them are of order 1/∆M ∼ 1/Γ. Acting ∂ t , we have (J.5) Then replacing q 0 by Ω i , the first relation of (J.4) is obtained. The second relation comes from the relation and 1(q 0 − Ω j ) → O(1/Γ). Hence, for both of the derivatives, (∂f )/f is smaller than (∂G)/G by a factor Γ/T . (∂Π)/Π is also smaller than (∂G)/G. It is because of the following reason. The self-energy Π(q 0 ) contains a factor 1/(q 0 − (ω p + ω k ∓ iΓ ℓφ /2)) where ω p and ω k are energies of the SM lepton and the Higgs propagating in the selfenergy diagram. After performing q 0 integration, q 0 is replaced by the position of poles Ω i,j in G(q 0 ). Due to (J.2), 1/(q 0 − (ω p + ω k ∓ iΓ ℓφ /2)) becomes of order 1/T ∼ 1/Γ ℓφ unlike 1/Γ for G. Therfore, acting derivatives, (∂Π)/Π becomes smaller than (∂G)/G by a factor Γ/T .

J.1 Leading contributions
We first calculate three terms that are obtained by taking variations of the Wightman-type functions. The other terms, which turn out to be subleading, are evaluated in Appendix J.3.
Since ∆G d ≷ is written as a sum of four terms of V 1,2,3,4 in (I.3), the 3rd term in (J.1) can be also written as a sum of the following four terms V ′ 1,2,3,4 respectively. By using (I.3), (G.14) and (F.8), it becomes The first term V ′ 1 is written as V ′ 2 and V ′ 3 are given by and Both of them remain in the above limit.
Finally the last term V ′ 4 is written as Similarly, by using (I.3), (G.14) and (F.7), the 4th term in (J.1) becomes The arguments are parallel to V ′ and we do not repeat the discussions. The last term containing a variation of the Wightman function comes from the 8th term in (J.1). It contains ∆Π ′ ≷ , and becomes (J.24) Only the second term of V ′′′ 2 becomes O(H/Γ) and remains in the limit. The first terms of V ′′′ 1 and V ′′′ 2 are combined to be (J. 43), and the second term of V ′′′ 1 is combined to be (J.46); these are negligibly small.

J.2 ∆G ′ ≷ in the time-representation
We now calculate the time-representation of ∆G ′ ≷ by Fourier transforming the leading contributions in the previous subsection J. 1. The other terms in (J.1) are evaluated in the next subsection and after being combined they are shown to become negligible. We consider the case x 0 > y 0 in the following. The calculations are performed in parallel to those in Appendix I, so we do not repeat detailed calculations.
We suppose that the reference time t and the arguments x 0 , y 0 in ∆G(x 0 , y 0 ) satisfy the conditions, (X xy − t) < ∼ 1/Γ ℓφ and s xy < ∼ 1/Γ ℓφ since such situations appear in the calculation of the lepton asymmetry (5.1). It is due to the fast damping of the self-energy Π ≷ . In performing the Fourier transformation, we drop higher order terms with respect to H/Γ compared to the leading terms of order O(H/Γ).
For V ′ 1 , the leading order contribution with respect to H/Γ comes from V ′ 14 . Performing q 0 integration, it becomes ). The first line is the contribution of the pole at q 0 = Ω j while the second term corresponds to the residue of the pole at q 0 = Ω i . We have dropped higher order terms, e.g., a residue at the pole of the self-energy Π ′ or the distribution function f since they are much more suppressed by a factor Γ/T . The second equality is obtained by identifying Ω i and Ω j in the leading approximation.
V ′ 2 and V ′ 3 are evaluated as (J.27) For V ′ 4 , the leading contribution comes from V ′ 41 . It becomes Similarly, we get the leading order terms from V ′′ as Finally, the leading order contributions in (J.24) comes from the second term of V ′′′ 2 . It becomes Summing up all these contributions (J.25)∼(J.33), we get the final expression (J.34)

J.3 Irrelevant contributions
In this Appendix we will see that after being combined, the other terms in (J.1) neglected in Appendices J.2 become subdominant in the limit of H ≪ Γ and (X xy − t) < ∼ 1/Γ ℓφ , s xy < ∼ 1/Γ ℓφ . The neglected terms are obtained by taking variations of the retarded or advanced-type functions. By using (G.14) and (F.7), the first term in (J.1) becomes The second term becomes (J.36) The 5th and 6th term become (J.38) The 7th and 9th term become .   (J.23) and the third line of (J.36) and (J.37), which contain ∂ t Π R,A , are combined to be Each of them is of order O(HT /Γ 2 ) but each line becomes of order O(H/Γ). Moreover, because of the cancellation between the first and second line, or the third and fourth line, it turns out to be of order O(H/T ) so that they are negligible again.
K Another derivation of ∆G ′ ≷ In this appendix we give another, quick and heuristic, derivation of ∆G ′ ≷ . The derivation use some of the results justified in the systematic derivation adopted in the paper. First we assume that the deviations of the Wightman functions from the thermal value at time t is given by the following form: The Wightman functions in the thermal equilibrium at t are given in (3.40) for the diagonal component and (3.44) for the off-diagonal component. Hence they are similarly written in terms of A and B as for the diagonal component and for the off-diagonal component.

≷
In the following we obtain the deviation of the off-diagonal component of the Wightman functions ∆A ′ > < directly by solving the KB equations using the above information. The KB equations for the off-diagonal Wightman functions are given by Setting x 0 = y 0 = t and take a difference of these two equations. Summing over the spinor indices, we have On the other hand, multiplying γ 0 and then summing over the spinor indices, we have We are now interested in the deviation from the thermal values at time t. The equations (K.54) and (K.55) are rewritten as (K.58) The first line of each equation is nothing but the l.h.s. of (K. 54 With this relation, (K.57) and (K.58) are simplified to be (K.61) (K.49) indeed satisfies (K.60). The second term of the l.h.s. of (K.61) vanishes in the leading order approximation (K.47), but using the next to leading order approximation of ∆G d , the second term becomes −iΓ q (|q| 2 /a 2 )/(M ω q ). Then eq. (K.61) is satisfied.
Then we study the off-diagonal component. Using Π(ω q , q) ∝ / q, (K.57) and (K.58) become Here we absorbed the real part of the self-energy Π h into the mass term in the l.h.s. by the mass renormalisation. For the off-diagonal component, we can expect that This comes from the relation (3.45) or equivalently (K.51). Since it vanishes in the thermal equilibrium, its variation due to the change of the local temperature is also expected to vanish in the leading order approximation. Indeed it is confirmed by using (J.43). On the contrary, since the equilibrium diagonal Wightman function survives in the same limit, its variation (or ∆Ḃ d = 0) does not vanish either. Furthermore, the second term in the l.h.s. of (K.63) is also expected to give no leading contribution like the first term ∆Ḃ The regulator M i Γ i + M j Γ j controls the enhancement of the solutions for ∆G ′ ≷ . This expression corresponds to (4.24).

K.5 ∆G
′ ≷ based on a wrong assumption G ′ ≷ = 0 Finally in this appendix, we discuss how we will obtain an erroneous answer with the regulator of the type M i Γ i − M j Γ j . Let us assume (which turns out to be wrong) that the off-diagonal component did not vanish and is given by Here Γ q = ΓM/ω q is of the same order as Γ iq = Γ i M/ω q (i = 1, 2). These are similar to the correct relations for the diagonal components, (K.60) and (K.61). The above equations (K.66) and (K.67) are based on a correct-looking assumption that the deviations of the off-diagonal Wightman functions out of equilibrium are obtained by taking a variation of the the equilibrium value with respect to the local temperature. In other words, it is assumed that there exists an "off-diagonal distribution function f ′ (eq) q " which does not vanish at s xy = x − y = 0 and its deviation from the equilibrium value satisfies the relation ∆f Under such incorrect assumptions, additional terms change the l.h.s of (K.62) and (K.63), and the regulator is modified to be This is the way we could obtain an erroneous enhancement factor.

L Effects of backreactoin
Backreactions of the generated lepton asymmetry to the RH Wightman functions are given by inserting (5.17) into (3.34) and (3.35): Here by using (D.6) ∆ µ Π> < (q) becomes . (L. 74) In the leading order approximation of small deviations from the thermal equalibrium, the square bracket in (L.72) becomes .

(L.75)
Let us write (L.69) and (L.70) as We have adopted the approximation e −iΩis ≃ e −iΩj s ≃ e −iΩs because we are especially interested in the case 1/Γ ℓφ > ∼ s > 0. Also we have dropped rapidly damping contributions ∼ e −sΓ ℓφ /2 from the pole of self-energy which are the higher order in Γ/Γ ℓφ . Plugging this form into (5.1) with the equilibrium value of π> < and performing the time integration, we get at their time representations, (3.28) and (3.29): For the out-of-equilibrium Green functions ∆G ′ R/A , the separation is a bit more involved. Let us start with the expression (4.11) of the off-diagonal retarded propagator: (M.82) By using identities in Appendix F, we can calculate these integrals. W 1 becomes The leading order terms containing waves with frequency Ω i is given by taking the residue of the Ω i in G d(eq)ii R and we have . (M.84) Here we used the identity (F.5). To extract a leading order contribution with frequency Ω j , it's convenient to rewrite W 1 into the following form: It is equal to (M.84) up to a total derivative. Taking the residue of the pole Ω j in ∆G djj R , we get the leading order contribution with frequency Ω j as .

(M.86)
Similarly, we can extract the following leading order contributions from W 2 : (M. 89) It turns out that all terms are suppressed by a factor Γ/T and there are no leading order contributions from W 3 : Summing up W i (i = 1, 2, 3), the retarded Green function ∆G ′ R is separated in the leading order approximation as . (M.91)

M.3 Useful identities
In order to perform the following calculations, we first introduce two useful identities dq 2π Note that the complex frequency Ω is introduced in the above identities unlike (F.5). In the following, we will use these identities together with (F.5).
M.4 G ′ ≷ By using the decomposition (M.80) and (M.81), we separate G ≷ ′ (eq) into "onshell" and"off-shell" terms as (M.94) In the "on-shell" terms, the same mass eigenstate i or j propagates 22 . On the other hand, the "off-shell" terms contain both of the mass eigenstates i and j simultaneously.
Summing the above two on-shell contributions, we have eq. (6.3). 22 Since Π ′ is flavor off-diagonal, differences between the flavor eigenstates and the mass eigenstates in G R/A are higher orders with respect to (h † h) ′ /(h † h) d .
M.6 Off-shell part of G ′ (eq) ≷ The first off-shell contribution comes from the i-th propagation in G ′ ij R of the first term of (3.35), and becomes (supposing x 0 > y 0 ) (M.97) The j-th propagation of the second term in (3.35) gives the off-shell contribution and becomes (for

(M.98)
From the last line of (M.94), we get the off-shell contribution by using the identity (F.5) (for x 0 > y 0 ): M.7 On-shell part of ∆G ′ ij ≷ Finally, we move on to ∆G ′ij ≷ . Taking a variation of (3.35), we get the following 7 terms: (M.100) 23 Note that the frequency of the wave function is given by Ω ǫi . It is because, for x 0 > y 0 , the pole of the j-th eigenstate of G ′ ij A does not contribute to the Cauchy integral because the pole is in the lower complex plane.
Let us apply the decomposition (M.91) into the first term. Using the identity (M.92) and performing q 0 integration to pick up the pole q 0 = Ω ǫi,j − Q 1 /2 in [∆G ′ ij R ] i,j , f (q 0 ) is replaced by f (Ω ǫi,j − Q 1 /2). The derivatives with respect to Q 1 or Ω ǫi,j act on it in the leading order approximation. As a result we have (M.104) On-shell contribution from the 5th term in (M.100) is similarly calculated as (M.108) From the second term in (M.100), we have the following off-shell contribu-tion: (M.110) The third term in (M.100) has no leading order contribution: ] j * Π d(eq)jj > <q * ∆G djj Aq = 0 . (M.111) To see this, perform the time and frequency integration. Using (M.92) with Q = 0, we have a factor f (q 0 + Q 2 /2) = f (Ω ǫ ) by picking up the pole at Ω ǫ − Q 2 /2. The derivative with respect to Ω would come from the double poles of ∆G A (G.14), but since we are interested in the region x 0 > y 0 , these poles on the upper complex plane do not contribute. Therefore no derivatives of the distribution function appear. Similarly the 6th term in (M.100) also vanishes in the leading order approximation: −G d(eq)ii Rq * Π d(eq)ii The 5th term in (M.100) gives the off-shell contribution as (M.115) The last line in (M.100) gives off-shell contributions. They are composed of (J.22),(J.40) and (J.41). Using (F.5), we get the following leading order contributions: Summing up these, we have