Flavour Covariant Transport Equations: an Application to Resonant Leptogenesis

We present a fully flavour-covariant formalism for transport phenomena, by deriving Markovian master equations that describe the time-evolution of particle number densities in a statistical ensemble with arbitrary flavour content. As an application of this general formalism, we study flavour effects in a scenario of resonant leptogenesis (RL) and obtain the flavour-covariant evolution equations for heavy-neutrino and lepton number densities. This provides a complete and unified description of RL, capturing three distinct physical phenomena: (i) the resonant mixing between the heavy-neutrino states, (ii) coherent oscillations between different heavy-neutrino flavours, and (iii) quantum decoherence effects in the charged-lepton sector. To illustrate the importance of this formalism, we numerically solve the flavour-covariant rate equations for a minimal RL model and show that the total lepton asymmetry can be enhanced by up to one order of magnitude, as compared to that obtained from flavour-diagonal or partially flavour off-diagonal rate equations. Thus, the viable RL model parameter space is enlarged, thereby enhancing further the prospects of probing a common origin of neutrino masses and the baryon asymmetry in the Universe at the LHC, as well as in low-energy experiments searching for lepton flavour and number violation. The key new ingredients in our flavour-covariant formalism are rank-4 rate tensors, which are required for the consistency of our flavour-mixing treatment, as shown by an explicit calculation of the relevant transition amplitudes by generalizing the optical theorem. We also provide a geometric and physical interpretation of the heavy-neutrino degeneracy limits in the minimal RL scenario. Finally, we comment on the consistency of various suggested forms for the heavy-neutrino self-energy regulator in the lepton-number conserving limit.


Introduction
The observed matter-antimatter asymmetry in the Universe and the observation of non-zero neutrino masses and mixing (for a review, see [1]) provide two of the strongest pieces of experimental evidence for physics beyond the Standard Model (SM). Leptogenesis [2] is an elegant framework that satisfies the basic Sakharov conditions [3], dynamically generating the observed matter-antimatter asymmetry. According to the standard paradigm of leptogenesis (for reviews, see e.g. [4][5][6][7]), there exist heavy Majorana neutrinos in minimal extensions of the SM, whose outof-equilibrium decays in an expanding Universe create a net excess of lepton number (L), which is reprocessed into the observed baryon number (B) through the equilibrated (B + L)-violating electroweak sphaleron interactions [8]. In addition, these heavy SM-singlet Majorana neutrinos N α (with α = 1, . . . , N N ) could explain the observed smallness of the light neutrino masses by the seesaw mechanism [9][10][11][12][13]. Hence, leptogenesis can be regarded as a cosmological consequence of the seesaw mechanism, thus providing an attractive link between two seemingly disparate pieces of evidence for new physics at or above the electroweak scale.
In the original scenario of thermal leptogenesis [2], the heavy Majorana neutrino masses are typically close to the Grand Unified Theory (GUT) scale, M GUT ∼ 10 16 GeV, as suggested by natural GUT embedding of the seesaw mechanism [10][11][12]. In a 'vanilla' leptogenesis scenario [14], where the heavy neutrino masses are hierarchical (m N 1 m N 2 < m N 3 ), the solar and atmospheric neutrino oscillation data impose a lower limit on m N 1 10 9 GeV [15][16][17][18]. As a consequence, such leptogenesis models are difficult to test in foreseeable laboratory experiments. Moreover, these high-scale thermal leptogenesis scenarios, when embedded within supergravity models of inflation, could potentially lead to a conflict with the upper bound on the reheating temperature of the Universe, T R 10 6 -10 9 GeV, required to avoid overproduction of gravitinos whose late decays may otherwise spoil the success of Big Bang Nucleosynthesis [19][20][21][22][23][24][25]. In general, it is difficult to build a testable low-scale model of leptogenesis, with a hierarchical heavy neutrino mass spectrum [4,26].
A potentially interesting solution to the aforementioned problems may be obtained within the framework of resonant leptogenesis (RL) [27][28][29]. The key aspect of RL is that the heavy Majorana neutrino self-energy effects [30] on the leptonic CP-asymmetry become dominant [31,32] and get resonantly enhanced, even up to order one [27,28], when at least two of the heavy neutrinos have a small mass difference comparable to their decay widths. As a consequence of thermal RL, the heavy Majorana neutrino mass scale can be as low as the electroweak scale [33], while maintaining complete agreement with the neutrino oscillation data [1].
A crucial model-building aspect of RL is the quasi-degeneracy of the heavy neutrino mass spectrum, which could be obtained as a natural consequence of the approximate breaking of some symmetry in the leptonic sector. In minimal extensions of the SM, there is no theoretically or phenomenologically compelling reason that prevents the singlet neutrino sector from possessing such a symmetry and, in fact, in realistic ultraviolet-complete extensions of the SM, such a symmetry can often be realized naturally. For instance, the RL model discussed in [27,28] was based on a U(1) L lepton symmetry in the heavy neutrino sector, motivated by superstringinspired E 6 GUTs [34][35][36]. The small mass splitting between the heavy neutrinos was generated by approximate breaking of this lepton symmetry via GUT-and/or Planck-scale-suppressed higher-dimensional operators. The RL model discussed in [29] was based on the Froggatt-Nielsen (FN) mechanism [37] in which two of the heavy Majorana neutrinos, having opposite charges under U(1) FN , naturally had a mass difference comparable to their decay widths. There is a vast literature on other viable constructions of RL models, e.g. within minimal extensions of the SM [38][39][40][41][42][43], with approximate flavour symmetries [44][45][46][47][48][49], with variations of the minimal type-I seesaw [17,[50][51][52][53][54], within SO (10) GUTs [55][56][57][58], within the context of supersymmetric theories [59][60][61][62][63], and in extra-dimensional theories [64][65][66][67][68]. There also exist other variants of the RL scenario, such as radiative RL [69][70][71][72] and soft RL [73,74].
In another important variant of RL, a single lepton-flavour asymmetry is resonantly produced by out-of-equilibrium decays of heavy Majorana neutrinos of a particular family type [75,76]. This mechanism uses the fact that the sphaleron processes preserve, in addition to B − L, the individual quantum numbers X i = B/3 − L i [77][78][79][80][81], where i = 1, 2, 3 is the SM family index and L i is the lepton asymmetry in the ith family. Therefore, it is important to estimate the net baryon number B created by sphalerons just before they freeze out. In particular, a generated baryon asymmetry can be protected from potentially large washout effects due to sphalerons if an individual lepton flavour is out of equilibrium. We refer to such scenarios of RL as resonant -genesis (RL ). In this case, the heavy Majorana neutrinos could be as light as the electroweak scale [33] and still have sizable couplings to other charged-lepton flavours = . This enables the modeling of minimal RL scenarios [76] with electroweak-scale heavy Majorana neutrinos that could be tested at the LHC [82], while being consistent with the indirect constraints from various low-energy experiments at the intensity frontier [83].
Flavour effects play an important role in determining the final lepton asymmetry in RL models. There are two kinds of flavour effects, which are usually ignored in vanilla leptogenesis scenarios, namely: (i) heavy neutrino flavour effects, assuming that the final asymmetry is produced dominantly by the out-of-equilibrium decay of only one (usually the lightest) heavy neutrino, with negligible contributions from heavier species; and (ii) charged-lepton flavour effects, assuming that the flavour composition of the lepton quantum states produced by (or producing) the heavy neutrinos can be neglected and all leptons can be treated as having the same flavour. Neglecting (i) can be justified in 'vanilla' scenarios, because the CP asymmetries due to the heavier Majorana neutrinos are usually suppressed in the hierarchical limit m N 1 m N 2,3 . Moreover, even if a sizable asymmetry is produced by these effects, it is washed out by the processes involving the lightest heavy neutrino [14]. 1 However, for quasi-degenerate heavy neutrinos, as in the RL case, the flavour effects due to the neutrino Yukawa couplings do play an important role [75,86]. In fact, a sizable lepton asymmetry can be generated through CP-violating oscillations of sterile neutrinos [87][88][89][90][91], which is then communicated to the SM lepton sector through their Yukawa couplings.
On the other hand, the lepton flavour effects, as identified in (ii) above, are related to the interactions mediated by charged-lepton Yukawa couplings [92]. Depending on whether these interactions are in or out of thermal equilibrium at the leptogenesis scale, the predicted value for the baryon asymmetry could get significantly modified, as already shown by various partially flavour-dependent treatments [93][94][95][96][97][98]. 2 The lepton flavour effects can be neglected only when the heavy neutrino mass scale m N α 10 12 GeV, in which case all the charged-lepton Yukawa 1 There is an exception to this case depending on the flavour structure of the neutrino Yukawa couplings, when the contribution from the next-to-lightest heavy neutrino decay could be dominant [84,85]. 2 Similar partial flavour effects have also been considered for other variants of leptogenesis models, e.g. with type-II seesaw [99][100][101] and soft leptogenesis [102].
interactions are out-of-equilibrium and the quantum states of all charged-lepton flavours evolve coherently, i.e. effectively as a single lepton flavour, between their production from N α → L l Φ and subsequent inverse decay L l Φ → N α . Here, L l = ( ν lL l L ) T is the SU(2) L lepton doublet (with flavour index l = e, μ, τ ) and Φ is the SM Higgs doublet. For 10 9 m N α 10 12 GeV, the τ -lepton Yukawa interactions are in thermal equilibrium, and hence, the lepton quantum states are an incoherent mixture of τ -lepton and a coherent superposition of electron and muon. Finally, for m N α 10 9 GeV, since the muon and electron Yukawa interactions are also in equilibrium, their impact on the final lepton asymmetry must be taken into account in low-scale RL models. Note that flavour effects also play an important role in the collision terms describing L = 1 scatterings that involve Yukawa and gauge interactions, as well as L = 0 and L = 2 scatterings mediated by heavy neutrinos [33]. Therefore, a flavour-covariant formalism is required, in order to consistently capture all the flavour effects, including flavour mixing, oscillations and (de)coherence. These intrinsically quantum effects can be accounted for by extending the classical Boltzmann equations for number densities of individual flavour species to a semi-classical evolution equation containing a matrix of number densities, analogous to the formalism presented in [103] for light neutrinos. Following this approach, a matrix Boltzmann equation in the lepton flavour space was obtained in [93,98]. Similar considerations were made in [104] to include heavy neutrino flavour effects in a hierarchical scenario. However, in RL scenarios, the interplay between heavy-neutrino and lepton flavour effects are important. With these observations, a fully flavour-covariant treatment of the quantum statistical evolution of all relevant number densities, including their off-diagonal coherences, is entirely necessary. This is the main objective of this long article.
To this end, we derive a set of general flavour-covariant transport equations for the number densities of any population of lepton and heavy-neutrino flavours in a quantum-statistical ensemble. This set of transport equations are obtained from a set of master equations for number density matrices derived in the Markovian approximation, in which quantum 'memory' effects are ignored (see e.g. [105]). We demonstrate the necessary appearance of rank-4 tensor rates in flavour space that properly account for the statistical evolution of off-diagonal flavour coherences. This novel formalism enables us to capture three important flavour effects pertinent to RL: (i) the resonant mixing of heavy neutrinos, (ii) the coherent oscillations between heavy neutrino flavours, and (iii) quantum (de)coherence effects in the charged-lepton sector. In addition, we describe the structure of generalized flavour-covariant discrete symmetry transformations C, P and T , ensuring definite transformation properties of the transport equations and the generated lepton asymmetries in arbitrary flavour bases. Subsequently, we obtain a simplified version of the general transport equations in the heavy-neutrino mass eigenbasis, but retaining all the flavour effects. We further check that these rate equations reduce to the well-known Boltzmann equations in the flavour-diagonal limit.
To illustrate the importance of the effects captured only in this flavour-covariant treatment, we consider a minimal low-scale RL scenario in which the baryon asymmetry is generated from and protected in a single lepton flavour [75]. As a concrete example, we consider a minimal model of resonant τ -genesis (RL τ ) [75], involving three quasi-degenerate heavy neutrinos, at or above the electroweak scale, with sizable couplings to the electron and muon, while satisfying all the current experimental constraints. We show that the final lepton asymmetry obtained in our flavour-covariant formalism can be significantly enhanced (by roughly one order of magnitude), as compared to the partially flavour-dependent limits.
We should emphasize that our flavour-covariant formalism is rather general, and its applicability is not limited only to the RL phenomenon. The flavour-covariant transport equations presented here provide a complete description of the leptogenesis mechanism in all relevant temperature regimes. In addition, this formalism can be used to study other physical phenomena, in which flavour effects may be important, such as the evolution of multiple jet flavours in a dense QCD medium in the quark-gluon plasma (see e.g. [106]), the evolution of neutrino flavours in a supernova core collapse (see e.g. [107]), or the scenario of CPT-violation induced by the propagation of neutrinos in gravitational backgrounds [108]. We have also developed a flavour-covariant generalization of the helicity amplitude technique, and a generalized optical theorem in the presence of a non-homogeneous background ensemble, which may find applications in non-equilibrium Quantum Field Theory (QFT).
It is worth mentioning here that there have been a number of studies (see e.g. [109][110][111][112][113][114][115][116][117][118][119][120]), aspiring to go beyond the semi-classical approach to Boltzmann equations in order to understand the transport phenomena from 'first principles' within the framework of non-equilibrium QFT. Such approaches are commonly based on the Schwinger-Keldysh Closed Time Path (CTP) formalism [121,122]. This real-time framework allows one to derive quantum field-theoretic analogues of the Boltzmann equations, known as Kadanoff-Baym equations [123], obtained from the CTP Schwinger-Dyson equation and describing the non-equilibrium time-evolution of the two-point correlation functions. The Kadanoff-Baym equations are manifestly non-Markovian, accounting for the so-called 'memory' effects that depend on the history of the system. These equations can, in principle, account consistently for all flavour and thermal effects. However, one should note that in order to define particle number densities and solve the Kadanoff-Baym equations for their out-of-equilibrium evolution (as e.g. in the context of leptogenesis), particular approximations are often made. These specifically include quasi-particle approximation and gradient expansion in time derivatives [124]. Moreover, the loopwise perturbative expansion of non-equilibrium propagators are normally spoiled by the so-called pinch singularities [125], which are mathematical pathologies arising from ill-defined products of delta functions with identical arguments. Recently, a new formalism was developed for a perturbative non-equilibrium thermal field theory [126], which makes use of physically meaningful particle number densities that are directly derivable from the Noether charge. This approach allows the loopwise truncation of the resulting transport equations without the appearance of pinch singularities, while maintaining all orders in gradients, thereby capturing more accurately the early-time non-Markovian regime of the non-equilibrium dynamics. An application of this approach to study the impact of thermal effects on the flavour-covariant RL formalism presented here lies beyond the scope of this article.
The rest of the paper is organized as follows: in Section 2, we review the main features of the flavour-diagonal Boltzmann equations. In Section 3, we derive a set of general flavour-covariant transport equations in the Markovian regime. In Section 4, we apply the formalism developed in Section 3 to a generic RL scenario and derive the relevant flavour-covariant evolution equations for the heavy-neutrino and lepton-doublet number densities. In Section 5, we present a geometric understanding of the degeneracy limit in minimal RL scenarios and also discuss an explicit model of RL τ . In Section 6, we present numerical results for three benchmark points, which illustrate the impact of flavour off-diagonal effects on the final lepton asymmetry. We summarize our conclusions in Section 7. In Appendix A, we comment on different forms of the self-energy regulator used in the literature to calculate the leptonic CP-asymmetry in RL models and check their consistency in the L-conserving limit. In Appendix B, we develop a flavour-covariant generalization of the helicity amplitude formalism and describe the flavour-covariant quantization of spinorial fields in the presence of time-dependent and spatially-inhomogeneous backgrounds. In Appendix C, we justify the tensorial flavour structure of the transport equations introduced in Section 3, by means of a generalization of the optical theorem. Finally, in Appendix D, we exhibit the form factors relevant for the lepton flavour violating decay rates discussed in Section 6.

Flavour diagonal Boltzmann equations
The time-evolution of the number density n a of any particle species a can be modeled by a set of coupled Boltzmann equations (see e.g. [127]). Adopting the formalism described in [128,129], this may be written down in the generic form dn a dt + 3H n a = − aX↔Y n a n X n a eq n X eq γ (aX → Y ) − n Y n Y eq γ (Y → aX) , (2.1) where the drift terms on the left-hand side (LHS) arise from the covariant hydrodynamic derivative and include the dilution of the number density due to the expansion of the Universe, parametrized by the Hubble expansion rate H . The right-hand side (RHS) of (2.1) comprises the collision terms accounting for the interactions that change the number density n a . Here, we have summed over all possible reactions of the form aX → Y or Y → aX, in which the species a can be annihilated or created, respectively. If the species a is unstable, it can occur as a real intermediate state (RIS) in resonant processes of the form X → a → Y , which must be properly taken into account in order to avoid double-counting of this contribution from the already considered decays and inverse decays in the Boltzmann equations [128]. At this point, it is important to note that the formalism leading to (2.1) neglects both the coherent time-oscillatory terms, describing particle oscillations between different flavours, and off-diagonal correlations in the matrix of number densities n ab , corresponding to the annihilation of a particle species b and the correlated creation of a particle species a. For this reason, we refer to (2.1) as a set of flavour-diagonal Boltzmann equations. It is useful to summarize the notation and definitions used in (2.1). Firstly, the Hubble expansion rate in the early Universe is given as a function of the temperature T by [127] H (T ) = 4π 3 45 where M Pl = 1.2 × 10 19 GeV is the Planck mass and g * (T ) is the number of relativistic degrees of freedom at temperature T . Throughout our discussions, all species are assumed to be in kinetic (but not necessarily chemical) equilibrium. In this case, the number density of a particle species a is given by where p ≡ d 3 p/(2π) 3 is a short-hand notation for the three-momentum integral, the − (+) sign in the denominator corresponds to particles obeying Bose-Einstein (Fermi-Dirac) quantum statistics, E a (p) = (|p| 2 + m 2 a ) 1/2 is the relativistic energy of the species a, m a being its rest mass, g a = g hel a g iso a is the total degeneracy factor of the internal degrees of freedom, g hel a and g iso a being the degenerate helicity and degenerate isospin degrees of freedom respectively, and μ a ≡ μ a (T ) is the temperature-dependent chemical potential, encoding the deviation from local thermodynamic equilibrium. It will prove convenient in our later discussions to define an in-equilibrium number density n a eq as the limit μ a → 0 in (2.3). We note however that the true equilibrium number density will depend on the equilibrium chemical potential μ eq a , which may not be zero in general.
There are two limits of (2.3) of interest here: (i) the Maxwell-Boltzmann classical statistical limit (E a (p) − μ a )/T 1 in which we can drop the ±1 term in the denominator of (2.3), giving where K n (x) is the nth-order modified Bessel function of the second kind; (ii) the relativistic limit (T m a , μ a ) in which case where σ χ = 1 (3/4) for bosons (fermions), and ζ(x) is the Riemann zeta function, with ζ(3) ≈ 1.20206. Following [29], we define the CP-conserving collision rate for a generic process X → Y and where we have used the shorthand superscript c to denote CP conjugation, and Here, the squared matrix element |M(X → Y )| 2 is summed, but not averaged, over the internal degrees of freedom of the initial and final multiparticle states X and Y . We have introduced an abbreviated notation XY in (2.7) for the phase-space integrals over X and Y . The phase-space measure for the multiparticle state X, containing N X particles, is defined as where δ(x) and θ(x) are the usual Dirac delta and Heaviside step functions, respectively, and N id ! is a symmetry factor in the case that the multiparticle state X contains N id identical particles. In a CPT-conserving theory, the CP-conserving collision rates must obey γ X Y = γ Y X . Analogous to (2.6), a CP-violating collision rate can be defined as [29] which obeys δγ X Y = −δγ Y X , following CPT invariance. The relevant Boltzmann equations for describing leptogenesis are those involving the number densities n N α (with α = 1, . . . , N N ) of the heavy Majorana neutrinos, n L l (with l = 1, . . . , N L ) of the lepton-doublets and n L l of their CP conjugates. When solving the coupled system of firstorder differential equations (2.1) for n N α , n L l and n L l , it is convenient to introduce a new variable z = m N 1 /T . In the radiation-dominated epoch, relevant to the production of lepton asymmetry, z is related to the cosmic time t via the relation t = z 2 /2H N , where is the Hubble parameter (2.2) at z = 1, assuming only SM relativistic degrees of freedom. We also normalize the number density of species a to the number density of photons, defining η a (z) = n a (z)/n γ (z), with n γ given by (2.5) for σ χ = 1 and g γ = 2, i.e.
With these definitions, we write down the flavour-diagonal Boltzmann equations (2.1) in terms of the normalized number densities of heavy neutrinos η N α and the normalized lepton asymmetries δη L l = (n L l −n L l )/n γ as follows [33]: (2.13) where η N eq ≈ z 2 K 2 (z)/2 is the normalized equilibrium number density of the heavy neutrinos, obtained using (2.11) and (2.4) with g N = 2. The various collision rates appearing in (2.12) and (2.13) can be readily understood from the general definitions in (2.6) and (2.9); their explicit expressions in terms of the Yukawa couplings will be given in Section 2.2. Here we have included only the dominant contributions arising from the 1 → 2 decays and 2 → 1 inverse decays of the heavy neutrinos, proportional to the rate γ N α L l Φ , and the resonant part of the 2 ↔ 2 L = 0 and L = 2 scatterings, proportional to γ L k Φ L l Φ and γ L k Φ L c l Φ c respectively. We ignore the sub-dominant chemical potential contributions from the right-handed (RH) charged-lepton, quark and the Higgs fields, as well as the L = 1 Yukawa and gauge scattering terms [33].
Note that for the collision rate pertinent to the heavy neutrino decay in (2.12), we have summed over the lepton flavours, and similarly, for the charged-lepton rate equation (2.13), we have summed over the heavy neutrino flavours; therefore, these are still designated as flavourdiagonal Boltzmann equations. The CP-odd collision rate δγ N α L l Φ in (2.13) can be expressed in terms of the flavour-dependent leptonic CP-asymmetries ε lα and the CP-even collision rate γ N α L l Φ , as follows: δγ N α L l Φ = ε lα k γ N α L k Φ , where ε lα is defined in terms of the partial decay widths (2.14) where Γ N α is the total decay width of the heavy Majorana neutrino N α . Since we are interested in the heavy neutrino decay for temperatures above the electroweak phase transition, where the SM Higgs vacuum expectation value (VEV) vanishes, only the would-be Goldstone and Higgs modes of the Φ-doublet contribute predominantly to the partial decay widths Γ lα and the total decay width Γ N α in (2.14). Fig. 1. The two types of CP-violation due to the heavy Majorana neutrino decay N → LΦ. The notation used here will be explained in Section 3.

Resummed effective Yukawa couplings
The physical CP-violating observable defined in (2.14) receives contributions from two different mechanisms (see Fig. 1): (i) ε-type CP violation due to the interference between the tree-level and absorptive part of the self-energy graphs in the heavy-neutrino decay, and (ii) ε -type CP violation due to the interference between the tree-level graph and the absorptive part of the one-loop vertex. This terminology is in analogy with the two kinds of CP violation in the K 0 K 0 -system (for reviews, see [1,130]), where ε represents the indirect CP violation through K 0 -K 0 mixing, while ε represents the direct CP violation entirely due to the decay amplitude.
The contribution of the self-energy diagrams to the CP-asymmetry can in principle be calculated using an effective Hamiltonian approach, similar to that applied for the K 0 K 0 -system [130]. However, the heavy neutrinos, being unstable particles, cannot be described by the asymptotic (free) in-and out-states of an S-matrix theory [131]. Instead, their properties can be inferred from the transition matrix elements of 2 ↔ 2 scatterings of stable particles, and by identifying the resonant part of the 2 ↔ 2 amplitude that contains the RIS contributions only. This allows one to perform an effective resummation of the heavy-neutrino self-energy diagrams contributing to the ε-type CP-asymmetry [28,29,132]. 3 Neglecting the charged-lepton and light neutrino masses, the absorptive part of the heavy Majorana neutrino self-energy transitions N β → N α can be written in a simple spinorial structure, as follows: where P L,R = (1 4 ∓ γ 5 )/2 are the left-and right-chiral projection operators respectively, and A αβ is the absorptive transition amplitude, summed over all charged-lepton flavours running in the loop: Here h lα is the Yukawa coupling of the heavy neutrino N α with the lepton-doublet L l , and the caret ( ) denotes the fact that (2.16) was derived in a basis in which the heavy Majorana neutrino mass matrix is diagonal. The tree-level decay width of the heavy Majorana neutrino N α is related to the diagonal transition amplitude A αα by To account for unstable-particle mixing effects between the heavy Majorana neutrinos, we define the one-loop resummed effective Yukawa couplings, denoted by (bold-faced Latin) h lα , and their CP-conjugates h c lα , related to the matrix elements M(N α → L l Φ) and M(N α → L c l Φ c ) respectively. This formalism captures all dominant effects of heavy neutrino mixing and CP-violation, and has been shown [29] to be equivalent to an earlier proposed resummation method [28] based on the Lehmann-Symanzik-Zimmermann reduction formalism [139]. Working in the heavy neutrino mass eigenbasis, the resummed effective Yukawa couplings are given by [29,140] where αβγ is the usual Levi-Civita anti-symmetric tensor, m 2 α ≡ m 2 N α is a shorthand notation used here for brevity, and All the transition amplitudes A αβ ≡ A αβ ( h) in (2.18) are evaluated on-shell with p 2 = m 2 N α . The respective CP-conjugate resummed effective Yukawa couplings h c lα can be obtained from (2.18) by replacing the tree-level Yukawa couplings h lα with their complex conjugates h * lα . 4 We will neglect the one-loop corrections to the proper vertices L l ΦN α , whose absorptive parts are numerically insignificant in RL. The partial decay widths Γ αl and Γ c αl appearing in (2.14) can now be expressed in terms of the effective Yukawa couplings h lα and h c lα , and the flavour-dependent absorptive transition amplitudes A l αβ ( h), as follows: Note the explicit dependence of the absorptive transition amplitudes on the effective Yukawa couplings h in (2.20). The total decay width of the heavy neutrino is thus obtained by summing over all lepton flavours: Replacing h by the tree-level Yukawa coupling h in (2.21), we can reproduce the tree-level decay width given by (2.17). Substituting (2.20) in (2.14), the flavour-dependent leptonic CP-asymmetry in RL can be written as αα . (2.22) Note that (2.22) encodes both εand ε -type CP asymmetries, although we simply denote it by ε for brevity. The analytic results for both types of CP-asymmetry and their L-conserving limits for a simplified case will be discussed in Appendix A.

Analytic solutions
It is instructive to derive approximate analytic solutions of the Boltzmann equations (2.12) and (2.29). To do this, we express (2.12) in terms of the non-equilibrium deviation parameter where the K-factors, defined by K α = Γ N α /ζ (3)H N , determine the depletion of the lepton asymmetry due to inverse decays. In deriving (2.23), we have used the analytic expression for the total collision rate γ N α LΦ pertinent to the heavy neutrino decay In the kinematic regime z > z α (2.25) independent of the initial conditions. The collision rates for the L = 0 and L = 2 scatterings are given by [76] where we have used the narrow-width approximation (NWA) for the resummed heavy neutrino propagators in the pole-dominance region, i.e.
since we are only interested in the resonant part of these 2 ↔ 2 scatterings in the RL case. Separating the diagonal α = β RIS contributions from the off-diagonal α = β terms in the sum, (2.13) can be rewritten as [33] n γ H N z where B lα = (Γ lα + Γ c lα )/Γ N α is the heavy neutrino decay branching ratio, and γ X Y ≡ γ X Y − (γ X Y ) RIS denote the RIS-subtracted collision rates, which can be obtained from (2.26) and (2.27) taking α = β. Including only the important RIS-subtracted collision rates proportional to δη L l , and neglecting the terms proportional to δη L k (for k = l) which are numerically small for the minimal RL scenarios [76], we can simplify (2.29) to where we have introduced a flavour-dependent parameter . (2.31) Using the attractor solution (2.25) in the kinematic regime z > z α 1 , (2.30) can be written as where ε l = α ε lα is the total leptonic CP-asymmetry stored in a given lepton flavour l, and K eff l = κ l α K α B lα ≡ κ l K l is the effective washout parameter due to 2 ↔ 2 scatterings mediated by heavy neutrinos. Note that if we only consider the diagonal α = β terms representing the RIS contributions in the sum in (2.31), κ l reaches its maximum value, i.e. κ l = 1 +O(ε 2 l ). On the other hand, in the L l -conserving limit, κ l vanishes at a rate at least equal to that of ε l (see Appendix A). In the regime z > z l 2 ≈ 2(K eff l ) −1/3 , the total lepton asymmetry, dominated by ε-type mixing effects, can be approximated by the analytic solution to (2.32): up to a point z = z l 3 ≈ 1.25 ln(25K eff l ), beyond which the lepton asymmetry freezes out and approaches a constant value δη

Observed lepton asymmetry
Having obtained the net lepton asymmetry l δη L l , the next step is to convert it to the asymmetry in the total baryon-to-photon ratio δη B ≡ (n B −n B )/n γ via (B + L)-violating sphaleron interactions. In a sphaleron transition, an SU(3) c and SU(2) L -singlet neutral object from each generation of the SM is created out of the vacuum [8,[141][142][143]. The operator responsible for sphaleron transitions can be written as where i is the family index; d, e, f are the SU(3) c colour indices; k, l, m, n are the SU(2) L isospin indices; and Q = ( u d ) T L is the SU(2) L quark doublet. The operator O B+L is invariant under both gauge transformations and U(3) flavour rotations. For the case of our interest, the latter freedom can be used to make the charged lepton Yukawa matrix positive and diagonal. Above the electroweak phase transition, all the SM processes, including the sphaleron interactions in (2.34), are assumed to be in full thermal equilibrium, which leads to the following relations among their chemical potentials [78]: where V stands for all vector bosons, u R and d R for up and down-type SU(2) L quark singlets, and e R for SU(2) L lepton singlets in the SM. The total chemical potentials of the baryonic and leptonic doublet fields are then given by Using the relations (2.36) into (2.4), in the linear approximation of μ a /T , we obtain the conversion of the total lepton asymmetry stored in the SM lepton-doublet to the baryon asymmetry 5 assuming a rapid sphaleron transition rate Γ sph H (z = 1). This is valid at temperatures T > T c , where T c is the critical temperature for the electroweak phase transition, given at one loop by [144]  For T < T c , the so-generated baryon asymmetry (2.37) gets diluted by standard photon interactions until the recombination epoch at temperature T 0 . Assuming that there is no source or mechanism for significant entropy release while the Universe is cooling down from T c to T 0 , the baryon number in a comoving volume, n B /s, is constant during this epoch. Here, s = (2π 2 /45)g s T 3 is the entropy density and g s is the corresponding effective number of relativistic degrees of freedom. Thus, the baryon-to-photon ratio, η B = π 4 g s /45ζ (3)s, at the recombination epoch is different from its value predicted by (2.37) due to the change in g s with temperature. The prediction for the current baryon asymmetry is thus given by where we have used g s (T c ) = 106.75 and g s (T 0 ) = 3.91 [127]. The theoretical prediction (2.40) has to be compared with the observed value today, which remains almost unchanged from the end of recombination epoch until the present. The latter can be expressed in terms of directly measurable quantities, namely, the baryon density Ω B h 2 and the primordial 4 He mass fraction Y P , as follows [147]:

Flavour covariant transport equations
As discussed in Section 2, the semi-classical Boltzmann equations (2.1) and, in particular (2.12) and (2.13), do not take into account quantum flavour effects such as the coherent oscillations between different flavours of heavy neutrinos and the quantum-statistical decoherence of flavour off-diagonal matrix number densities. In order to capture these effects consistently, we will derive a set of fully flavour-covariant transport equations for the matrix number densities describing the statistical content of the system. In the next section and as an application of our general formalism, we will consider the specific case of RL, subsequently demonstrating the importance of the flavour effects captured here, but missed in earlier treatments of the subject. Keeping this particular application in mind, we will consider a specific system of N L Dirac lepton isospin doublets L l , N N heavy Majorana neutrinos N α , and an SU(2) L Higgs doublet Φ. Nevertheless, we emphasize that the analysis of this section can be easily generalized to other physical situations involving flavour effects, such as the evolution of jet flavours in a quark-gluon plasma or of neutrino flavours in a core-collapse supernova.
In a general flavour basis, the relevant part of the Lagrangian, involving the heavy Majorana neutrinos, is given by where N R,α ≡ (P R N) α are the RH heavy neutrino fields and Φ ≡ iσ 2 Φ * . Unless otherwise stated, the Einstein summation convention is implied henceforth in the summations over the lepton flavour indices (lower-case Latin) l, m, . . . and the heavy-neutrino flavour indices (lower-case Greek) α, β, . . . . In order to familiarize the reader with the covariant convention used here and throughout this text, we will first discuss the flavour-covariant transformations of the field operators appearing in (3.1). In addition, we will illustrate the flavour-covariant generalizations of the discrete symmetry transformations C, P and T . Subsequently, we will derive general Markovian master equations for the matrix number densities, and use them to describe the statistical evolution of our model system due to the out-of-equilibrium decays and inverse decays of the heavy neutrinos.

Flavour transformations
Under the U(N L ) ⊗ U(N N ) flavour transformations, the lepton fields transform as follows in their fundamental representation: In the physical mass eigenbasis, the Dirac field can be expanded in a basis of plane waves: where b l (p, s, t; t i ) and d † l (p, s, t; t i ) are respectively the interaction-picture particle annihilation and antiparticle creation operators evaluated at the time t = 0. Hereafter, for notational convenience, we will suppress the dependence of the operators on the boundary time t i at which the three pictures of quantum mechanics, viz. Schrödinger, Heisenberg and interaction (Dirac) pictures, are coincident [126]. The index s = ± denotes the two helicity states with the unit spin vector n = ss = sp/|p| aligned parallel or anti-parallel to the three momentum p, respectively. Herein, we have suppressed the isospin index of the lepton doublet. Notice also that we have chosen the normalization of the creation and annihilation operators, such that they have mass dimension −3/2. This choice of normalization is made so that the bilinears of these creation and annihilation operators have the dimension of the number density operator, i.e. mass dimension −3. A discussion of the flavour-covariant Bogoliubov transformations relating this choice of normalization to the more common and manifestly Lorentz-invariant normalization, in which the creation and annihilation operators have mass dimension −1, is included in Appendix B.
It is now important to note that b k (p, s, t) and d † k (p, s, t) transform under the same representation of U(N L ) and so also do b k (p, s, t) ≡ (b k (p, s, t)) † and d †,k (p, s, t) ≡ (d † k (p, s, t)) † . 6 The equal-time anti-commutation relations for these operators are obtained by a flavourtransformation from the mass eigenbasis of the corresponding flavour-diagonal anti-commutators, i.e.
where the rank-2 tensors in (3.7) and (3.8) are defined by means of the flavour transformations (3.2) from the mass eigenbasis, i.e. (3.9) in which M L is the charged-lepton mass matrix and δ l m is the usual Kronecker delta. Since For the Dirac spinors, our notation is such that Full details of the flavour-covariant spinor algebra are given in Appendix B.
In order to write down the flavour covariant expansion of the Majorana field, we recall that in the mass eigenbasis the expansion of a Majorana fermion can be obtained from a Dirac one by imposing the Majorana condition d †,α (k, −r,t) = b α (k, r,t) ≡ a α (k, r,t). (3.11) Since b α (k, r, t) and d †,α (k, − r, t) transform differently under the transformations given in (3.3), this condition cannot be imposed in a general flavour basis. Instead, writing the mass-eigenbasis operators in terms of the ones in a general basis, we obtain U β α b β (k, r, t) = U α γ d †,γ (k, − r, t), where U α β is the flavour transformation that connects the mass eigenbasis to the basis under consideration. We thus obtain the flavour-covariant Majorana condition (3.12) where G αβ denote the elements of a unitary and symmetric matrix. It can be shown that this matrix G is a rank-2 contravariant tensor, 7 which is equal to the identity matrix 1 in the mass eigenbasis. Combining the constraint (3.12) with the expansions (3.7) and (3.8), the flavour covariant expansion of the RH part of the Majorana neutrino field and its Dirac conjugate are given by Notice that the helicity of the v spinors is different from those of the corresponding creation and annihilation operators (see e.g. [149]). The rank-2 tensors in (3.13) and (3.14) can be defined using the flavour transformations (3.3) from the mass eigenbasis, e.g.
The anti-commutation relation for the heavy-neutrino creation and annihilation operators are given by a α (k, r,t), a β k , r ,t = (2π) 3 δ (3) k − k δ rr δ α β . (3.16) From (3.13) and (3.14), we see that the elements G αβ play the role of generalized Majorana creation phases. 8 Finally, the complex scalar field in (3.1) can be expanded as where the interaction-picture creation and annihilation operators for the scalar field satisfy the commutation relations In a general flavour basis, the free Hamiltonians of the lepton doublet and heavy neutrino fields are 7 Performing a flavour transformation U on G αβ defined in (3.12), we get which is the transformation law of a rank-2 contravariant tensor. 8 With the necessary introduction of the matrix G, we may decompose the Majorana field in terms of its RH and LH components as N α = N R,α + G αβ N β L , having definite transformation properties. However, we see that there remains an ambiguity in any such decomposition, since we could equally well have written N α = G αβ N R,β + N α L . (3.20) as can readily be verified by flavour transformations from the mass eigenbasis in which the Hamiltonians are flavour-diagonal. The flavour occupancies and coherences in the evolution of our multiparticle system can be described in terms of flavour-covariant matrix number densities, analogous to the ones for light neutrino flavours introduced in [103]. For the lepton doublets, we define is the infinite coordinate three-volume of the system and the macroscopic time t =t −t i is the interval of microscopic time between the specification of the initial conditions (t i ) and the observation of the system (t). We note in particular the relative reversed order of indices in the lepton and anti-lepton number densities, which guarantees that the two quantities transform in the same representation and thus can be combined to form a flavourcovariant lepton asymmetry. Similarly, we define the heavy-neutrino number densities G αγ a γ (k, r 1 ,t)G βδ a δ (k, r 2 ,t) t , (3.24) and the scalar number densities The total number densities n X (without three-momentum argument) are obtained by integrating (3.21)-(3.25) over the corresponding three-momenta and tracing over helicity and isospin, i.e.
where we have identified explicitly that the traces are taken in isospin space. Analogous definitions are assumed for the antiparticle number densities. Note that all the matrix number densities defined above, as well as the energy matrices (3.9) and (3.15), are Hermitian in flavour space.

Flavour covariant discrete symmetries
It is useful to derive the transformation properties of the discrete symmetries C, P , and T in the flavour-covariant formalism. We assume that the action of these operators on the fermion fields in the mass eigenbasis is the standard one (see e.g. [150]), and find its generalization to an arbitrary flavour basis by means of the appropriate flavour transformations discussed in Section 3.1. In the mass eigenbasis, the action of the unitary charge-conjugation operator U C on elements of the Fock space is given by [150] b l (p, s,t) C . (3.27) Note that the phase convention for the operators used here is in accordance with those used for the spinors in Appendix B. By writing the mass-eigenbasis operators in terms of those in a general basis, i.e. b l (p, s, t) = V m l b m (p, s, t), d †,l (p, s, t) = V l n d †,n (p, s, t), we find the flavour- where we have been required to introduce the matrix G for the charged leptons, analogous to G in (3.12) for heavy neutrinos. Thus, we see that in a flavour-covariant formulation the action of C necessarily involves the appropriate flavour rotation. We find it useful to define the generalized which is a combination of the C-transformation and the appropriate flavour rotation. 9 Thus we see that, in an arbitrary flavour basis, the particle and antiparticle operators are related by a C-transformation, which reduces to the usual charge-conjugation operation in the mass eigenbasis. The action of C on the fermion fields is obtained analogously, i.e.
with C = iγ 0 γ 2 (in the helicity basis), whereas the action of C gives the more familiar result using the fact that G lm G mn = δ l n . Similarly, the parity transformation, given by the unitary operator U P in Fock space, can be generalized straightforwardly from the mass eigenbasis, i.e.
with P = γ 0 . Under CP, the action of the heavy-neutrino interaction Lagrangian (3.1) transforms as where we have introduced the short-hand notation x ≡ d 4 x for the integration over spacetime. Eq. (3.34) defines the CP-transformations of the Yukawa couplings: For a general matrix element, the relation (3.35) can be generalized to (3.36) where Xc ≡ X CP is the generalized CP-transformation of the state X, which can, for instance, be obtained from (3.31) and (3.33).
The action of the time-reversal transformation T is described by an anti-unitary operator A T in the Fock space. Again, starting from the mass-eigenbasis relation we find, because of the anti-linearity of with T = iγ 1 γ 3 . Thus, we may introduce the generalized T -transformations T as follows: Hence, in a general basis, incoming and outgoing states are exchanged by a T operation. This can be seen by transforming the interaction Lagrangian (3.1) under A T from which we obtain Generalizing the above transformations to the matrix elements gives From (3.31) and (3.41), we obtain an important equivalence relation: As a consequence, we have the identity Combining the results (3.36) and (3.43), and using the CPT-invariance of the Lagrangian (3.1), the identity (3.45) allows us to relate the matrix elements as The number density matrices defined in (3.21)-(3.24) have simple transformation properties under C. Since d †,l (p, s, t) = ib l (p, s, t) C , for the lepton number densities (3.21) and (3.22), we have where the transposition on the far RHS of (3.47) is on both flavour and helicity indices. Similarly, for the heavy neutrinos we have a α (k, r) C = −ia α (k, r), and hence, the transformation relation Thus, n N (k, t) T is the C-conjugate of n N (k, t). For Majorana heavy neutrinos, the two C-conjugate quantities are not independent, and are related by Using the C-transformation relations (3.47) and (3.48), we can construct the following quantities with definite C transformation properties: which transform under C as Thus, the quantities δn X are C-"odd", and n N is C-"even", where the quotation marks refer to the fact that this is not meant to be in the usual sense due to the transposition involved. The definite C-properties of the above quantities can be extended to CP, once the total unpolarized number densities defined by (3.26) are considered. Note that we did not define a C-even quantity for lepton number densities (analogous to n N ), since this can be approximated by the equilibrium number density n L eq , i.e.
However, this is not always true for the heavy neutrinos, i.e. n N (t) = n N eq 1, since we must have a departure from thermal equilibrium in order to satisfy the basic Sakharov conditions [3] for the generation of a non-zero lepton asymmetry.
In the heavy-neutrino mass eigenbasis the transformation matrix G reduces to the identity matrix 1, and hence, the transformations C and C are identical for the heavy neutrinos. In this basis, the heavy Majorana neutrino number densities (3.51) and (3.52) reduce to It should be noted that both n N (t) and δ n N (t) are even under the usual charge conjugation operation in the mass eigenbasis, as expected for Majorana fermions 10 : In addition, we note that the total lepton asymmetry δn L (t) ≡ Tr[δn L (t)] is CP-odd in any basis: 10 This is consistent with the C transformations in a general basis, as given by (3.53), due to the transposition involved.
However, under a naive T -reversal transformation N α ↔ N β , we have due to the fact that, in the mass eigenbasis, n N is a symmetric matrix, while δ n N is anti-symmetric.
We will use these definitions to write down the flavour-covariant rate equations for RL in Section 4.

Markovian master equation
In this section, we derive a master equation governing the time evolution of the matrix number densities n X (p, t) in a Markovian approximation. We will work in the interaction picture, beginning from the (picture-independent) definition of the number density in terms of the quantum-mechanical density operator ρ(t; t i ): where the trace is over the Fock space and, for notational simplicity, we leave the momentum dependence implicit. Here we have used the accent (ˇ) to distinguish the quantum-mechanical number density operator ň X (p, t; t i ) from its expectation value n X (p, t), where recall that t = t −t i is the macroscopic time. In addition, for the purposes of this section, it is necessary to reintroduce the explicit dependence of the operators on the microscopic boundary time t i . Differentiating (3.59) with respect to time, we have where we have used the fact that d/dt = d/dt for fixed t i . As we work in the interaction picture, the time evolution of the quantum-mechanical operator ň X (t; t i ) is governed by the free Hamiltonian H X 0 given by (3.19) or (3.20) depending on whether X = L or N . Hence, we use the Heisenberg equation of motion to write the first term in (3.60) as This term generates flavour oscillations in the case of a non-diagonal energy matrix. The second term in (3.60) involves the interaction Hamiltonian, e.g.
As we shall see below, this term will generate the collision terms for the generalized Boltzmann equations (in addition to dispersive corrections). The starting point is the Liouville-von Neumann equation (see e.g. [153]) Rewriting (3.63) in the form of a Volterra integral equation of the second kind, iterating once and subsequently differentiating with respect to time, we obtain the integral form Inserting (3.64) into (3.60), we obtain The first term on the RHS of (3.65) vanishes for the H int term given by (3.62), since it involves the product of an odd number of fields. For the second term on the RHS of (3.65), we may use the cyclicity of the trace to obtain When used in (3.60), this gives an exact and self-consistent time-evolution equation, which captures the entire evolution of the system, including any non-Markovian memory effects. We now perform a set of Wigner-Weisskopf approximations [154] to obtain the leading-order Markovian form of (3.60). Let us define the t-dependent function Inserting the Fourier transforms of H int (t; t i ) and F (t; t i ) with respect to t −t i and t −t i in (3.66), we obtain Making the change of variables ω = ω − ω , this may be recast in the form As long as F (ω − ω ) remains dynamical on inverse Fourier transformation, i.e. ω = ω , the dominant contribution to the integral (3.69) occurs for t ∼t . We may therefore replace . We now make the change of variables t =t −t and take the limit t i → −∞. Herein, we assume that the statistical evolution is slow compared to the quantum-mechanical evolution, such that the system remains out-of-equilibrium in spite of the quantum-mechanical boundary time being in the infinitely distant past. With this approximation, we replace the interaction-picture creation and annihilation operators in H int (t; t i ) and ň X (t; t i ) by their asymptotic 'in' counterparts via where Z = 1 + O(h) is the wavefunction renormalization factor. Notice that in the replacement (3.70), we must account for the free phase evolution of the interaction-picture operators.
The contribution I 2 now takes the form Performing the t integration, we find where P denotes the Cauchy principal value. We are then able to write the result where objects constructed from asymptotic operators depend only on the time t and the O(h) corrections from the wavefunction renormalization have been neglected. The first term in (3.73) can be identified as the collision term C X and the second term represents the dispersive selfenergy corrections arising from vacuum contributions (Lamb shift) and thermal contributions (Stark shift), which we neglect in the following. A discussion of our resummation scheme in relation to self-energy corrections is included at the beginning of Section 4.
Restoring the implicit momentum dependence of the number densities, we finally obtain the leading-order master equation in the Markovian approximation 11 : This is our central equation governing the time-evolution of the particle number densities with arbitrary flavour content. 12 With (3.74) no longer in integro-differential form [cf. (3.64)], we are now free to specify the initial conditions at any finite macroscopic time t 0 . It remains the case however that the macroscopic time t 0 = 0 corresponds to the microscopic time t =t i → −∞. We note that although a similar form of (3.74) was also used in [103] to describe the timeevolution of active neutrinos in a thermal bath, the full flavour structure contained in (3.74) to describe the simultaneous time-evolution of multiple species, e.g. heavy neutrinos and SM leptons as in the context of leptogenesis, was not discussed before in the literature. In order to 11 There are two major assumptions in this approximation: (i) separation of time scales, i.e. the QFT processes governed by H int occur at much smaller time scales as compared to the coarse-grained statistical evolution governed by ρ(t); and (ii) molecular chaos, i.e. the velocity correlations that may form between different species in the QFT processes are lost on time-scales relevant for the statistical evolution, so that the background can be factorized at all times. 12 This formalism is sometimes called the 'density matrix formalism' in the literature. In our opinion, this terminology is misleading, since n X (k, t) is actually a matrix of densities [103], which should be distinguished from the quantumstatistical density matrix ρ(t). Such confusion could potentially lead to incorrect results, since there is a crucial sign difference in the time evolution equations for the two quantities, as can be seen by comparing (3.63) with the first term on the RHS of (3.74). Therefore, we avoid referring to our approach as the 'density matrix formalism'. elucidate these flavour effects, we will explicitly derive fully flavour-covariant transport equations for the system described by the Lagrangian (3.1) in the following subsection.

Transport equations
Using the expressions (3.19) and (3.62) for the free and interaction Hamiltonians respectively, we can explicitly calculate the oscillation and collision terms in (3.74). Specifically, we obtain the following evolution equations for the charged lepton and anti-lepton number densities, specifically where commutators carrying flavour indices are understood to be commutators in flavour space. The collision terms are given by where we have suppressed the overall momentum dependence and introduced the compact notation It is important to emphasize that our flavour-covariant formulation requires new rank-4 tensors in flavour space: (i) the statistical number density tensors and (ii) the absorptive tensors In (3.83) and (3.84), the rank-4 delta function of on-shell four-momenta originates from the integration of tensor exponentials, such as and is defined as a linear combination of ordinary delta functions in terms of the appropriate flavour transformation from the mass eigenbasis, as defined in Section 3.1, i.e.
The absorptive tensors (3.83) and (3.84), obtained from the Markovian master equation (3.74), represent the contributions from decays and inverse decays of the heavy Majorana neutrinos to the statistical evolution of the system (see Section 4.1 and Figs. 2(a)-2(b)). As shown in Appendix C, these rank-4 objects can be interpreted in terms of the unitarity cut of the partial one-loop heavy-neutrino self-energies, using a generalized optical theorem (see Fig. 3). This justifies the necessity of the tensorial structure of these objects, and also the form of the 2 ↔ 2 scattering terms that will be included later in the evolution equations when directly applied to the RL scenario in Section 4. This formalism can be generalized to include higher order processes involving multiple flavour degrees of freedom, e.g. LN ↔ Le R and LN ↔ LN , by introducing the corresponding rate tensors of rank 6 and higher. Proceeding analogously for the heavy neutrino number densities, we find the evolution equa- with the heavy-neutrino collision terms given by where the tensor contraction is analogous to (3.79) and (3.80), with the role of charged-lepton and heavy-neutrino indices exchanged. Note the appearance of the matrix G in the transport equations (3.89) and (3.90), and the transposition of both flavour and helicity indices in the corresponding collision terms. One should however remember that (3.89) and (3.90) are not independent, because of the relation (3.48). Note also that the transport equations have an internal structure in isospin space, which has been suppressed here for brevity. In Section 4, when we derive the rate equations for the total number densities, we will explicitly trace over these degenerate isospin degrees of freedom. As a final remark, we point out that the flavour-covariant formalism developed so far can also be applied more generally to describe quantum coherences between species with different SM quantum numbers, e.g. L l and N α . In this case, we may need to study the evolution of the number densities [n LN ] l α , which transforms as a rank-2 tensor in the flavour space U(N L ) ⊗ U * (N N ), corresponding to the correlated creation of L l and the annihilation of N α . The evolution equation would still have the generic form (3.75), with collision terms of the form (3.77). However, we will neglect these effects in our discussions, since they are not expected to play an important role for the RL scenarios under study.

Application to resonant leptogenesis
As already discussed in Section 2.1, there are two types of CP-violation possible in the outof-equilibrium decay of the heavy Majorana neutrino. In the limit when two (or more) heavy Majorana neutrinos become degenerate, the ε-type CP-violation can be resonantly enhanced by several orders of magnitude, as compared to the ε -type CP-violation [27,28] (see also Appendix A). In this case, finite-order perturbation theory breaks down and it is necessary to perform a consistent field-theoretic resummation of the self-energy and vertex corrections, as shown schematically in Fig. 1. However, such resummation can only be performed in a closed algebraic form when working in the mass eigenbasis and at full thermodynamic equilibrium in a Markovian approximation. This fact is illustrated explicitly in Appendix B, where we derive the most general flavour-covariant propagators in a non-homogeneous background within the Schwinger-Keldysh CTP formalism. Therein, we show that, when out-of-equilibrium (where flavour-coherences must be permitted), canonical quantization necessarily leads to off-diagonal propagators in flavour space. As a result, the inversion of the Schwinger-Dyson equation for the resummed propagators contains an infinite nesting of convolution integrals, which does not collapse to the usual algebraic equation of resummation.
As identified in [126], the loopwise truncation of quantum transport equations is two-fold, i.e. the transport equations may be truncated both spectrally and statistically. The spectral truncation corresponds to choosing those objects that will be counted as particulate degrees of freedom; in our case, we will choose to count the number densities of spectrally-free on-shell particles. The statistical truncation, on the other hand, corresponds to the restriction of the set of processes that drive the macroscopic evolution of the system. It is the latter that we will treat non-perturbatively in order to consistently account for the εand ε -type CP-violation effects. Recalling the argument above, we ensure that this resummation can be performed algebraically by first flavour-rotating to the heavy neutrino mass eigenbasis, maintaining the Markovian approximation used in Section 3.3 and resumming only zero-temperature contributions, thereby neglecting thermal loop effects [155]. As we shall see later in this section (see Section 4.4), the omission of these thermal effects is appropriate in the classical statistical limit, as long as one accounts systematically for both the zero-temperature and thermal RIS contributions. Lastly, we will assume that the heavy-neutrino helicity states are fully decohered and equally populated [29].
Hence, we proceed by replacing the tree-level neutrino Yukawa couplings by their resummed counterparts in the transport equations given in Section 3. 4 where h m β ≡ h mβ in the mass eigenbasis. We emphasize however that the resummation itself must be performed in the mass eigenbasis only. Further justification for the choice of this basis will be given in Section 5.1, where we show that in the degenerate symmetry limit of the minimal RL model, it is important to choose the correct basis in order to get meaningful results.

Rate equations for decay and inverse decay
In order to obtain the rate equations from the general transport equations (3.75), (3.76), (3.89) and (3.90), we need to impose kinetic equilibrium. This can be ensured throughout the evolution of the system by assuming that the elastic scattering processes rapidly change the momentum distributions on time-scales much smaller than the statistical evolution time of the particle number densities. This approximation is valid as long as the mass splittings between different flavours inside thermal integrals are much smaller than the average momentum scale, i.e. |k| ∼ T |m N α − m N β |. In this regime, the momentum distributions governed by the elastic processes are flavour singlets [156,157]. Using this approximation, we introduce an average mass for N N quasi-degenerate heavy neutrinos: to be used within the thermal integrals. Correspondingly, we may introduce an average energy E N (k) = (|k| 2 + m 2 N ) 1/2 . Furthermore, to simplify the general transport equations given in Section 3.4, we take the classical statistical limit in which (3.81) and (3.82) can be approximated as The spinor traces appearing in (3.83) and (3.84) can be simplified if we assume the chargedleptons to be massless, neglecting their thermal masses. In this limit, n L is the number density matrix for the LH lepton doublets L l and one helicity index for the charged leptons can be dropped in the spinor traces, thus yielding (see also Appendix B) where − (+) corresponds to the helicity of the massless LH lepton (RH anti-lepton).
With these approximations, the heavy-neutrino kinetic-equilibrium number density, given by the flavour-covariant generalization of (2.4) (with the chemical potential now being a rank-2 tensor), can be approximated as where n N (t) is the total heavy neutrino number density, as defined in (3.26), and n N eq is the equilibrium number density given by (2.4) setting μ N = 0 and g N = 2 for the two heavy-neutrino helicity states. An analogous expression can be obtained for the charged lepton number density where E L is the average energy of the lepton doublets, and in the massless limit, n L eq is the equilibrium number density given by (2.5) with g L = 2 for the two isospin states. Notice that the factor g L is not present in (4.6), and will appear explicitly only after the trace over isospin is performed [cf. (3.26)]. Finally, for the Higgs number density, we assume an equilibrium distribution n Φ (q, t) = e −E Φ (q)/T . Throughout the remainder of this article, we suppress the t -dependence of the number densities for notational convenience.
We may now integrate both sides of (3.75) and (3.76) over the phase space and sum over the degenerate helicity and isospin degrees of freedom. The resulting rate equations for the total number densities of the charged lepton and anti-lepton doublets, accounting for the decay and inverse decay of the heavy neutrinos, can be written in the form where E L is the thermally-averaged effective lepton energy matrix with M L (T ) being the thermal mass matrix of the lepton doublets. Note that the 3T -term on the RHS of (4.9) is isotropic in flavour space, commutes with the number densities, and therefore, does not give any contribution to the rate equations (4.7) and (4.8). The 1 → 2 and 2 → 1 collision rates appearing in (4.7) and (4.8) are derived from the rank-4 absorptive tensors (3.83) and (3.84). Replacing the tree-level Yukawa couplings h l α appearing in (3.83) and (3.84) with the resummed ones h l α , the 1 → 2 collision rates can be explicitly written as which are the flavour-covariant generalizations of the collision rates defined in (2.7). The 2 → 1 collision rates in (4.7) and (4.8) are related to the 1 → 2 rates (4.10) and (4.11) by virtue of CP T The corresponding rank-2 collision rates within the anti-commutators in (4.7) and (4.8) are obtained from the corresponding rank-4 tensors by contracting the heavy-neutrino flavour indices, e.g. (4.14) In Appendix C, we present an alternative derivation of these collision rates by considering a flavour-covariant generalization of the optical theorem in the presence of a statistical background. Therein, the necessity of the rank-4 flavour structure of these collision rates is further justified. This is illustrated in Fig. 2 for the in-medium production of heavy neutrinos in a spatially-homogeneous statistical background of lepton and Higgs doublets. The production rates in the thermal plasma can be better understood from the unitarity cut of the partial one-loop heavy-neutrino self-energy graph, as shown in Fig. 3. Imposing kinetic equilibrium as given by (4.6), we obtain tree-level thermally-averaged heavy-neutrino production rates for the processes L c Φ c → N and LΦ → N [cf. (C.51) and (C.52)], which are exactly the same as those obtained in (4.10) and (4.11), respectively.
Analogous to the charged-lepton case, we obtain the flavour-covariant rate equations for the total number densities of heavy neutrinos from the general transport equations (3.89) and (3.90), as follows: 16) where E N is the thermally-averaged effective heavy-neutrino energy matrix, defined analogous to (4.9). The thermally averaged collision terms C and C in (4.15) and (4.16) can be defined analogous to E N , and are explicitly given by  where the rank-4 collision rates are given by (4.10) and (4.11), and the corresponding rank-2 objects appearing in (4.17) and (4.18) are obtained by contracting the charged-lepton indices, e.g.
Using the expressions (4.10) and (4.11), we can define the flavour-covariant generalizations of the CP-even and CP-odd quantities in (2.6) and (2.9), which now have definite transformation properties under CP: The corresponding rate equations for the CP-"even" and -"odd" number densities [cf. (3.53)] are derived from (4.7), (4.8), (4.15) and (4.16): where we have kept terms only up to O(μ a /T ) and O(h 4 ), except the last term on the RHS of (4.24), which is the only washout term for the lepton asymmetry, and appears at O(h 6 ). In (4.22) and (4.23), we have defined, for a given Hermitian matrix A = A † , its generalized real and imaginary parts, as follows: Observe that in the heavy-neutrino mass eigenbasis, the definitions (4.25) and (4.26) reduce to the usual real and imaginary parts: In obtaining (4.22) and (4.23), we have used the relations LΦ ] l m , appear with the wrong sign and do not lead to the correct equilibrium behaviour, since, by construction, there are no washout terms, induced by the 2 ↔ 2 scatterings, introduced so far. It is well known that the inclusion of these scattering terms (with the RIS contribution subtracted), changes the sign of these inverse decay terms and gives the correct equilibrium limit [29,128]. In Section 4.4 we will explicitly show this result in a flavour-covariant manner by including in the rate equations the RIS-subtracted collision rates for scattering in the presence of a statistical background, as illustrated in Fig. 4. For the moment, we take this result at face value, and correct the signs 'by hand', to be able to qualitatively discuss some important physical phenomena, before we include the scattering terms. Finally, we also take into account the Hubble expansion of the Universe and change the independent variable t in favour of z = m N /T , to write down the rate equations in terms of the normalized number densities, introduced in Section 2: In the last term on the RHS of (4.31), we have used η L eq = 3/4, which follows from (2.5) and (2.11).

Lepton asymmetry via heavy neutrino oscillations
The rate equations (4.29)-(4.31), contain two physically distinct mechanisms for the generation of lepton asymmetry. One is the standard T = 0 εand ε -type CP violation given by (2.22) due to the mixing and decay of the heavy Majorana neutrinos. We have taken this into account by means of the one-loop resummed effective Yukawa couplings defined in (2.18), which appear in the collision rates (4.10) and (4.11). This is the only source of lepton number asymmetry in the flavour-diagonal Boltzmann equations (2.12) and (2.13). However, the flavour-covariant rate equations (4.29)-(4.31) also include a second source for the asymmetry, due to the CP-violating oscillations of the on-shell heavy neutrinos in the thermal bath. This originates from the sequence of an on-shell production of heavy neutrinos in the bath due to inverse decays, which can oscillate between different flavours in the bath and then decay back into charged-leptons. Formally, this process has the same structure as the scattering diagrams in Fig. 4. While the T = 0 QFT processes are taken into account by the resummation of the Yukawa couplings [28,29], the oscillation phenomenon corresponds to the thermal part of the intermediate heavy-neutrino propagator. Thus, our flavour-covariant approach captures the leading order effect of a complete thermal resummation procedure that would generalize the analysis in [29].
In this section, we will present a qualitative analysis of the heavy-neutrino oscillation phenomenon in the RL scenario, and show that this mechanism contributes to the total lepton asymmetry at order O(h 4 ) around z = 1. Note that, although conceptually similar, this is qualitatively as well as quantitatively different from the phenomenon first proposed in [87], and studied in [88][89][90] for the light sterile neutrino case, where the final lepton number asymmetry is of order O(h 6 ), as also recently stressed in [158]. In order to extract the heavy-neutrino oscillation effect from the flavour-covariant rate equations (4.29)-(4.31), we consider a simplified case with the tree-level Yukawa couplings (instead of the resummed ones), thus artificially 'switching off' all T = 0 εand ε -type CP violation effects. In this case, we can drop the CP-"odd" δγ collision terms in the rate equations (4.30)-(4.31), which further simplify in the mass eigenbasis to give The rate equation for η N is still given by (4.29), specialized to the mass eigenbasis. Notice that, in the mass eigenbasis, the flavour rotation matrices G = G = 1, and therefore, we do not need to distinguish between upper and lower flavour indices in (4.32) and (4.33). It is useful to perform a time-stepping analysis to see the infinitesimal time evolution of the total lepton asymmetry.
We start with an incoherent diagonal heavy-neutrino number density matrix, and a zero initial lepton asymmetry at some z = z in : Interference in the inverse decays will generate off-diagonals in η N at O(h 2 ). However, we see from (4.29), (4.32) and (4.33) that only η N , and not δ η N or δ η L , receives these contributions. This is due to the fact that the processes (the c α being complex coefficients in the linear combination) have to proceed at the same rate for CP-conserving inverse decays, if no initial lepton asymmetry is present, and hence, only Re( η N ) = η N can be generated. Therefore, after a small time interval, at z = z in + δz 1 we have Now from (4.32), we see that heavy-neutrino oscillations will induce imaginary parts of η N , and hence, a non-zero δ η N [cf. (4.32)] due to the off-diagonals of η N in (4.36). Thus, at a later time z = z in + δz 2 , the number densities will be Finally at z = z in +δz 3 , interference in the O(h 2 ) decays will create a non-zero lepton asymmetry of order O(h 4 ) from the δ η N term in (4.33): The physical reason for this is that the CP-conjugated processes α c α N R,α → LΦ, α c * α Nc R,α → LcΦc, (4.39) are respectively proportional to the number densities η N and η N = ( η N ) * , which now differ by Thus, the O(h 4 ) contribution to the total lepton asymmetry is due to a sequence of the coherent heavy-neutrino inverse decays, oscillations and decays. These effects get enhanced in the same regime as the resonant T = 0 ε-type CP violation, namely, for z ≈ 1 and m N ∼ Γ N α [28]. For z 1, this effect is suppressed by the small mass of the Majorana neutrinos, and for z 1 the inverse decays are frozen out and ineffective to create an asymmetry. Similarly, if the heavyneutrino mass-splitting m N is too large compared to Γ N α , the oscillations are averaged out during the decay time scale, whereas if it is too small the oscillations proceed too slowly and δ η N produced thereof is constantly washed out.
In Section 6.2, we will show quantitatively that, for the RL τ model under consideration there, the lepton asymmetry generation via the heavy-neutrino oscillation phenomenon discussed above is of the same order as the one due to their mixing in the vacuum, and hence, leads to an enhanced total lepton asymmetry (even in the charged-lepton flavour diagonal case) compared to that predicted by the flavour-diagonal Boltzmann equations discussed in Section 2.

Decoherence in the charged lepton sector
In this subsection, we will include in the rate equations the effect of the charged-lepton Yukawa couplings, described by the interaction Lagrangian L y = y l kL k Φe R,l + H.c., (4.40) where e R,l ≡ l R (with l = e, μ, τ ) and the Yukawa couplings are real and diagonal in the charged-lepton Yukawa eigenbasis, i.e. y l k = y l δ l k . The number densities of the RH leptons and anti-leptons, n R and n R respectively, can be defined analogous to n L and n L [cf. (3.21) and (3.22)]. The processes involving the charged-lepton Yukawa couplings will be responsible for the decoherence of the LH leptons into the mass eigenbasis. However, due to the interaction of the charged-leptons with the heavy-neutrinos [cf. (3.1)], non-zero off-diagonal elements are necessarily induced in the charged-lepton number density matrix, which tend to recreate a coherence between the charged-lepton flavours. Thus, there is a competition between the coherence effect induced by the heavy-neutrino Yukawa couplings and the decoherence effect due to the chargedlepton Yukawa couplings, and for reasonably large neutrino Yukawa couplings, the coherence effect could be significant, as we will see explicitly in Section 6.2. In particular, we will find that the decoherence is not complete in the temperature range relevant for the production of the asymmetry in the RL τ scenarios with 200 GeV m N 2 TeV. Hence, it is important to include the off-diagonal lepton number densities, and the dynamics of their decoherence effects, in the rate equation for the lepton asymmetry.
Let us first consider the contribution of thermal Higgs decays and inverse decays In the above, we have defined the charged-lepton decoherence and back-reaction rates where the flavour-singlet term A(p, k, q), whose explicit form is not needed here, contains the relevant kinematic factors.
In the would-be mass eigenbasis for the charged-leptons in which the charged-lepton Yukawa coupling matrix is diagonal, the diagonal entries of (4.42) have the form where the index l is not summed over, and we have assumed n R (k) to be diagonal, neglecting higher-order phenomena involving heavy-neutrino Yukawa couplings. Since the evolution equations of n R and n Φ contain the same term (4.45), and since the rate of the process (4.41) is much larger than the Hubble rate for the relevant time period, we can safely assume that their fast evolution always guarantees the chemical (as well as kinetic) equilibrium for this reaction, i.e. μ Φ = μ L,l +μ R,l , in addition to the decoherence of n R . Therefore, the evolution equations for Φ and ē R need not be considered explicitly, and instead, we can use the relevant detailed balance condition, or the KMS relation [159][160][161], to solve the rate equation. For the process (4.41) under consideration, the KMS relation is simply given by The first condition comes from the fact that Γ back dec and Γ dec are simultaneously diagonal in the charged-lepton mass eigenbasis. The second condition is the generalized KMS relation (4.46) involving any species in chemical and kinetic equilibrium with the LH charged leptons in the mass eigenbasis, whose fast evolution ensures the vanishing of the overall statistical factors that multiply the large diagonal rates in the rate equation (4.42).
The charged-lepton Yukawa contributions to the rate equation for anti-lepton number density n L will be analogous to that given in (4.42). To obtain the corresponding contribution to the rate equation for the lepton asymmetry, we use the same set of approximations as in Section 4.1, and in particular, the kinetic-equilibrium number density (4.6). Taking into account the expansion of the Universe, we finally obtain where γ dec and δγ back dec are the CP-even and -odd thermally-averaged decoherence and backreaction rates, respectively. Here we have ignored the sub-dominant {η L eq , δγ dec } term, which depends on the asymmetry in the RH charged-lepton sector that is assumed to be small compared to the asymmetry in the LH sector. The CP-even rate can be expressed in terms of the charged-lepton thermal width as γ dec = Γ T n L eq . In the mass eigenbasis, the thermal width is given by Γ T = diag{Γ T ,l }, and has been calculated explicitly in [144], taking into account the inverse Higgs decays and the relevant fermion and gauge scatterings: where y t is the top quark Yukawa coupling, and x = M H (T )/T = zM H (T )/m N , M H (T ) being the Higgs thermal mass. Note that while calculating the final rates for the processes involving the charged-lepton Yukawa couplings, it is important to take into account their thermal masses, which control the phase space suppression for the decay and inverse decay of the Higgs boson [144]. Additionally, note that all the chemical potentials can be consistently neglected in the calculation of the rate, as long as we satisfy the generalized KMS relations given by (4.47). After thermal averaging, (4.47) can be written as (4.50) These equations ensure that the detailed balance conditions are satisfied, without having to resort to following the evolution of all the SM species involved in the charged-lepton Yukawa-mediated processes.

Scattering terms
In this section, we will describe the flavour-covariant generalization of the subtraction of the so-called RIS contributions present in the L = 2 and L = 0 scattering terms (see Fig. 4). Specifically, we will show how the sign of the CP-"odd" inverse decay terms in (4.23) and (4.24) is flipped, so that the correct approach to equilibrium is restored. Moreover, we will illustrate that it is necessary to account for thermal corrections in the RIS contributions when considering off-diagonal flavour correlations.
Following [29], we first consider the case of a single scalar sneutrino N , and write down the thermally corrected scattering amplitude T (LΦ → LcΦc) as where we have used on-shell renormalization scheme and have neglected thermal dispersive corrections. Notice that, since the Lorentz-covariance of thermally-corrected self-energies is broken, the absorptive parts of the time-ordered and retarded equilibrium CTP self-energies Im By virtue of the fluctuation-dissipation theorem for the equilibrium self-energies, we have the relation where ε(p 0 ) ≡ θ(p 0 ) − θ(−p 0 ) is the generalized signum function. In the pole-dominance region, using the NWA given by (2.28), we obtain Generalizing the result in (4.53) to fermions, we see that the RIS contribution to the LΦ → LcΦc scattering terms contains the following combination of statistical factors: (4.55) In the classical statistical limit, (4.55) becomes The first two terms in (4.56) arise from the T = 0 part of the RIS subtraction and contribute to the charged-lepton transport equations (3.75) and (3.76). The third term arises from the thermal correction to the RIS subtraction and contributes to the heavy-neutrino transport equations (3.89) and (3.90). The rank-4 scattering rates introduced here can be derived by means of the generalized optical theorem given in Appendix C. For consistency with the previous calculations, only their resonant parts must be kept and these are listed in (C.63)-(C.66). Using the NWA for the heavy neutrino propagator, and the same approximations as for the decay and inverse-decay terms in Section 4.1, we obtain the CP-"even" collision rates to the unitarity of the scattering matrix [128], we obtain the RIS-subtracted CP-"odd" collision rates

Contributions to lepton transport equations
and the corresponding RIS-subtracted CP-"odd" quantity (4.70) The flavour structure of the γ -terms in (4.62) and (4.69) can be understood diagrammatically, as shown in Figs. 3 and 5, from the unitarity cuts of partial self-energies, obtained by virtue of a generalized optical theorem (see Appendix C). (4.75) Note that, thanks to the first identity (4.72), the sign of the inverse-decay term (first term in the RHS) is now flipped with respect to that in (4.24), as anticipated at the end of Section 4.1, and it guarantees the correct approach to equilibrium. The remaining identities (4.73) and (4.74) are important to guarantee the consistency of the formalism: following [128], we have described processes like LΦ → N → LΦ (with an on-shell N ) statistically, i.e. as the successive statistical evolution of the number density n N first due to an inverse decay and then a decay process. The RIS-subtracted scattering term is considered in order to avoid double-counting. However, (4.73) and (4.74) allow us to write the washout term in terms of a complete (including RIS) scattering rate, with no inverse decay rate at all, thus describing resonant processes like LΦ → N → LΦ field-theoretically. Both these descriptions lead to the same result, as shown in [76].

Contributions to heavy-neutrino transport equations
The L = 0 and L = 2 scattering contributions to the heavy-neutrino rate equations are given by where, for notational simplicity, we have defined and the charged-lepton indices are understood to be traced over [cf. (C. 67)]. As explained in the beginning of this subsection, these terms arise due to the last term on the RHS of (4.56). Following the same procedure as in the charged-lepton case discussed above, we obtain the relevant contributions to the rate equations for the CP-"even" and -"odd" number densities: This contribution to δ n N , when added to the decay and inverse decay contribution given by (4.23), flips the sign of the inverse decay term with respect to that given in (4.23), as expected in order to achieve the correct equilibrium behaviour.

Final rate equations
Here we put together the various contributions from heavy-neutrino decays and inverse decays discussed in Section 4.1, from processes involving charged-lepton Yukawa interactions [cf. Section 4.3], and from L = 0 and L = 2 scatterings via heavy neutrino exchange [cf. Section 4.4]. Finally, taking into account the expansion of the Universe, the following set of manifestly flavour-covariant rate equations is obtained for the CP-"even" number density matrix η N and CP-"odd" number density matrices δη N and δη L : in the first line of (4.84). We should stress here that this phenomenon of coherent oscillations is an O(h 4 ) effect on the total lepton asymmetry, and so differs from the O(h 6 ) mechanism proposed in [87] (see Section 4.2). (iii) Decoherence effects due to charged-lepton Yukawa couplings, described by the last line of (4.84). Our description of these effects goes along the lines of [93], which has been generalized here to an arbitrary flavour basis.
As an application, we will use the rate equations (4.82)-(4.84) in Section 6.2 for the numerical evaluation of the lepton asymmetry in the RL τ model under consideration there. We will also derive approximate analytic solutions of these general rate equations in Section 5.3. Finally, we note that taking the limit in which the number densities are diagonal, i.e.

Minimal resonant -genesis model
In this section, we discuss the basic theoretical framework underlying the minimal RL model in which the lepton asymmetry is dominantly generated and stored in an individual lepton flavour [75]. We start with the heavy Majorana neutrino sector Lagrangian given by (3.1). Note that above the scale of the electroweak phase transition, only the singlet neutrinos are massive, whose origin must lie in some ultraviolet-complete extension of the SM. Within the minimal RL setup, the masses of all these heavy neutrinos N α (α = 1, . . . , N N ) The compatibility of the light neutrino masses generated via the seesaw mechanism with the solar and atmospheric neutrino oscillation data requires that, for electroweak-scale heavy neutrinos with m N ∼ O(100) GeV, the norm of the Yukawa coupling matrix must be much smaller than unity. Given (5.2), this implies that the norm of the O(N N )-breaking matrix M N must be small compared to m N , i.e. M N /m N 10 −7 . As shown in [76], a M N of the required order can indeed be generated radiatively for RL models.
A non-zero total lepton asymmetry summed over all flavours can be created, if and only if the CP-odd quantity [28] does not vanish for a finite non-zero interval of RG scales. In general, the total number of all nontrivial CP-violating phases in a model with N L weak iso-doublets and N N neutral iso-singlets is N CP = N L (N N − 1) [162,163]. This results in N CP CP-odd quantities, analogous to the one defined in (5.3), which generally mix under RG effects. However, after summing over all final state lepton flavours occurring in the heavy neutrino decays, only one CP-asymmetry remains, which is odd under the generalized CP-transformations discussed in Section 3.2. Using the definition of CP in (5.3) and the heavy neutrino mass matrix (5.1), and taking into account the RG evolution, one can find the necessary and sufficient condition for the generation of a non-zero total CP-asymmetry in the minimal RL setup. The results of this exercise, including RG effects except for the charged-lepton Yukawa couplings, are summarized in Table 1 for different choices of N L and N N . Table 1 Total CP-asymmetry in the minimal RL model with N L lepton iso-doublets and N N neutral iso-singlets, including the RG effects, except for the charged-lepton Yukawa couplings. Here / 0 means CP = 0.

Geometry of the degeneracy limit
In To illustrate this point, let us first consider a simple case with one charged-lepton flavour (N L = 1) and two heavy-neutrino flavours (N N = 2). In this case, the tree-level heavy-neutrino mass matrix M N is symmetric under O(2). Let us therefore define the tree-level Yukawa coupling as a two-dimensional complex vector in the heavy neutrino mass eigenbasis {N 1 , N 2 }: h ≡ ( h l1 , h l2 ), which can also be written in terms of two real vectors, i.e. h = Re( h) + i Im( h) ≡ a + i b. Using the O(2)-invariance of the {N 1 , N 2 } parameter space in the degenerate limit, one can always rotate the vectors a and b, such that b is along the N 1 -axis, i.e. the N 2 -component of b vanishes, which in turn implies Im( h l2 ) = 0 (see Fig. 6).
For N N = 2, the R αβ -dependent terms in (2.18) are absent [cf. (A.1)]. Thus, in the degenerate limit m N α = m N β , the resummed Yukawa couplings given in (2.18) become Note that in the O(N N )-symmetric limit, any basis in the heavy-neutrino flavour space is a mass eigenbasis, and hence, the resummed Yukawa couplings can be defined consistently in any O(N N )-rotated basis. From (5.4), we find that in general, h l1,2 = 0 for h l1,2 ∈ C. However, if both h l1 and h l2 are real, i.e. b = 0, we can rotate a to align with either N 1 or N 2 direction in  Table 2 The resummed Yukawa couplings in the degeneracy limit, without including RG effects.
Similarly for three heavy-neutrino flavours (N N = 3), we can define a three-dimensional complex vector in the {N 1 , N 2 , N 3 } mass eigenbasis: h = ( h l1 , h l2 , h l3 ) = a + i b. In this case, using the O(3)-invariance of the parameter space in the degenerate limit, one can always rotate the vectors a and b in such a way that b points in the N 1 -direction and a lies on the (N 1 , N 2 )-plane, as shown in Fig. 7. Thus, the N 3 -components of both a and b identically vanish, i.e. Re( h l3 ) = 0 = Im( h l3 ). For the simple case with N L = 1, as considered above, the resummed Yukawa couplings flow to the O(3)-symmetric limit, i.e. h lα = 0 for h lα ∈ C. However, for h lα ∈ R, h lα is undetermined (0/0 form) in the degenerate limit. One can similarly work out the degeneracy limits for other values of N L ; the results of this analysis are summarized in Table 2.
In a realistic situation, the degeneracy of the heavy-neutrino parameter space in the O(N N )-symmetric limit M N → 0 will be broken by RG effects. Specifically, the RG evolution from a high scale μ X , at which the heavy neutrino masses are degenerate, i.e. m N α = m N β , to the scale m N induces a non-zero mass-splitting given by (5.2). For instance, for the case (N L , N N ) = (1, 3) discussed above, the inclusion of RG effects yields h lα = 0 for h lα ∈ C (α = 1, 2) and h l3 = 0. This is consistent with the fact that h l3 does not run in the mass eigenbasis, since h l3 (μ X ) can be rotated to zero, as shown in Fig. 7. As a consequence, the mass parameter m N 3 does not evolve under the action of the RG.
It is worth noting that the degeneracies encountered above are reminiscent of the singular degeneracies occurring in ordinary Quantum Mechanics, where one has to apply carefully degenerate time-independent perturbation theory in order to obtain meaningful results. For illustration, let us consider a three-state system whose time-evolution is governed by the following perturbed Hamiltonian (see [153], p. 348): where the unperturbed energy levels are E 1 , E 1 , E 2 (with E 2 > E 1 ), and a, b (E 2 − E 1 ) are treated as small perturbations. Since the perturbation matrix H is off-diagonal, the first-order perturbation does not change the energy eigenvalues, i.e. does not remove the degeneracy between the first two energy eigenstates. At second-order, applying non-degenerate perturbation theory would lead to an undetermined result for the corrections to the degenerate eigenvalues, i.e. (2) Instead, applying degenerate perturbation theory, in a suitable O(2)-rotated basis, leads to welldefined second-order energy shifts, thus lifting the degeneracy: (2) Note that the first energy eigenvalue will remain unperturbed to all orders. This is in close analogy with the (N L , N N ) = (1, 3) case discussed above, where m N 3 remains invariant, irrespective of the RG effects. Consequently, the resummed Yukawa couplings in the minimal RL model are finite and consistently flow to the O(N)-symmetric limit of the theory. The role of RG flow in lifting the degeneracy of heavy neutrino masses is akin to that of degenerate time-independent perturbation theory in ordinary Quantum Mechanics (e.g. electric field in a linear Stark effect). Just as one needs to choose carefully a basis in which to apply perturbation theory, one must define the resummed Yukawa couplings only in the mass eigenbasis, in which the RG effects consistently break the degeneracies of the heavy-neutrino parameter space.

A model of resonant τ -genesis
As an explicit example of the RL scenario, we consider an RL τ model with O(3) symmetry imposed on the heavy-neutrino sector at the GUT scale, μ X ∼ 2 × 10 16 GeV, which is explicitly broken to the U(1) L e +L μ × U(1) L τ subgroup of lepton-flavour symmetries by a neutrino Yukawa coupling matrix of the following form [76]: where δh vanishes in the flavour symmetric limit. In this symmetric limit, the light neutrinos remain massless to all orders in perturbation theory [164], while a and b are arbitrary complex parameters. In order to give masses to the light neutrinos, we consider the following form of δh as a minimal departure from the flavour-symmetric limit [76]: where | l |, κ 1,2 |a|, |b|, and γ 1,2 are arbitrary phases. To leading order in the symmetrybreaking parameters M N and δh, the tree-level light neutrino mass matrix is given by the seesaw formula In deriving this expression, we have assumed that The light neutrino mass matrix in (5.10) is diagonalized by the usual PMNS mixing matrix where the m ν i 's are the light neutrino mass eigenvalues. As we will see in Section 6.1, for the benchmark points considered therein, the non-unitarity of the 3 × 3 PMNS mixing matrix due to the light-heavy neutrino mixing [165] is very small. Hence, we can assume that U PMNS in ( with c ij ≡ cos θ ij , s ij ≡ sin θ ij . Assuming a particular mass hierarchy between the light neutrino masses m ν i 's and for given values of the CP-phases δ, ϕ 1,2 , we can fully reconstruct the light neutrino mass matrix using (5.12). Substituting (5.12) in (5.10), we can determine the following model parameters appearing in the Yukawa coupling matrix (5.8): In this way, the Yukawa coupling matrix (5.8) in the RL τ model can be completely fixed in terms of the heavy neutrino mass scale m N and the symmetry-breaking parameters κ 1,2 and γ 1,2 . Similar models can be constructed for RL e and RL μ scenarios [76].

Approximate analytic solutions
In this section, we will find some analytic solutions of the evolution equations in different regimes. In the first part, we qualitatively study the role of the heavy-neutrino coherences, and the generation of lepton number asymmetry via heavy-neutrino oscillations (see also Section 4.2). Here, we will obtain an approximate analytic solution for the rate equations with the εand ε -type CP asymmetries artificially switched off, and also neglecting the γ effects and charged-lepton off-diagonal number densities, for a simplified case with two heavy neutrinos.
Even though this is a rough approximation, this will allow us to estimate the relative magnitude of the different effects present in the statistical evolution of the system. In the second part, we will obtain a quantitatively accurate analytic solution for the case of diagonal heavy-neutrino number densities (hence no oscillations), but retaining the full off-diagonal number-density matrix for the charged leptons, thereby capturing the charged-lepton decoherence effects. In Section 6.2, we will show that the analytic solution presented here reproduces quite accurately the exact numerical solution of the full charged-lepton rate equation (4.84) in the attractor limit.

Qualitative estimate of the asymmetry via oscillations
Let us find the approximate solution to the simplified set of Eqs. (4.29), (4.32) and (4.33), with the tree-level Yukawa couplings, instead of the resummed ones, where the εand ε -type CP violating sources have been artificially switched off. For simplicity, we take the N N = 2 case and neglect the charged-lepton off-diagonal number densities. Promoting the vector η N in Section 2.2 to a matrix in the heavy-neutrino flavour space: and combining (4.29) and (4.32), we find where the matrix K N is defined as the tree-level generalization of K α 's appearing in (2.23): In the strong washout regime, i.e. for [ K N ] αβ 1, the system evolves towards the attractor solution, obtained by setting the RHS of (5.16) to zero: where we have neglected η η η N compared to 1. From (5.18), it is clear that all the elements of η η η N will have the usual 1/z behaviour, as expected in the attractor limit [cf. (2.25)]. The exact numerical solution of the fully flavour-covariant rate equations (4.82)-(4.84) also exhibit this behaviour, as shown explicitly in Section 6.2 (see Fig. 8).
To compute the charged-lepton asymmetry we are interested in the value of [δ η N ] 12 = 2η N eq i Im([ η N ] 12 ) (see the discussion in Section 4.2). From (5.18), we get  (5.20) where the index l is not summed over and The attractor solution is obtained by setting the RHS of (5.20) to zero: , we see that the total lepton asymmetry due to heavy-neutrino oscillations around z ∼ 1 is of the same sign and of the same order in magnitude, as compared to that obtained from the standard ε-type CP asymmetry due to heavy-neutrino mixing. Even though some of the approximations leading to this result may not be realistic in certain cases, we expect (5.21) to be qualitatively correct, and in Section 6.2 we will numerically verify this (see Figs. [9][10][11] for the RL τ model under consideration.

Analytic results for the charged lepton decoherence effect
We will now obtain the attractor analytic solution for the charged-lepton asymmetry neglecting the heavy-neutrino off-diagonal number densities and performing a number of approximations valid for the RL τ scenario discussed in Section 5.2. Neglecting the heavy-neutrino coherences and the sub-dominant γ LΦ LcΦc − γ LΦ LΦ term, but retaining the full flavour structure for the charged leptons, the rate equation for the asymmetry (4.84) can be written in the would-be mass eigenbasis for charged-leptons as follows: (5.22) where the various effective K-factors are defined as ). The CP asymmetries in the charged-lepton flavour space can be defined as . (5.24) Notice that, in general, this is a tensor in the charged-lepton flavour space, even though here we are working in the would-be mass eigenbasis for charged leptons. In the 2 heavy-neutrino mixing case, it can be approximated as This expression is analogous to that of the ε-type CP-asymmetry ε lα (see Appendix A) in the quasi-degenerate heavy neutrino limit. More precisely, for m N m N α , we have the relation ε ll = α ε lα , where ε lα is given in (A.2).
Using the analytic solution for the diagonal heavy-neutrino evolution equation (2.25) [ η N ] αα 1/([ K N ] αα z), we find the attractor solution in the strong washout regime by setting the RHS of (5.22) to zero, thus obtaining In the RL τ model discussed in Section 5.2, the dominant contribution to the total lepton asymmetry comes from the τ -sector involving [δ η L ] kτ (with k = e, μ, τ ) for which the third column of (5.26) provides a closed set of equations. Imposing the detailed balance condition (4.50), the third column of (5.26) can be explicitly written as ⎛ ⎜ ⎝ We can safely neglect [ K dec ] ee and [ K dec ] μμ which are much smaller compared to Assuming that the imaginary part of [ K eff ] τ k [δ η L ] kτ is small compared to its real part, (5.29) has the form of a closed linear system of equations for [δ η L ] kτ : Notice that, in the limit [ K dec ] τ τ → ∞, we recover the solutions of the diagonal rate equations in the RL τ model, i.e.
where we have only taken the real part of the solution, since the imaginary parts in the third row of (5.29) were assumed to be small in this derivation.
In the next section, we will show that the approximate analytic solution (5.32) reproduces the exact numerical solution of the rate equations remarkably well for the case with diagonal heavy-neutrino number densities, and provides a fairly good estimate for the total lepton number asymmetry predicted by the fully flavour-covariant rate equations.

Numerical examples
In this section, we present some numerical results for the evolution of the lepton asymmetry governed by the flavour-covariant rate equations given by (4.82)-(4.84). For definiteness, we choose to work within a minimal RL τ model presented in Section 5.2. For illustration, we take a set of neutrino Yukawa couplings satisfying the neutrino oscillation data for a normal hierarchy of light neutrino masses with the lightest neutrino mass m ν 1 = 0. We use the best-fit values of the light neutrino oscillation parameters from a recent three-neutrino global analysis [166]: For definiteness, we choose the leptonic CP phases δ = 0, ϕ 1 = π and ϕ 2 = 0, and reconstruct the light neutrino mass matrix using the definition (5.12). From this, we can determine the parameters a, b, e,μ,τ of the RL τ model using the relations (5.14) for a given value of the heavy neutrino mass scale m N , after taking into account the mass splitting between the three heavy neutrinos due to the RG effects given by (5.2). Note that for a given light neutrino mass matrix M ν , the solutions for a and b obtained using (5.14) are unique up to a sign factor, and the sign discrepancy could be eliminated only if the sign of Re(a) were known. There is a similar sign freedom for e,μ,τ , but this is irrelevant since it applies to a whole column of the Yukawa coupling matrix δh [cf. (5.9)] and can be rotated away. Thus, the only free parameters we have in this model, apart from the neutrino-sector CP phases, are κ 1,2 , γ 1,2 , sign[Re(a)] and m N . Below we present some benchmark values for these free parameters.

Benchmark points
We choose three benchmark scenarios with the heavy neutrino mass scales m N = 120 GeV, 400 GeV and 5 TeV, in order to illustrate the flavourdynamics of the RL model in different temperature regimes. The other free model parameters are chosen such that they satisfy all the relevant experimental constraints, as discussed below. The values of the free model parameters, viz. m N , γ 1,2 , κ 1,2 for our benchmark scenarios are given in Table 3, along with the corresponding values of the derived model parameters, viz. a, b, e,μ,τ , with Re(a) > 0. The low-energy Table 3 The numerical values of the free (m N , γ 1,2 , κ 1,2 ) and derived parameters (a, b, e,μ,τ ), with Re(a) > 0, in the RL τ model for three chosen benchmark points.  Table 4, along with the current experimental limits.

LFV observables
The mixing between the light and N N heavy Majorana neutrinos induces LFV processes such as → γ [9,[167][168][169][170][171][172], → 1¯ 2 [171,173,174] and μ → e conversion in nuclei [175][176][177][178][179], through loops involving the heavy neutrinos. In general, the light-heavy neutrino mixing is parametrized in terms of an arbitrary 3 × N N matrix ξ [163], which depends on the Yukawa coupling matrix h and the heavy Majorana neutrino mass matrix M N . In the mass eigenbasis, and to leading order in ξ , the mixing is given by 13 which governs the rare LFV decay rates, as discussed below. The branching ratio for the μ → eγ process is given by [170] BR where m μ and Γ μ are respectively the mass and width of the muon, s w ≡ sin θ w is the weak mixing parameter, and α w ≡ g 2 w /(4π) is the weak coupling strength, all evaluated at the weak scale M Z . The form factor G μe γ is defined in Appendix D. The other kinematically allowed rare decay modes of this type, namely, τ → μγ and τ → eγ , can be defined similar to (6.3), in terms of the mass and width of the τ -lepton. The branching ratio predictions of these rare decay modes for our chosen benchmark points are given in Table 4. For comparison, we also give the current experimental upper limits [1,184]. The MEG limit on μ → eγ branching ratio [184] is the most stringent one, and the model prediction for BP2 is within reach of the upgraded MEG sensitivity [185].
The branching ratio for the μ − → e − e + e − process is given by [174]  where the various form factors are defined in Appendix D. The predictions for BR(μ → eee) for the three benchmark points chosen here are shown in Table 4 and, for comparison, we have also shown the current experimental upper limit [1]. We find that the model predictions are well within the current limit. One can similarly define the LFV decay rates involving the τ -lepton [174]; however, the numerical values for these rates turn out to be several orders of magnitude smaller than the current experimental limits and hence we do not show them in Table 4. The μ → e conversion rate in an atomic nucleus A Z X is given by [179] where e = g w s w is the magnitude of the electron charge, Γ capt is the nuclear capture rate and V (p),(n) , D are various nuclear form factors, whose numerical values for some typical nuclei of interest are given in Appendix D [cf.   Table 4 for certain isotopes of titanium, gold and lead nuclei, along with their experimental upper limits from SINDRUM-II [186][187][188]. It is worth mentioning here that the next generation experiments, such as COMET [189] and Mu2e [190] have planned sensitivities around 10 −16 , which could easily test the first two benchmark points. The distant future proposal PRISM/PRIME [191] could probe μ → e conversion rates down to 10 −18 , thus testing the third benchmark point as well.
Apart from these LFV observables, a non-zero light-heavy neutrino mixing also leads to a non-unitary PMNS mixing matrix, which can be parametrized as [163,180,181] where U PMNS is the unitary matrix given by (5.13), and the non-unitarity effects are captured by the Hermitian matrix Ω, which is a function of the light-heavy neutrino mixing parameter ξ [cf. (6.2)]. To leading order in ξ 2 , the non-unitarity parameters are given by Ω α B * α B α . For the benchmark points given in Table 3, the predictions for |Ω| eμ are shown in Table 4, along with the current experimental limit at 90% CL from a global fit of neutrino oscillation data, electroweak decays, universality tests and rare charged-lepton decays [165]. The predictions for other elements of Ω are much below the current experimental  [195] limits, and hence, are not shown here. Note that an upgraded MEG limit on BR(μ → eγ ) could reach a sensitivity of |Ω| eμ < 10 −6 which includes the first two benchmark points in Table 4.
Since the non-unitarity parameter Ω is very small for all the benchmark points chosen here, we use U PMNS [cf. (5.13)] instead of U PMNS [cf. (6.7)] as the diagonalizing matrix in (5.12) for our phenomenological purposes.

LNV Observables
The Majorana nature of the light and heavy neutrinos in the type-I seesaw models violate lepton number by two units, which can manifest in the neutrinoless double beta decay (0νββ) process at low-energy (for a review, see e.g. [192]). In the minimal seesaw model, the 0νββ process gets contributions from diagrams mediated by both light and heavy Majorana neutrinos, and the corresponding half-life is given by where G 0ν 01 is the decay phase space factor, M 0ν ν,N 's are the nuclear matrix elements (NMEs) for 0νββ mediated by light and heavy neutrinos respectively, and the dimensionless parameters A ν,N are defined as where m e and m p are the electron and proton masses, respectively. For all the benchmark points given in Table 3, the heavy neutrino contribution A N to the half-life (6.8) turns out to be negligible compared to the light neutrino contribution A ν , and hence, we can ignore the second term on the RHS of (6.8), to rewrite it in the canonical form 1 T 0ν where m ≡ | i (U PMNS ) 2 ei m ν i | is known as the effective neutrino mass. For a normal hierarchy of neutrino masses with m ν 1 = 0, and using the three-neutrino oscillation parame-ter values given in (6.1), we obtain m = 3.8 meV. For comparison, we note that the current 90% CL experimental upper limits are m < (0.3-0.9) eV from the NEMO-3 limit on T 0ν 1/2 ( 100 Mo) [193], < (0.2-0.4) eV from the GERDA+Heidelberg-Moscow+IGEX combined limit on T 0ν 1/2 ( 76 Ge) [194], and < (0.12-0.25) eV from the KamLAND-Zen+EXO-200 combined limit on T 0ν 1/2 ( 136 Xe) [195], where the range of limits is due to the NME uncertainties involved. 14 Apart from the low-energy observables discussed above, the Majorana nature of the heavy neutrinos as well as their mixing with the light neutrinos could manifest simultaneously via their 'smoking gun' signature of same-sign dilepton plus two jets with no missing energy [197] at the LHC [164,[198][199][200][201][202][203]. Note that even for quasi-degenerate heavy Majorana neutrinos (as occurs in the RL models discussed in Section 5), the LNV signal can be sizable when the mass splitting m N ≡ |m N α − m N β | is comparable to the average width Γ N ≡ (Γ N α + Γ N β )/2 [204]. Within the minimal seesaw framework, both CMS [205] and ATLAS [206] experiments have derived limits on the mixing parameters |B lα | 2 (for l = e, μ) between 0.01-0.1 for m N = 100-300 GeV. Including infrared enhancement effects due to t -channel processes involving photons, these limits can be improved (by at least a factor of 5) and extended to higher heavy neutrino masses at √ s = 14 TeV LHC [82]. The improved limits will encompass the first two benchmark points in Table 3.

Results for lepton flavour asymmetries
Using the parameter values given in Table 3, we numerically solve the flavour-covariant rate equations (4.82)-(4.84) for the evolution of the charged-lepton and heavy-neutrino number densities. For definiteness, in this section we work in the basis in which the heavy-neutrino mass matrix as well as the charged-lepton Yukawa coupling matrix are diagonal. First, we discuss the results for the heavy-neutrino number densities, as shown in Fig. 8 for the three benchmark points given in Table 3. Here we have chosen the initial conditions with zero lepton asymmetry, i.e. δη L in = 0, and the heavy neutrinos in thermal equilibrium, i.e. η N in = η N eq 1. As we will see below, other choices of initial conditions lead to similar results. The vertical dotted line indicates the critical value z c = m N /T c , where T c is the critical temperature [cf. (2.39)]. The number densities are shown in terms of the deviation from their equilibrium values: η N αβ = η N αβ /η N eq − δ αβ [cf. (5.15)]. The evolution of the diagonal elements are shown as solid lines and that of the off-diagonal elements as dashed lines. As discussed in Section 5.3 [cf. (5.18)], both diagonal and off-diagonal heavy-neutrino number densities rapidly follow the attractor solution. Numerically, the value of η N 11 is several orders of magnitude larger than the other elements, whereas the values of η N 22 and η N 33 overlap for all the benchmark points. In all three cases, unlike η N 23 , the off-diagonal elements η N 12 and η N 13 are larger than the diagonal elements η N 22 , η N 33 . Therefore, the effect of the off-diagonal contributions of η N αβ to the lepton asymmetry cannot be neglected, as we will illustrate below.
Our results for the asymmetries in the SM lepton-doublet sector for the three benchmark points given in Table 3   (i) Zero lepton asymmetry δη L in = 0, heavy neutrinos in thermal equilibrium η N in = η N eq 1 (thick black line); (ii) Zero lepton asymmetry δη L in = 0, heavy neutrinos strongly out-of-equilibrium η N in = 0 (thick grey line); Fig. 9. Lepton flavour asymmetries as predicted by the BP1 parameters given in Table 3. The top panel shows the comparison between the total asymmetry obtained using the fully flavour-covariant formalism (thick solid lines, with different initial conditions) with those obtained using the flavour-diagonal formalism (dashed lines). Also shown (thin solid line) is the analytic result discussed in Section 5.3. The bottom panel shows the diagonal (solid lines) and off-diagonal (dashed lines) elements of the total lepton number asymmetry matrix in the fully flavour-covariant formalism. For details, see text. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.) (iii) Extremely large lepton asymmetry with opposite sign compared to the observed one δη L in = 1; heavy neutrinos strongly out-of-equilibrium η N in = 0 (thick yellow line).
It is clear that the final lepton asymmetry δη L (z 1) is independent of the initial conditions, which is a general consequence of the RL mechanism in the strong washout regime. Even starting with extremely large initial lepton asymmetry and/or with the wrong sign, this primordial asymmetry is rapidly washed out and the final asymmetry, with the right sign as required by (2.43), is set by the RL mechanism itself for z ∼ 1. Fig. 10. Lepton flavour asymmetries as predicted by the BP2 minimal RL τ model parameters given in Table 3. The labels are the same as in Fig. 9. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.) We may now compare the results of our fully flavour-covariant rate equations with their diagonal [29,33] and partially flavour-dependent limits. To this end, we show in the top panels of Figs. 9-11 limiting cases, where either the heavy-neutrino number density (red dashed line) or the charged-lepton number density (green dash-dotted line) or both (blue dotted line) are treated as diagonal in flavour space. For these cases, we have chosen the initial conditions δη L in = 0 and η N in = η N eq 1 for concreteness. As mentioned above, the final lepton asymmetry is insensitive to the initial conditions in each case. We also show the analytic solution, discussed in Section 5.3, for the case with diagonal heavy-neutrino number densities, which agrees well with the asymptotic limit of the corresponding exact numerical solution.
From Figs. 9-11, we see that the final asymmetry including all flavour effects is significantly enhanced in the fully flavour-covariant formalism, as compared to the predictions from the par- Fig. 11. Lepton flavour asymmetries as predicted by the BP3 minimal RL τ model parameters given in Table 3. The labels are the same as in Fig. 9. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.) tially flavour-dependent limits. The enhanced CP asymmetry in the flavour-covariant formalism can be understood as arising predominantly from two physically-distinct phenomena: (i) Coherent oscillations between different heavy-neutrino flavours, which create an O(h 4 ) asymmetry in the charged-lepton sector in the RL scenario, as discussed in Sections 4.2 and 5.3. This leads to an enhancement by a factor of two over the flavour-diagonal limit in all the benchmark points shown in Figs. 9-11. (ii) The evolution of flavour coherences in the charged-lepton sector, which are generated through the heavy-neutrino Yukawa couplings and destroyed through the charged-lepton Yukawa couplings, as discussed in Sections 4.3 and 5.3.
The latter charged-lepton decoherence effects give rise to the distinctive 'plateau' in the offdiagonal τ e and τ μ elements of the lepton asymmetry matrix δη L in the RL τ model. This can be seen from the bottom panels of Figs. 9-11, which show the individual charged-lepton flavour contributions to the total asymmetry. Since we are considering an RL τ model, the asymmetry produced around z = 1 is dominantly in the τ -flavour (δη L τ τ ), which has relatively smaller couplings to the heavy neutrinos [cf. (5.8)], and hence, a smaller washout factor. On the other hand, the asymmetries generated in electron (δη L ee ) and muon (δη L μμ ) flavours, with relatively larger couplings to the heavy neutrinos, are suppressed due to larger washout rates. In addition to the lepton asymmetry in the diagonal element δη L τ τ , the coherence effect in the charged-lepton sector generates an extra asymmetry in the off-diagonal number densities involving the τ -flavour, i.e. δη L τ e and δη L τ μ . These could be large compared to those involving other flavours, i.e. δη L ee , δη L eμ and δη L μμ , depending on the values of the input parameters. This effect is more prominent around z = 1, since with increasing z values, the off-diagonal lepton-flavour coherences decay, leading to a complete decoherence of the system to the charged-lepton Yukawa eigenbasis. This explains the distinctive plateau at intermediate z values in the evolution of δη τ e and δη τ μ in Figs. 9-11 (bottom panels). This additional source contributes to the total lepton asymmetry δη L , which can exhibit a similar feature, depending on the model parameters. In the case where the sphalerons freeze out within the 'plateau' region, which occurs for 200 GeV m N 2 TeV, an additional enhancement of a factor ∼ 5 in the final asymmetry is observed, as can be seen for BP2 in Fig. 10 (top panel). For BP3 with a much higher heavy-neutrino mass scale, the coherence effect is already subdued, and the system has completely decohered to the charged-lepton Yukawa eigenbasis, well above the critical temperature, thus giving no additional enhancement, as shown in Fig. 11 (top panel). On the other hand, for BP1, the neutrino Yukawa couplings are smaller than those in BP2 and BP3, and hence, the coherence effects are not pronounced in the total asymmetry, as can be seen from Fig. 9 (top panel).
The impact of the additional enhancements, discussed above, is exemplified in BP2 (see Fig. 10), where the total lepton asymmetry obtained in the fully flavour-covariant formalism is above the observed value, whereas the predictions obtained in various partially flavour-dependent limits all fall below the observed asymmetry. Note that the total lepton asymmetries predicted for all of the benchmark points exceed the observed asymmetry. Nevertheless, due to the freedom in the choice of the CP phases γ 1,2 , these benchmark points represent viable choices of model parameters for successful leptogenesis.
Before concluding this section, we note that there is no decoherence effect in the heavyneutrino sector, and hence, both diagonal and off-diagonal heavy-neutrino number densities decay coherently, as shown in Fig. 8. This is due to the fact that the evolution of the heavy-neutrino number densities are entirely governed by the heavy-neutrino Yukawa couplings [cf. (4.82) and (4.83)], ignoring sub-dominant collision terms, such as L = 1 scattering processes. In the heavy-neutrino mass eigenbasis, the corresponding decay rates are not diagonal, which leads to the occurrence of coherences in the heavy-neutrino sector.

Conclusions
We have presented a fully flavour-covariant formalism for transport phenomena, in which we have derived Markovian master equations describing the time-evolution of particle number densities in a quantum-statistical ensemble with arbitrary flavour degrees of freedom. In particular, we have obtained a flavour-covariant generalization of the semi-classical flavour-diagonal Boltzmann equations.
In order to explicitly demonstrate the importance of the effects captured only in the flavourcovariant formalism, we have discussed a particular application to the phenomenon of resonant leptogenesis (RL). It is known that the RL scenario offers a unique opportunity for testing the connection between the origin of neutrino mass and matter-antimatter asymmetry by the ongoing LHC experiments as well as by various low-energy experiments probing lepton flavour and number violation. For this reason, it is essential to capture all the flavour effects due to the heavy neutrinos as well as SM leptons in a consistent manner, in order to obtain a more accurate prediction for the baryon asymmetry in this scenario. As we have shown in this paper, including all flavour off-diagonal effects could enhance the predicted lepton asymmetry as much as one order of magnitude in certain RL models, as compared to predictions obtained from partially flavourdependent treatments. Thus, our flavour-covariant formalism allows us to access an enlarged parameter space of the RL models, which could be tested in ongoing and planned experiments at both the high energy and intensity frontiers.
The main new results of our fully flavour-covariant formalism for RL scenarios, as contained in the final rate equations (4.82)-(4.84), are the following: (i) The appearance of new rank-4 tensors in flavour space in transport equations (see Section 3.4). These are necessary to describe the time-evolution of the number density matrices for leptons and heavy neutrinos in a flavour-covariant manner. One can extend this formalism, by introducing even higher rank rate tensors, to describe sub-dominant processes, such as LN ↔ Le R , involving more flavour degrees of freedom. The existence of the tensorial structure in the rate equations is firmly supported by an explicit calculation of the transition matrix elements, by virtue of a generalization of the optical theorem (see Appendix C). To further elucidate the consistency of our treatment, we develop a flavour-covariant generalization of the helicity amplitude technique, applied to spinorial fields in the presence of time-dependent and spatially-inhomogeneous backgrounds (see Appendix B). (ii) A systematic treatment of two intrinsically quantum effects, i.e. oscillations between different heavy neutrino flavours (see Section 4.2) and quantum decoherence between the charged-lepton flavours (see Section 4.3). Numerical studies for a particular RL τ model reveal that these flavour off-diagonal effects could enhance the total lepton asymmetry by up to one order of magnitude, as compared to the flavour-diagonal case (see Figs. 9-11). (iii) The approximate analytic solutions (see Section 5.3) to the fully flavour-covariant transport equations, which capture the two relevant flavour effects discussed above. Taking this into account in the strong washout regime, the total lepton asymmetry at z 1 may be estimated by the sum of the contributions from flavour mixing, oscillation and decoherence effects: Aside from these main results specific to the flavour-covariant formalism, we have also given a geometric and physical understanding of the degeneracy limit in the heavy neutrino parameter space (see Section 5.1) for the minimal RL scenario in which the quasi-degeneracy of the heavy neutrino masses at low scale can be naturally explained as a small deviation from the O(N N )-symmetric limit at some high scale through RG effects. We also point out that the role of RG effects in lifting the degeneracies encountered here is reminiscent of the role of timeindependent perturbations in the degenerate perturbation theory of ordinary Quantum Mechanics.
We also comment on the various existing forms of the self-energy regulator used to calculate the ε-type CP-asymmetry and make a comparative study in a simple toy model, in order to demonstrate their behaviour in certain lepton number conserving limits (see Appendix A). We find that only the regulator given by (A.3) gives a valid and well-defined CP-asymmetry in the entire parameter space possessing the correct L-conserving limit, whereas the other regulators are not well-defined in certain regions of the parameter space.
In conclusion, our flavour-covariant formalism provides a complete and unified description of transport phenomena in RL models, capturing three relevant physical phenomena: (i) the resonant mixing between the heavy neutrino states, (ii) coherent oscillations between different heavy neutrino flavours, and (iii) quantum decoherence effects in the charged-lepton sector. The formalism developed here is rather general and may also find applications in various other transport phenomena involving flavour effects.
Using (A.1) in (2.22) and neglecting higher-order Yukawa couplings of O(h 4 lα ) at the amplitude level, we obtain the following expression for the ε-type CP-asymmetry 15 : where α, β = 1, 2 (α = β), and the self-energy regulator is given by [28,29] Note that h and Γ (0) N β are the tree-level Yukawa couplings and decay width, respectively. In the degenerate heavy neutrino mass limit, m N ≡ (m N 1 − m N 2 ) → 0, the would-be singular behaviour of the CP-asymmetry is regularized by the absorptive term (m N α Γ (0) N β ) 2 in the denominator on the RHS of (A.3).
Based on the simplified expression (A.2), the following two necessary conditions for a resonant enhancement of the leptonic CP-asymmetry may be derived [28]: The generic feature of these resonant conditions remains valid even in the presence of flavour effects [33]. The condition (i) is exactly met when the unitarity limit on the resummed heavyneutrino propagator gets saturated [27], i.e. when the total CP-asymmetry ε l = α ε lα takes its maximum possible value of unity. Note that the limit ε l ≤ 1, similar to the Lee-Wolfenstein bound for the K 0 K 0 -system [207], gets saturated when the regulator in (A.2) takes its maximum possible value of |f max reg | = 1/2. Since the exact location of the pole of the propagator determines the maximum enhancement of the CP asymmetry as well as its behaviour in some limiting cases, it is worth commenting on other analytic forms of the regulator for the resonant part of the lepton asymmetries in (A.2), existing in the literature [118,132,133]. In particular, their results differ by the way in which the singularity m N → 0 in (A.2) is regulated once the heavy neutrino mixing effects are taken into account. For instance, using a perturbative quantum-mechanical approach, [133] obtained a regulator of the form which has a pathological behaviour for certain flavour-dependent RL scenarios, for which Re(A αβ ) = 0, but Re(A l αβ ) = 0. In this limit, the unitarity upper bound ε l ≤ 1 is violated in the degenerate heavy-neutrino mass limit m N → 0 and the individual CP asymmetries ε l become singular.
Following a modified version of the quantum field-theoretic approach introduced in [28], a different regulator was obtained by [132,138]: This regulator diverges in the doubly degenerate limit m N → 0 and Γ N ≡ |Γ N 2 even though A l 11 = A l 22 for a given lepton flavour l. For instance, such a scenario can occur naturally in approximate L-conserving RL models [54,58]. In these cases, using the regulator (A.6) might lead to an overestimation of the leptonic CP asymmetries by several orders of magnitude [76], as we will explicitly demonstrate below for a toy model.
Finally, using an effective Kadanoff-Baym approach with a particular quasi-particle ansatz [118,120], derived an effective regulator of the form The predictions for the CP-asymmetry for the regulators (A.6) and (A.7) are comparable so long as the widths Γ N 1,2 are hierarchical. Although the regulator (A.7) does not have any pathological behaviour in the doubly degenerate limit m N → 0, Γ N → 0, it predicts a CP-asymmetry 4 times smaller than that predicted by the regulator f reg in (A.3). As a result, the CP-asymmetry ε ≤ 1/2 never gets saturated to unity, unlike the general expectations based on unitarity arguments [27].  [54,208] where M, μ 1,2 are real parameters. The full neutrino mass matrix, in the basis {ν C L,l , N 1 , N C 2 }, is given by To first order in the lepton number breaking parameters Y l and μ 1,2 , the light neutrino mass matrix is given by Note that, although at tree-level the light neutrino mass matrix does not depend on μ 1 , it receives a one-loop contribution from the electroweak radiative corrections [164,180], which is directly proportional to μ 1 . This is taken into account by the effective μ-parameter in (A.10): μ eff = μ 2 + x N f (x N )μ 1 , where x N = M 2 /M 2 W and the loop function f (x N ) is given by [180] f (x N ) = α w 16π Therefore, in order to satisfy the neutrino mass constraints, we require the lepton-number breaking parameters to be small, i.e. |Y l | |Y l | and μ 1,2 M.
In the basis in which the heavy neutrino mass matrix is diagonal with real and positive eigenvalues, the Lagrangian (A.8) can be recast into the general form given by (3.1), with α = 1, 2 and the heavy neutrino masses M 1,2 M ∓ μ/2. Here, we have defined μe iφ = μ 1 + μ 2 e iθ . To leading order in μ/M and Y , the Yukawa couplings in this mass eigenbasis are given by [54]  where λ = sin θ(μ 1 μ 2 /μM). Using the tree-level Yukawa couplings given by (A.12) and the corresponding decay widths given by (2.17), we calculate the ε-type CP-asymmetry in a particular lepton flavour ε l = α ε lα , with ε lα given by (A.2), and compare the magnitudes of ε l obtained using the two regulators given by The singular behaviour of the regulator is (A.6) is further illustrated in Fig. A.2 where we plot the ratio ε l /κ l , with κ l defined by (2.31). Note that in the L-conserving limit, κ l vanishes by construction, since in this limit, there are no γ LΦ L c Φ c terms in (2.31). Now in the same L-conserving limit, the CP-asymmetry obtained using the regulator (A.3) also vanishes, as shown in Fig. A.1  (left panel). For the toy model under consideration, both ε l and κ l go to zero at the same rate, while keeping ε l /κ l constant, in the limit μ 1 → μ 2 , as shown in Fig. A.2 (left panel). Note that this does not mean the final lepton number asymmetry given by (2.33) is non-zero in the L-conserving limit, since the expression (2.33) is valid only in the strong washout regime. In the weak washout regime with small κ l , the final lepton number asymmetry will instead be proportional to ε l alone [cf. (2.32)], and therefore, vanishes as long as ε l → 0. However, for the regulator given by (A.6), due to the fact that ε l does not vanish in the μ 1,2 → 0 limit, the ratio ε l /κ l blows up, as shown in Fig. A.2 (right panel).
It was argued in [138] that, due to the limitations of the perturbative approach used to derive (A.6), the range of validity of the analytic expression (A.2) with the regulator structure given by (A.6) is restricted to the parameter space in which the degree of degeneracy of heavy neutrinos must be much larger than the pole expansion parameter, determined by the neutrino Yukawa couplings. In the toy model presented above, this validity condition can be written as from which it is clear that for any choice of θ = (2n + 1)π/2 (with n = 0, 1, 2, . . .), the condition (A.13) is always satisfied, as long as |Y l | 2 ≤ 1. Thus, the pathological behaviour of the regulator (A.6) cannot be avoided simply by imposing the validity condition (A.13).
On the other hand, for the analytic solution (A.7) obtained in the Kadanoff-Baym approach, the expansion for small Yukawa couplings requires an additional validity condition [118] Re  Comparing (A.14) and (A.18), we see that the condition (A.17) is not satisfied in the limit μ 1 → μ 2 and hence, the regulator (A.7) is also not valid in the pathological limit. These issues do not arise for the regulator (A.3), which is well-defined in the entire parameter space shown in Fig. A.1.
For comparison with the ε-type CP-asymmetry discussed above, let us also compute the ε -type CP-asymmetry due to one-loop vertex corrections. For the two heavy neutrino case, this is given by [32] is the Fukugita-Yanagida loop function [2]. The magnitude of the ε lα contribution is shown numerically in Fig. A.3 for our toy model with the same parameter values as chosen for Fig. A.1. We find that, as expected, the ε -type CP-asymmetry vanishes in the L-conserving limit. Moreover, the ε -part of the transition amplitude squared for the N 1 decay becomes equal and opposite in sign to that for the N 2 decay, and thus, these two contributions cancel each other to give a vanishing total CP-asymmetry ε l = α ε lα , which goes to zero faster than the individual contributions in the limit μ 1,2 → 0, similar to the ε-case with the regulator (A.3). We also note that the ε -type CP violation given by Fig. A.3 is smaller than the ε-type CP violation given by Fig. A.1 (left panel) for all values of μ 1,2 . This is consistent with the general expectation that the ε -type contribution can become comparable to the ε-type term only in the hierarchical limit m N 2 m N 1 , when it approaches the asymptotic relation l ε lα = (1/2) l ε lα [32].

Appendix B. Flavour covariant helicity amplitude formalism
In this appendix, we describe the pertinent details of the fully flavour-covariant quantization of spinorial fields in the presence of time-dependent and spatially inhomogeneous backgrounds. Thereby, we exemplify the consistency of our treatment of flavour mixing in the transport equa-tions derived in Sections 3 and 4. To this end, we begin by describing the flavour-covariant generalization of the helicity amplitude formalism [209,210], in the context of an N -flavour model of Dirac fermions. After highlighting the generalized discrete symmetry transformations of the Dirac helicity four-spinors, we illustrate the inter-dependence of the spinorial and flavour structure. 16 Subsequently, we derive the flavour-covariant propagators of our toy N -flavour model, thereby generalizing the non-homogeneous propagators described in [126]. In addition, we introduce the spatially inhomogeneous and time-dependent flavour-covariant statistical distribution functions relevant to a complete quantum field theoretic treatment of flavour-coherent and Gaussian statistical backgrounds. Finally, we highlight a consequence of this treatment that is anticipated to impact upon the flavour-dependent quasi-particle approximations currently employed in the literature in the application of the CTP formalism and resulting Kadanoff-Baym equations to transport phenomena. Specifically, we show how the time-translational invariance of flavour-covariant propagators is necessarily broken in the presence of flavour-coherent statistical backgrounds.

B.1. Flavour covariant spinor algebra
We begin by introducing the four-component Dirac spinor ψ k and its Dirac conjugate ψ k in the Weyl basis The U(N )-symmetric Dirac Lagrangian may be written in the following form: where the gamma matrices are defined in the Weyl basis: with σ μ ≡ (1, σ ) and σ μ ≡ (1, −σ ), σ i 's being the usual 2 × 2 Pauli matrices. The mass matrix m l k transforms as a rank-2 tensor under U(N ). We may rotate to the Dirac basis, independent of the flavour structure, by means of the orthogonal transformation where the gamma matrices in the Dirac representation are 16 Previously, the issue of spin coherence in quantum kinetic equations has been considered using a truncated gradient expansion of Wigner functions [211], which we do not follow here.
with 0 2 and 1 2 being the 2 × 2 null and identity matrices, respectively. The relevant flavour transformations for the Dirac fields and the mass matrix are where U k l ≡ (U k l ) * and U k l U k m = U l k U m k = δ l m . Varying the Lagrangian with respect to the fields, we obtain the flavour-covariant Dirac equations where, in the latter, the derivative acts to the left and ψ k (x) = [ψ k (x)] † γ 0 is the Dirac-conjugate spinor.
The Dirac field operators in the interaction picture may be written as follows: where s = ± is the helicity index, denoting the two helicity states with the unit spin vector n = ss aligned parallel and anti-parallel to the three momentum p, respectively, i.e. s = p/|p| = (sin θ cos φ, sin θ sin φ, cos θ). (B.10) The three momentum p is obtained by boosting from the rest frame along the direction specified by s. Notice that the four-component Dirac spinors u and v transform as rank-2 tensors under U(N ). In addition, we draw attention to the fact that b k and d † k (b k and d †,k ) transform under the same fundamental (anti-fundamental) representation of U(N ). The particle and anti-particle creation and annihilation operators satisfy the anti-commutation relations The creation and annihilation operators of mass dimension −1 in (B.11) are related to the corresponding ones of mass dimension with analogous expressions for the antiparticle creation and annihilation operators, obtained by the replacement b → d † and ω b → ω d of the arbitrary Bogoliubov angles. The matrix G lm is defined in (3.28).
In momentum space, the Dirac field operators take the forms The matrix product of two spinors is given by where T = iγ 1 γ 3 . All of these identities may readily be verified using those in (B.24) for the two component spinors. Notice that in the mass eigenbasis, where the G matrices are proportional to the identity, we recover the expected results.

B.3. Chiral field operators
We write the chiral projection operators by introducing a dummy index χ = ±, such that P + ≡ P R and P − ≡ P L , as appearing in (2.15).
With the chiral field definitions the Dirac Lagrangian (B.2) takes the form Notice that in the relativistic limit E m, the helicity and chirality states coincide and only the ξ χ spinors survive. Hereafter, we neglect the ϕ s u(v) -dependent phases for notational convenience.

B.4. Spinor traces
We may now define the following spinor trace involving the objects defined in (B.53): Tr P p, s; p , s lm k n P χ P q, r; q , r l m k Summing over the helicities s and r, we find s,r Tr P (p, s; p, s) ml k m P χ P (q, r; q, r) m l k (B.77) The statistical distribution functions f and g in (B.80) satisfy the following properties: f ss p, p , t k l = f s s p , p, t l k * = f ss p, p , t k l C * , (B.81) g ss p, p , t k l = ḡ ss p, p , t k l C * .
We assume a Gaussian density operator, so that we must specify only the bilinear EEVs of operators. Notice that in this case only the four combinations in (B.80) are permitted by the associated spinorial structure. In the homogeneous limit, we have the correspondence The positive and negative-frequency fermionic Wightman propagators are defined as By evaluating the EEV of field operators directly, we obtain the explicit forms Rotating to the mass eigenbasis, we obtain If the system is out-of-equilibrium, flavour coherences cannot be neglected. As a consequence, the statistical distribution function [f ] k l will in general have non-zero off-diagonal entries. In this case, we see from (B.91) that the Wightman propagator depends explicitly on two zeroth component momenta p 0 and p 0 , which need not be equal. In this case, the time-translational invariance of the flavour-covariant propagators is necessarily broken.
Assuming that f is diagonal in flavour space and helicities in thermodynamic equilibrium, we have f hom, ss p 0 , p 0 , p, t k l →f eq, k p 0 , p 0 , p, t δ l k δ ss , is the Fermi-Dirac particle distribution, and f F (E) is the corresponding anti-particle distribution with μ → −μ. In this equilibrium limit (B.92), the Wightman propagators (B.90) reduce to iS eq,≷ p, p ,t k l = 2π |2p 0 | 1/2 δ p 2 where, following (B.53), we have defined Shifting the boundary time t i → 0, so that t →t −t i = t in (B.91), we see that the degree of violation of time-translational invariance is controlled by the phase e i(p 0 −p 0 )t . Working in the mass eigenbasis, we find that the violation of translational invariance is maximal when where we have assumed |p| m 2 kl and m 2 kl = m 2 k − m 2 l is the mass splitting. This violation of time-translational invariance is expected in the presence of flavour coherences, since the non-equilibrium propagators describe correlations of coherent superpositions of states, and these superpositions are not eigenstates of the Hamiltonian. In the CTP formalism, these coherence effects will be captured in the memory integrals occurring in the collision terms of the resulting master equations. We anticipate that this violation of time-translational invariance due to flavour mixing may have significant impact on the quasi-particle resummations currently employed in the applications of Kadanoff-Baym equations to such phenomena [118]. A detailed discussion may be given elsewhere.
Before concluding Appendix B, let us consider the Schwinger-Dyson equation of the fermion propagators. Working in the double momentum representation discussed above, this reads as follows: where [S 0 ab ] k l and [S ab ] k l (with a, b = 1, 2 being the CTP indices) are respectively the free and resummed 2 × 2 matrix CTP propagators in the doublet notation employed by [213,214] Due to the violation of time-translational invariance, the resulting Feynman-Dyson series will contain an infinite nesting of momentum integrals, which will not collapse to the usual algebraic equation of resummation, i.e. for the time-ordered (1, 1) component of the equilibrium CTP propagator for a single flavour: where Σ(p 0 , p) is the time-ordered self-energy. This result can readily be verified by rotating to the CTP eigenbasis (or so-called Feynman basis) [216] (see also [217]). In (B.99), the function Σ (p 0 , p) = Re[Σ ret (p 0 , p)] +iε(p 0 ) Im[Σ ret (p 0 , p)] and its complex conjugate Σ * (p 0 , p), written in terms of the retarded self-energy Σ ret (p 0 , p), are the eigenvalues of the CTP self-energy matrix, having no physical significance at finite temperature (see e.g. [126]). The self-energy function Σ (p 0 , p) (and likewise the time-ordered self-energy Σ(p 0 , p)) permits the following decomposition in terms of the Dirac gamma matrices γ μ = (γ 0 , γ ) [cf. (2.15)]: where we emphasize the violation of Lorentz-covariance due to thermal effects. In (B.100), the form factors Σ S L,R correspond to the scalar components of the self-energy and Σ 0,V L,R the chiral components. In the zero-temperature limit, Σ coincides with the time-ordered self-energy Σ and (B.99) reduces to the expected form in which Lorentz-covariance is restored, since the form factors are functions only of p 2 and . Maintaining the thermodynamic equilibrium in (B.92) and rotating to the mass eigenbasis simultaneously, we can reduce the Wightman propagators (B.94) to in which time-translational invariance is restored in the free propagators. If one also makes the Markovian approximation that energy is conserved in the interaction vertices, i.e. the interactions take place over an infinite time domain, time-translational invariance is restored globally and the Schwinger-Dyson equation (B.97) reduces to the relatively simpler form for which an exact matrix inversion both in the CTP and flavour structure is possible.
In summary, we have presented a fully flavour covariant formalism of helicity amplitudes. In addition, we have shown that the resummation of self-energy corrections in the case of flavour mixing may be performed in closed algebraic form only in a flavour-diagonal thermodynamic equilibrium. Thus, we justify the resummation scheme employed in the derivation of the Markovian master equations of Section 4 in which the resummed Yukawa couplings are obtained at zero-temperature, whilst the aforementioned flavour coherence effects discussed above have been included at the level of the quantum statistics.

Appendix C. Generalized optical theorem
In this appendix, we justify the tensorial flavour structure of the rates introduced in Sections 3 and 4 by means of an explicit calculation of transition matrix elements. To this end, we derive a background-dependent analogue of the optical theorem that is able to account for off-diagonal flavour coherences.
We begin by writing the scattering operator S in terms of the transition operator T as S = 1 + iT . Subsequently, using the unitarity of the scattering operator S † S = SS † = 1, one can easily show that In the usual derivation of the optical theorem, we would now proceed by multiplying (C.1) from the right and left by a given initial Fock state and insert a complete set of final states between the transition operators on the RHS of (C.1). For our purposes, the completeness of the Fock space can be written in the short-hand notation where the sum over A runs over all possible multi-particle states, e.g. Higgs, heavy-neutrinos and charged-leptons in our case, and also contains the helicity summations, isospin traces and momentum integrals. Writing explicitly, (C.2) has the form g=0 s g ,I g p g p g , s g , I g , L l g p g , s g , I g ,  (2) isospin. The heavy-neutrino states |k, r, N α are isospin singlets. In (C.3), the wedge product denotes the anti-symmetrized tensor product, e.g.
The zeroth term in each of the summations over g, h, i and j in (C.3) is understood to be the outer product of vacuum states, e.g. In order to incorporate these off-diagonal flavour coherences, we must take into account the contribution of the background ensemble, which we assume to be described by a factorizable density operator ρ = ρ N ⊗ ρ L ⊗ ρ Φ , where coherences between different particle species have been neglected. Note that the density operator itself is also a singlet under U(N N ) ⊗ U(N L ). Throughout this section, we suppress the time-dependence of the density operator ρ ≡ ρ(t) for notational convenience.
In the presence of a background ensemble, we proceed by taking the EEV of the unitarity relation in (C.1), yielding 2 Im T = 2 Tr[ρ Im T ] = Tr ρT † T . (C.7) Inserting the completeness of the Fock space in (C.2) and using the cyclicity of the trace operation, we obtain The sets of states A and B consist of all other possible multi-particle heavy-neutrino, lepton, anti-lepton and Higgs spectator states. 18 Using the orthonormality of the Fock states, we have F r ,β k F α r (k) = B ∧ k , r , N β k, r, N α ∧ A = (2π) 3 δ β α δ rr δ (3) k − k δ A B . (C.10) Consequently, we obtain from (C.8) the following equality: 2 Im T a β k , r a α (k, r) = 2 The creation and annihilation operators appearing in this appendix are understood to be asymptotic 'in' operators, see (3.70). By inspection, we see that the LHS of (C.11) defines the total in-medium heavy-neutrino production rate Γ rr {X} → N ; k, k α β ≡ 2 Im T a β k , r a α (k, r) , (C. 12) in which {X} represents the set of all possible initial states and at least one heavy-neutrino appears in the final state. All internal degrees of freedom in the set of initial states {X} are summed over, except for isospin, i.e. the general rate is a tensor in isospin space. We may convince ourselves that (C.12) necessarily gives the production rate (and not the decay rate), since in the zero-temperature limit, where ρ = |0 0|, the EEV 18 These spectator states should not be confused with the SM spectator processes which could enhance the washout of the lepton asymmetry in thermal leptogenesis [218].
on the RHS of (C.12) is identically zero. Similarly, the total in-medium LΦ production rate is given by Γ s 1 s 1 {X} → LΦ; p 1 , p 1 , q 1 , q 1 k l ≡ 2 Im T b l p 1 , s 1 b k (p 1 , s 1 )c † q 1 c(q 1 ) , (C.14) in which at least one lepton and one Higgs appears in the final state. In order to study specific contributions to the total heavy-neutrino and LΦ production rates, defined in (C.12) and (C.14), we must isolate particular initial states and truncate the transition operator to a given order in perturbation theory. Truncating the transition operator to first order in the interaction Hamiltonian density, T → T 0 = x H int (x), we find from (C.12) that only states in {X} that contain at least one lepton and one Higgs will contribute to the heavy-neutrino production rate. Isolating the initial states |p , s , L l ; q , Φ † ∧ |C ⊂ |C and |p, s, L k ; q, Φ † ∧ |D ⊂ |D , we obtain from (C.11) and (C.12) the tree-level heavy-neutrino production rate In (C.15), the spectator states A do not contribute to the transition matrix elements k, r, N α |T 0 |p , s , L l ; q , Φ † and p, s, L k ; q, Φ † | T 0 † |k , r , N β at first order in perturbation theory. The ρ-dependent term in (C.15) is A A ∧ p , s , L l ; q , Φ † ρ p, s, L k ; q, Φ † ∧ A = Tr ρb k (p, s)b l p , s c † (q)c q = n L s s p , p l k n Φ q , q , (C. 16) where we have taken the leptonic and Higgs ensembles to be spatially inhomogeneous in general, depending explicitly on two three-momenta. The time-dependence of the distribution functions n X (p, p ) is assumed implicitly. In the spatially-homogeneous limit, the lepton and Higgs distribution functions satisfy the correspondence n L ss p, p k l → n L s (p) k l (2π) 3 δ ss δ (3) p − p , (C.17) The general inhomogeneous distribution functions n X (p, p ) are related to the number densities defined in Section 3.1 by a Wigner transformation [219] and integration over all coordinate space: The Wick contraction of field operators may be performed using the following set of flavourcovariant field-particle duality relations: where the creation and annihilation operators have the mass dimension −3/2 and satisfy the commutator algebra given in Section 3.1. We emphasize that the operators appearing in this appendix are understood to be asymptotic 'in' operators.
In the case of the tree-level heavy-neutrino production rate, we obtain the Hermitian conjugate pair of transition matrix elements As in Section 4.1, we assume that the charged-leptons are massless, and hence, only one of their helicity states are populated, which for concreteness, we choose to be s = − (+) for leptons (anti-leptons) in (C.42). Summing over the heavy-neutrino helicities, tracing over lepton and Higgs isospins, and performing the momentum integrals, we derive the tree-level thermallyaveraged rates for the processes N → LΦ and LΦ → N , as follows: where we have divided through by the coordinate-space four-volume V 4 = (2π) 4 δ (4) (0) in order to remove the factor arising from the product of identical energy-momentum conserving delta functions from the two transition matrix elements in (C.42). After performing the summations over the heavy-neutrino helicities and the traces over lepton and Higgs isospin, we find from (C.50) the heavy-neutrino decay and inverse decay rates where we recall that g L = g Φ = 2 count the degenerate isospin degrees of freedom of the lepton and Higgs doublets. In addition, we have used the notation defined in (2.7) for the thermallyweighted phase-space integrals. We note that a relabeling of indices k ↔ l and α ↔ β has been performed in (C.52) relative to the T -transform of (C.51), that is (C.53) The decay rates in (C.51) and (C.52), derived using the generalized optical theorem are exactly the same as those given in (4.10) and (4.11).

C.2. Scatterings
We may write the homogeneous limit of the heavy-neutrino Feynman propagator in a general basis by rotating (C.40) to the mass eigenbasis, i.e. Rotating the matrix elements and four-momentum delta functions to the heavy-neutrino mass eigenbasis, and working in the massless limit for the charged-leptons and Higgs, we obtain where the four-momentum k = p 1 + q 1 = p 2 + q 2 . Finally, we substitute for the resummed Yukawas and propagators and extract the resonant part of the amplitude by writing i F,N (k) α P R (/ k + m α )P L i F,N (k) res,α = P R / k i F,N (k) res,α . (C.62) After performing the helicity summations via appropriate Fierz rearrangement, we obtain the following set of L = 0 and L = 2 thermally-averaged scattering rates:    in which the indices for the T -transformed processes have been relabeled, as in (C.53). Tracing over lepton flavour indices, the scattering rates relevant to the heavy neutrino transport equations are, for instance,

Appendix D. Form factors for LFV decay rates
For completeness, we list below various form factors appearing in the expressions (6.3), (6.4) and (6.5), which follow from the results given in [174,178]: with the limiting values H Box (0, 0) = 4 and F Box (0, 0) = 1. The nuclear form factors D, V (p) , V (n) and the capture rate Γ capt appearing in (6.5) are summarized in Table D.1 for various nuclei relevant to μ → e conversion searches. Table D.1 The nuclear form factors and capture rates for various nuclei of interest. The numbers were taken from [179] (see also [220,221]