Reconciling resonant leptogenesis and baryogenesis via neutrino oscillations

Right-handed neutrinos offer an elegant solution to two well established phenomena beyond the Standard Model (SM) - masses and oscillations of neutrinos, as well as the baryon asymmetry of the Universe. It is also a minimalistic solution since it requires only singlet Majorana fermions to be added to the SM particle content. If these fermions are nearly degenerate, the mass scale of right-handed neutrinos can be very low and accessible by the present and planned experiments. There are at least two well studied mechanisms of the low-scale leptogenesis: baryogenesis via oscillations and resonant leptogenesis. These two mechanisms were often considered separate, but they can in fact be understood as two different regimes of one and the same mechanism, described by a unique set of quantum kinetic equations. In this work we show, using a unified description based on quantum kinetic equations, that the parameter space of these two regimes of low-scale leptogenesis significantly overlap. We present a comprehensive study of the parameter space of the low-scale leptogenesis with the mass scale ranging from $0.1$ GeV to $\sim 10^6$ GeV. The unified perspective of this work reveals the synergy between intensity and energy frontiers in the quest for heavy Majorana neutrinos.


I. INTRODUCTION
The origin of the light neutrino masses and the baryon asymmetry of the Universe remain some of the burning problems for physics beyond the Standard Model (SM). Adding heavy neutrinos to the SM offers an economical solution to both problems (throughout this work we will use the terms heavy neutrinos and HNLs-Heavy Neutral Leptons interchangeably).
The masses of active neutrinos are coming from the type-I seesaw mechanism [1][2][3][4][5][6], the baryon asymmetry of the Universe (BAU) is generated through combined action of anomalous processes with fermion number non-conservation [7] and lepton number and flavor violating reactions involving heavy neutrinos. The scale of the heavy neutrino masses however remains an open question. The seesaw mechanism on its own does not imply any specific scale for the heavy neutrino masses. One may wonder if the answer can come from leptogenesis. And indeed, the specific patterns of Majorana masses of HNLs single out different scales. For example, if the spectrum of HNLs is hierarchical, M I M J , a lower bound of 10 9 GeV was found for the mass of the lightest heavy neutrino [8]. The mechanism leading to BAU in this case is known as thermal leptogenesis. If the spectrum is nearly degenerate, i.e. there is a pair of HNLs such that ∆M M , the leptogenesis may take place for M as small as O(100) GeV or even few MeV [9]. The corresponding leptogenesis mechanisms are called resonant leptogenesis [10][11][12][13][14][15][16][17][18][19], and baryogenesis via oscillations [20,21].
In the recent years these last scenarios have received significant attention, from both the theoretical (see e.g. [9,) and experimental (see, e.g. ) perspectives. This interest is primary related to the potential testability of the model in the present and near future experiments. HNLs are a canonical example of feebly interacting particles (see, e.g. [95,96]), which could have avoided discovery not because they are heavy but because their interactions are very weak.
For the observed BAU to be generated, the three Sakharov conditions need to be met [97]: 1) Efficient baryon number violation.
2) Sizeable C and CP violation.
Heavy Majorana neutrinos can provide both CP violation via the complex Yukawa couplings, as well as a deviation from equilibrium. In leptogenesis scenarios the asymmetry is generated in the lepton sector (hence the name) and transferred to the baryon sector by non-perturbaitve sphaleron processes which violate baryon number conservation [7]. In the thermal leptogenesis the deviation from equilibrium takes place during the decays of HNLs, when the decay rate cannot catch up the expansion rate of the Universe. The asymmetry between decays into leptons and anti-leptons comes about because of the interference between tree and one loop processes like N → H. See, e.g. [98] for details. This decay asymmetry is further enhanced if the two HNLs are degenerate in mass, as it is the case in the resonant leptogenesis. 1 We will discuss this mechanism in greater detail in section IV A.
In the baryogenesis via neutrino oscillations, the deviation from equilibrium happens during freeze-in. As a result of the seesaw mechanism, the Yukawa couplings of light HNLs must be tiny. Thus the equilibration rate is much lower than the Hubble rate, so the HNLs remain out of equilibrium. The asymmetry generated in the processes such as scatterings, decays, and inverse decays of the HNLs is further enhanced by their oscillations.
Several studies of resonant leptogenesis suggested that the scale of heavy neutrinos can be as low as ∼ 100 GeV. At the same time, leptogenesis via oscillations was primarily studied in the few-GeV regime. In Ref. [99] it was suggested that the maximal HNL mass in this mechanism is around the W -boson mass, where decays into W and Z bosons become kinematically allowed, and enhance the equilibration of the lepton asymmetries. The argument came about as follows. If the HNLs are kinematically allowed to decay into W (or W and Z), the rate of this process will exceed the Hubble rate well before the moment of sphaleron freeze-out. This means that the HNLs will be in thermal equilibrium and all asymmetries will be washed out. However, this is not the end of the story. For T < ∼ M , the heavy neutrinos begin to freeze-out, and their abundance, as well as the washout of the lepton asymmetries become Boltzmann suppressed.
Due to the different approximations that were applied to these mechanisms, the overlap between them remained an open question.
These mechanisms may appear quite different at first glance-but, as we show in this work, it turns out that the same equations can be used to describe both mechanisms. 1 In fact, the resonant enhancement has been already noted in [7].
In this work we systematically study generation of the BAU in the model with two HNLs with masses in the range 10 − 10 4 GeV. We find that the parameter space of the baryogenesis via neutrino oscillations is seamlessly connected with the parameter space of resonant leptogenesis. So there is just one mechanism, with different regimes. These regimes are characterized by whether the majority of the asymmetry is produced during the freeze-in or freeze-out of the HNLs (c.f. Fig. 1). We identify three main reasons why leptogenesis remains viable for M ∼ M W,Z .
• Although the equilibration rate of the heavy neutrinos generically exceeds the Hubble rate at the temperature of the sphaleron freeze-out, the lepton flavor washout rate can be suppressed, thereby preserving the BAU in one of the leptonic flavors. As was pointed out in [100], this allows for freeze-in leptogenesis with HNL masses above the electroweak scale. 2 • Due to the finite mass of the heavy neutrinos, their equilibrium distribution changes with temperature. This effect acts as a source for a deviation from equilibrium. We find that this effect remains important even for heavy neutrinos with masses at the GeV scale. This confirms the possibility of GeV-scale freeze-out leptogenesis as suggested in [36].
• We find that another factor preventing washout of the lepton asymmetries is approximate generalized lepton number conservation. 3 The processes that would erase the lepton asymmetry are suppressed by a factor ∼ ∆M/Γ, as identified in [102,103].
The short account of the results was presented in [54]. The paper is organized as follows. We start from an overview of the evolution of the leptogenesis calculations in section II. Then we define our notations and introduce the seesaw mechanism in section III. In section IV we describe resonant leptogenesis and baryogenesis via neutrino oscillations. We present the quantum kinetic equations for matrices of densities which at the core of our numerical study. We also show how the usual Boltzmann equations can be obtained as a limit of the quantum kinetic equations. Section V is dedicated to the determination of the heavy 2 A more detailed study of the effect of freeze-in on hierarchical leptogenesis followed in [101], which confirmed the importance of freeze-in in the conventional leptogenesis scenario. 3 In principle there are several ways of assigning lepton number to the heavy neutrinos. We discuss them in section C. neutrino production rates entering the kinetic equations. With all ingredients at hand, we perform the study of the parameter space in section VI. This section also contains detailed discussion of the obtained results. Next, in section VII we discuss several other works which considered the HNLs with masses around M W . We conclude in VIII. Technical details of the implementation, relation to the pseudo-Dirac basis, conserved lepton numbers, and fine tuning are discussed in the appendices.
FIG. 1. A sketch of the evolution of the HNL abundance in the early Universe. If we assume that the initial HNL abundance vanishes, there are two opportunities to generate the observed BAU.
The first is during a period of freeze-in, while the first HNLs are being produced and they approach equilibrium. The second opportunity is when the Universe cools down to temperatures below the HNL mass, and the HNLs decay out-of-equilibrium simultaneously with a freeze-out of the SM lepton number caused by the Boltzmann-suppressed washout rates.

II. CONVERGENCE TOWARDS A UNIFIED PICTURE
In this section we present a brief overview of the development of the calculations and methodology in low-scale leptogenesis.
The importance of a resonant enhancement for leptogenesis [10][11][12][13][14][15]18] was realized soon after leptogenesis was proposed as a baryogenesis mechanism. 4 Such a resonantly enhanced decay asymmetry offered an exciting opportunity-the mass scale of the HNLs could in principle be lowered to the electroweak scale [15,19,105] (it is worth noting that both the effects from a non-instantaneous freeze-out of sphalerons and spectator effects [106][107][108] were included in [105]).
However, the fact that the finite-order perturbation theory breaks down in the limit of degenerate HNL masses was already clear in the earliest papers [14][15][16], and different approaches of resolving these issues followed soon thereafter [18,109,110], with the general conclusion that the decay asymmetry is regulated by the width of the HNLs. This prompted further investigations into resonant leptogenesis using techniques of non-equilibrium Quantum Field Theory (see e.g. [111][112][113][114][115][116][117][118][119][120][121]). One of the most successful formalisms for in this framework is the CTP (closed-time-path) formalism, also known as the Schwinger-Keldysh formalism [111,112]. In this formulation of non-equilibrium QFT one typically deals with n-point functions of the different particles and uses them to calculate observables, such as particle numbers and distributions, and their time-evolution can be obtained by solving the Schwinger-Dyson equations on the CTP. One of the biggest differences compared to equilibrium QFT is that time translation invariance is explicitly broken (both by the boundary condition and by the expansion of the Universe), and each two-point function is a function of two time coordinates (and one momentum coordinate assuming translation invariance).
Consequentially, the Schwinger-Dyson equations are integro-differential equations which include both derivatives and integration over the time variables.
The main difference between the existing approaches is in the strategies used to solve these equations. Perhaps the most straightforward way is to solve the Schwinger-Dyson equations directly. However, this is only possible numerically or in specific limits. Nonetheless, this can be used to cross-check the results of other methods [60,61]. Another approach is to perform a Wigner transformation of the equations (it is often used to describe transport phenomena [122]), which leads to density-matrix like equations in [62]. This method was further developed to estimate the resonant enhancement in resonant leptogenesis [65,123], and was also used in [40] to derive the quantum Boltzmann equations used in leptogenesis via oscillations. 5 One of the shortcomings of this approach is that the off-diagonal matrix elements lie on an unphysical "average" energy shell E = (E 1 + E 2 )/2, and cannot be used for arbitrarily large HNL mass differences. 6 Recently, equations valid for arbitrary HNL 5 One should note that quantum Boltzmann equations are also used to describe decoherence effects in flavored leptogenesis [124][125][126][127]. 6 One should note that as soon as ∆M Γ, the standard Boltzmann equations may be used again, ensuring overlap between the two approximations since M Γ is typically satisfied in leptogenesis.
mass spectra were derived in [52]. The formalism developed in [128] instead works in the socalled two-momentum picture, also leads to similar quantum-Boltzmann equations [63,129], however, with an important difference-an mixing term that acts as an additional source of lepton asymmetries. The differences between these results prompted further investigations in [130,131].
Meanwhile, another line of research start to develop. The idea of baryogenesis via oscillations was put forward in work [20]. It was further developed in Ref. [21]. In particular, it was shown there, that the lepton back-reactions (missed in [20]) play an important role and the BAU generation is possible with two singlet fermions. This work introduced the kinetic equations for three lepton chemical potentials and two 2 × 2 HNL density matrices.
This equations were used in all consequent works on the baryogenesis via neutrino oscillations. Further clarifications of the mechanism were presented in [23,30]. In particular, the role of the process without helicity flip (we will refer to them as fermion number violating) was pointed out in Ref. [23]. It was noted that once such processes are included, the total asymmetry may be generated at the fourth order in Yukawa couplings. However, the relevant rates were not estimated correctly at that time. The neutrality of cosmic plasma has been accounted for in [31] by introducing the so-called susceptibility matrices relating the leptonic chemical potentials to the number densities. The temperature dependence of the susceptibility matrices and corrections coming from fermion masses were introduced in [35].
An important step has been performed in refs. [35,44], where the role of the fermion number violating processes was systematically accounted for and the equations were extended to the Higgs phase. Effects related to the non-instantaneous freeze-out of sphalerons were discussed in refs. [46,48]. The equations of [35] were further improved in [50,51,53].
It is interesting that the oscillations found in [63] were initially assumed to be a genuinely different source of asymmetry than the one in the scenario of leptognesis via oscillations, since they already appeared at fourth order in the Yukawa couplings (O(F 4 )), instead of O(F 6 ) as reported in the initial papers [20,21]. As discussed above, this apparent discrepancy was a result of the relativistic approximations, where LNV terms suppressed by a factor (M/T ) 2 were neglected. Early investigations [36] into non-relativistic corrections suggested that these terms are already important for GeV-scale HNLs, and may allow for freeze-out leptogenesis with GeV-scale HNLs.
In spite of the gradual convergence of the methods and results, a conclusive study showing that the parameter space of the two regimes of leptogenesis are connected was missing.
Strong hints of such a possibility were already provided in [100], where it was shown that freeze-in leptogenesis remains important for HNL masses above the TeV scale. Similarly, the study in [41] suggests that leptogenesis via HNL oscillations extends to HNL masses as large as 100 GeV. On the other hand in [36], it was suggested that already moderate decay asymmetries allow for freeze-out leptogenesis at the GeV scale. More recently this was confirmed in a similar study including flavor effects [132]. In this work we perform a unified study including all the effects relevant for mass scales up to a few TeV.

III. THE SEESAW FORMULA AND THE LIGHT NEUTRINO MASSES
For completeness of the paper in this section we briefly review the type-I seesaw mechanism and how it constrains the properties of the heavy neutrinos. The Lagrangian of the model reads which gives us three light mass eigenstates and the heavy eigenstates with a mass matrix The number of the heavy eigenstates is the same as the number of the right-handed fields ν R I . At least two HNLs are need to explain the two observed mass splittings in the active neutrino sector. In this work we focus on the minimal scenario with two heavy neutrinos.
We label the physical states of heavy neutrinos as N 2 and N 3 7 and denote their masses M 2 and M 3 . Throughout this work we will be interested in the case when N 2,3 have close masses, i.e. |M 2 + M 3 | |M 2 − M 3 |. Therefore it will be convenient to use the average mass M and the mass splitting ∆M . In order to match the notations of Ref. [49], we define them through M 2 = M − ∆M, So strictly speaking ∆M is a half of the mass splitting.

A. Parametrization of the Yukawa couplings
The masses of the light neutrinos m ν are constrained by the neutrino oscillations experiments (we use the global fit [133]). Out of the 9 parameters in the light neutrino mass matrix, 5 are already measured: two mass differences, and three mixing angles. The remaining unknown parameters are the mass of the lightest neutrino, two Majorana phases, and the CP -violating phase δ. 8 In the model with two HNLs the lightest neutrino is massless (up to tiny loop corrections [137]). Therefore it makes sense to speak about the neutrino mass hierarchy rather than ordering. In what follows we refer to normal (inverted) mass hierarchy as NH (IH).
The measured low-energy parameters mean that the choice of heavy neutrino masses M and the Yukawa couplings F is not completely free. To take this into account, we can parametrize the neutrino Yukawa couplings using the Casas-Ibarra parametrization [138]: where the matrix m diag ν is the diagonal neutrino mass matrix (M M is already diagonal in our basis), U ν is the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix, and R is a complex orthogonal matrix RR T = 1. For the PMNS matrix we use the standard parametrization [139]: where U ±δ = diag(1, e ∓iδ/2 , e ±iδ/2 ), and the non-vanishing entries of V (αβ) for α = e, µ, τ In the case of two heavy neutrinos there is only one relevant Majorana phase in the PMNS matrix. We parametrize it as η = 1 2 (α 21 −α 31 ) for normal, and η = 1 2 α 21 for inverted neutrino mass hierarchy with η ∈ [0, 2π]. The light neutrino mass matrix m diag with m 1 = 0 for NH, and m 2 = 0 for IH.
In the model with two right-handed neutrinos the matrices R depend on the neutrino mass hierarchy are given by with a complex angle ω = Re ω + i Im ω, and the discrete parameter ξ = ±1. The change of the sign of ξ can be compensated by ω → −ω along with N 3 → −N 3 [125], so we fix ξ = +1.
It is sufficient to constrain, Re ω ∈ [0, π], as larger angles only change the overall sign of the Yukawa couplings.
To summarize, there are six free parameters of the theory. They are listed in  As a consequence of the seesaw mechanism, the heavy neutrino mass eigenstates are mixed with the doublet neutrinos, and can interact with the rest of the standard model, the interaction states ν Lα are superpositions of the mass eigenstates ν i and N I : The mixing angle between a heavy neutrino N I and the active neutrinos ν L α is expected to be small, and is approximately given by This mixing angle appears in amplitudes for heavy neutrino production, which means that the heavy neutrino production rate is suppressed by For phenomenological applications it is convenient to introduce the quantities which quantify the total mixing angle of a particular heavy neutrino, the total mixing to a particular flavor, and the overall mixing between the heavy and light neutrinos.
This total mixing is useful to characterize the overall suppression of the interactions of the HNLs.

IV. BARYOGENESIS THROUGH LEPTOGENESIS
As stated in the introduction, any leptogenesis mechanism relies on HNLs to satisfy two of the Sakharov conditions which are not met in the SM-CP violation, and a sufficient deviation from thermal equilibrium. In thermal leptogenesis the CP violation arises via the loop corrections to the HNL decay rates: These decay rates can differ for decays into leptons and anti-leptons, which leads to an overall lepton asymmetry. This difference is often parametrised by the decay asymmetry: The evolution of the lepton asymmetries and the HNL number densities is an out-ofequilibrium process which can be described by the following Boltzmann equations (see e.g. [108,124,140]): We used the notation z = T ref /T (where T ref is some reference temperature -either the HNL mass M or the sphaleron temperature T sph ), Γ N I is the HNL decay rate, W α is the lepton number washout rate, µ α are the chemical potentials in the lepton flavors, and Y X = n X /s are the yields-the ratios of the number density to the entropy density. In the case of hierarchical HNL masses, the decay asymmetry α I is bounded from above [8], which translates to a lower bound on the mass of the lightest HNL: known as the Davidson-Ibarra bound.
Leptogenesis with HNLs below this mass scale can either be realized through resonant leptogenesis [19] or through leptogenesis via neutrino oscillations [20,21]. These two mechanisms were discovered independently for the different possible masses of the heavy neutrinos, leptogenesis via oscillations for GeV-scale HNLs and resonant leptogenesis for HNLs with masses from the TeV scale to the Davidson-Ibarra bound. In the following sections we briefly overview the main features of the two mechanisms and discuss the similarities and differences between them.

A. Resonant Leptogenesis
The Davidson-Ibarra bound on the decay asymmetry from [8] relaxes (as does the bound on the HNL mass) if the HNL masses are degenerate [10][11][12][13][14][15]18]. This is a consequence of the resonant enhancement in the CP -violating decays, which was discussed even before the idea of leptogenesis in [104].
This resonant contribution to the decay asymmetry comes from the interference of the tree level and wave-function amplitudes (the first and the last diagram in (14)).
which is enhanced when M 3 → M 2 . In the exactly degenerate limit the decay asymmetry diverges in the one-loop approximation. This apparent divergence is as an artifact of applying the S-matrix theory to unstable HNLs, which cannot be asymptotic states. To find the correct size of such a decay asymmetry, in the past two decades there were several studies of leptogenesis using out-of-equilibrium QFT methods [57-66, 128, 141-148].
The decay asymmetry then obtains a finite regulator of order of the decay width of the heavy neutrinos A ∼ M I Γ I , and The different approaches of deriving the Boltzmann equations and the decay asymmetries from first principles do not always agree and can lead to O(1) different results in the degenerate limit, i.e. when ∆M ∼ Γ I . 9 In practice the decay asymmetries combined with the Boltzmann equations (even if derived using out-of-equilibrium QFT methods) do not accurately describe the HNL dynamics. Instead one has to rely on quantum-kinetic equations that include the HNL coherence terms as dynamical degrees of freedom as we discuss in the following section.

B. Leptogenesis via neutrino oscillations
In leptogenesis via HNL oscillations [20,21] the asymmetry is not produced in the decays of the HNLs, but instead while the HNLs are produced and approach equilibrium in the early Universe. The HNLs generated via the SM particle interactions are produced in their interaction basis, which does not necessarily coincide with their mass basis. Due to this 9 One should however note that the two-time formulation [60,61,149] (which relies on very few approximations), and the so-called Wigner space approach [62,65] (where all interaction rates are evaluated on the same HNL average mass shell) lead to excellent agreement when the masses are not hierarchical [150]. misalignment, the HNLs begin to oscillate and through these CP violating oscillations they generate the lepton asymmetry.
Since the GeV-scale HNLs in this mechanism remain relativistic at temperatures relevant for baryogenesis T ≥ T sph ≈ 130 GeV, the caluclation of the BAU does not rely on calculating their decay rates as in resonant leptogenesis. Instead, a more apropriate picture is that of neutrino oscillations. To correctly take these processes into account, in [20,21], the equations used to describe the kinetic evolution of neutrinos due to Raffelt and Sigl [151] were modified to include oscillations between HNLs.
a. Evolution equations for baryogenesis via neutrino oscillations. Compared to the initial developments [20,21], there have been several systematic improvements to the kinetic equations for HNLs. One of the most significant improvements in recent years is the inclusion of corrections caused by the finite HNL mass [44,45]. On the other hand, it was also found that the same equations arise in the non-equilibrium formulation of QFT [40,47]. The key difference compared to the Boltzmann equations (17) is that besides the HNL number density Y N I , one also has to keep track of the HNL correlations. This information is encoded in the (modified from [40,44,45,47,62] to be valid in both relativistic and non-relativistic limits, c.f. [51][52][53]) including both the positive (negative) HNL ρ N (ρ N ) helicities, and the leptonic asymmetries n ∆α are given by: This distribution function is temperature-dependent. Its time derivative acts as a source of the 10 The equations in the current form-with lepton chemical potentials-are valid as long as leptons stay in equilibrium. In particular, for temperatures above T > 85 TeV , the right-handed electrons are not in equilibrium [152], and susceptibility matrices relating chemical potentials with number densities need to be modified accordingly.
deviation from equilibrium, therefore in what follows we will refer to dρ eq /dt as the source term. Note that we omit any Hubble expansion terms as we implicitly consider the comoving densities. The effective Hamiltonian describing the coherent oscillations of the HNLs is where The quadratic combinations of Yukawa couplings Y ± appearing in Eq. (22) and in the rates are defined below. The damping rates are The communication terms, describing the transitions from HNLs to active neutrinos, arẽ In the expressions above the subscripts + and − refer to the fermion number conserving and violating quantities correspondingly. The functions h ± and γ ± depend only on kinematics (i.e. on the common mass of HNLs). These functions have to be determined over the whole temperature region of the interest, which includes both symmetric and Higgs phases.
The dependence on the Yukawa coupling constants factorizes out from the rates (22)-(24), and we have where the matrix G encodes the generalized Majorana condition between the Majorana fields If we transform N I from the basis where the Majorana condition is given byÑ I =Ñ C I by a unitary transformation N I = U IJÑJ , this matrix is given by G = U † U * . The relation between the number densities and the chemical potentials to leptons in Eq. (21) has to take into account the neutrality of plasma. When the system is in equilibrium with respect to sphaleron processes, this relation reads where ω αβ (T ) is the so-called susceptibility matrix, see, e.g. [35,46].

C. Boltzmann equations as a limit of the quantum kinetic equations
As mentioned above, one of the key quantities in thermal leptogenesis are the decay asymmetries . These quantities are typically calculated via Feynman diagrams as shown in Eq. (14), but can also be obtained as a limit of the full density matrix equations. In this section we show how these quantities arise from the density-matrix approach, and discuss when such approximations can be applied. To find approximate expressions, it is instructive to explicitly write the equations governing the HNL correlations ρ IJ for I = J (note that we omitted the subscript N for brevity): If the diagonal elements δρ II , δρ JJ change on a time-scale much slower than the oscillation period, we can approximate the off-diagonal correlations from (21) by the steady state limit (obtained by setting dδρ IJ dt → 0): where the analogous solutions for δρ are be obtained by replacing H → −H, and ρ →ρ.
We may insert this solution into the source term for the lepton asymmetry from Eq. (21a): where we neglect higher order terms in the Yukawa couplings, including the helicity asymmetryρ II − ρ II . The two terms can be identified with the lepton flavor violating and lepton number violating source terms for the asymmetry.
Substituting this solution into Eq. (21a), and inserting the vacuum expressions for H gives us an approximate expression for the two decay asymmetries and Notice that α LFV I only violates lepton flavor, as the sum over α leads to a vanishing lepton asymmetry.
A few comments on these expressions are in order: • in contrast to the asymmetries used in [36,132], the decay rates entering (31) and (32) include explicitly the helicity-dependence of the rates, and lead to a result similar to [101] in the hierarchical limit, • in the non-relativistic limit γ ± → M/(16π) and |k 0 | → M we recover the relation (20), • in the same limit these expressions qualitatively agree with the regulators from [60], as well as those from [65] • although the decay asymmetry from Eqs. (31) and (32) remain finite when ∆M 2 the approximation of fast oscillations which we used here breaks down, and one has to rely on quantum kinetic equations (21).
a. Applicability of the decay asymmetries. We conclude this section with a short comment on the different ways of regulating the decay asymmetry LFV . To explore the breakdown of Boltzmann equations in the GeV-mass regime, we recover the result from [20,21], following the procedure described in [29]. Taken at face value, Eqs. (31) and (32) appear well-behaved in the limit of ∆M N → 0, but once we take the full picture of heavy neutrino oscillations into account we find that this is not the whole story. In the relativistic limit T M , the rate γ + ∼ T , and γ − ≈ 0, equation (32) simplifies to: In contrast to the temperature independent decay asymmetries often used for resonant leptogenesis, this decay asymmetry has a dramatic temperature dependence, and becomes enhanced when M/T → 0: where z ≡ M/T . The question is what prevents the initial lepton asymmetries from diverging, we might expect that the asymmetry is already made finite by the width in Eq. (32). On the other hand, if we solve the kinetic equations, we find that another effect suppresses the decay asymmetries much sooner: the finite time required for HNL oscillations. As shown in figure  2. An example time evolution of the off-diagonal correlations which enters the decay asymmetry (blue, full). In the limit of fast oscillations, this highly oscillatory behavior is well described by the static limit (green, full). The fast oscillations lead to a numerically stiff system. In (red, dashed) we show how the full solution approaches the static limit once we average it over some a finite time interval. This example shows that using the static solution before a single oscillation period can lead to a significant overestimate of the BAU, even when it is regulated by a finite width.
section, a systematic treatment of the low-scale leptogenesis in all corners of the parameter space is only possible within the framework of quantum kinetic equations (21). The crucial ingredient of these equations are the rate coefficients entering equations (22)- (24). The next section is dedicated to the determination of these coefficients in all relevant regimes.

V. THE HEAVY NEUTRINO PRODUCTION RATES
In recent years significant progress has been made towards determining the production rate of the heavy neutrinos in the early universe [25,27,101,[153][154][155]. In particular, the production rate of GeV-scale heavy neutrinos at temperatures T M has been studied in great detail [35,45,48]. For GeV-scale neutrinos important effects arise at temperatures below the electroweak crossover, where the mixing between the heavy and light neutrinos can significantly affect the production rate [44,45]. The active neutrino interaction ratewhich affects the heavy neutrino production at these regime-has even been calculated at NLO [156].
At present, there are no readily available estimates of the full heavy neutrino production rate for T ∼ M . In [101], the rate was calculated only including the naive leading order 1 ↔ 2, and 2 ↔ 2 processes, thereby neglecting the 1 + n → 2 + n processes enhanced by multiple soft scatterings, also known as the LPM (Landau-Pomeranchuk-Migdal) effect.
On the other hand, earlier calculations [25,27,153], only provide the helicity-averaged rate, which is not sufficient to track the full evolution of heavy neutrinos as the universe cools down from T M to T ∼ M . Below we describe our approach to the calculation of the heavy neutrino production rate.
The production rate of the heavy neutrinos can be expressed through the spectral (antihermitian 11 ) part of their self-energy Σ A N . To simplify the calculation, we factor out the Yukawa couplings from the self energies and introduce where g w = 2 counts the SU (2) degrees of freedom, and the symmetric matrix G encodes the generalized Majorana condition N = GN c , and must be included to ensure flavor covariance of the equations.
The rate coefficients γ ± are related to the self-energy as where we define Σ ± and k ± as a ± = a 0 ± a i k i /|k|. The results of [153], as well as [35] provide a calculation of the spin-averaged heavy neutrino production rate, which is related to the 11 In literature this is also known as the imaginary part of the self-energy, which is strictly speaking not accurate for heavy neutrinos whose self-energy is a complex matrix in Dirac and flavor spaces.
rates γ ± as where Im Π-in the notations of [35]-is the antihermitian part of self-energy.
For M T , when neutrinos are relativistic, and k − ∼ M 2 2k the γ − rate is suppressed by M 2 /T 2 , and Im Π ≈ γ + corresponds to the total heavy neutrino production rate. In [45], these rates were expressed in terms of Q ± with On the other hand, for M T , the two rates γ ± become equal, as the heavy neutrinos become non-relativistic and both k + ≈ k − ≈ M , andΣ N ± ≈ M 32π . In the intermediate regime, for temperatures T ∼ M we estimate the full rate by extrapolating the relativistic rate from [45], and combining it with the rate of heavy neutrino decays above the "thresholds" 12 corresponding to decays into the Higgs, W ± and Z bosons for M N > M H,Z,W .
We perform the extrapolation as follows, we note that the heavy neutrino production rate Γ N = γ + + γ − appears to have a rather weak dependence on the heavy neutrino mass in the results from [27,153] However, we note that the dependence does not directly appear in the self-energyΣ, but mainly through the prefactor k − /|k 0 | ∼ M 2 /k 2 .
To summarize, we estimate the self-energy of the heavy neutrino as the sum of a Mindependent term, and a M -dependent term for masses above the heavy neutrino decay If we compare this approximation to the naive 1 ↔ 2 decay rate, we find that there is a clear M -dependence in the 1 ↔ 2 rate even for M < m φ . This dependence is caused by the finite size of the effective lepton mass m , and the fact that the decay of the heavy neutrino would be forbidden for m φ − m < M < m φ + m . Note however that this kinematically forbidden region is an artifact of our approximation when we treat the effective mass of the lepton as a physical mass term. Such kinematically forbidden regions in reality disappear once the 2 ↔ 2 scatterings, as well as the LPM effect are included.
The contributions from the 1 ↔ 2 processes to the antihermitian part of the heavy neutrino self-energy (c.f. Fig 3) are given bŷ where ∆ A φ (k) andŜ A (k) are the spectral functions. To obtain the naive result for the 1 ↔ 2 heavy neutrino production rate, we replace them by the tree-level approximation where m 2 = T 2 (g 2 1 + 3g 2 2 )/16 and m 2 φ = T 2 (g 2 1 + 3g 2 2 + 4h 2 t + 8λ)/16. Note that we ommit the chiral projector inŜ , as we factored it out in the definition ofΣ N .
When the on-shell condition is imposed through the delta functions from (41), the integral in Eq. (40) reduces to a one-dimensional integral with a well known result where the functions I n correspond to the integrals with the integration boundaries A. Production rates after electroweak symmetry breaking Following Ref. [35], we can identify two classes of processes that lead to the production of the heavy neutrinos. Direct processes, where the heavy neutrino interacts with Higgs and lepton doublets, as well as indirect processes in which the heavy neutrino interacts through the mixing with the light neutrinos.
One of the biggest differences compared to the symmetric phase arises in the decay of the heavy neutrino, as the Higgs doublet is replaced by a real scalar field. At the same time, decay channels into the Z and W ± bosons open up. There is a gauge dependence in how we separate the direct and indirect processes, but the total rate in their sum remains gauge invariant. This gauge invariance can be lost if we do not keep all the direct and indirect diagrams at the same order in perturbation theory. As the resonant mixing between the heavy and light neutrinos is best described by resumming the lepton propagators, we implicitly include diagrams that are not matched by the direct rate, and therefore introduce gauge dependence to our estimate of the heavy neutrino production rate. We resolve this issue by choosing a specific gauge, in this case the 't Hooft-Feynman gauge, due to the fact that it does not introduce any new scales.

Direct Heavy Neutrino Production
In this section we perform a simple estimate of the direct heavy neutrino production rate. As discussed at the beginning of this section, there are various contributions to the The situation slightly changes in the broken phase of the standard model. In the 't Hooft-Feynman gauge, the goldstone bosons also contribute to the heavy neutrino production rate.
We can approximate this rate by the same 1 ↔ 2 integrals from the previous section however, with m φ replaced by the mass of the appropriate goldstone boson (c.f. Fig. 6) i.e. the gauge  boson to which the goldstone mode correspondŝ The individual rates were divided by factors of 2g w , to take into account that only one isospin runs in the loop, as well as the 1/ √ 2 factor that appears in the coupling to the field h compared to φ.
Despite the fact that our calculation of the 1 ↔ 2 rate misses the 1 + n ↔ 2 + n processes which are included in the LPM resummation, our results-after averaging over helicitiesshow a nice agreement with the full helicity averaged computation, see figure 5.
The remaining direct processes, such as 2 ↔ 2 scatterings can also be included in such An accurate estimate of the non-perturbative contribution to masses of the Higgs and vector bosons typically requires a lattice calculation [157,158]. To ensure consistency with the rate calculation from [50], we follow their approach and rely on a 1-loop calculation.

Indirect Heavy Neutrino Production
At tree level the sum of the direct and indirect rates can be shown to be gauge invariant.
However, if we resum the propagator of the doublet leptons, we implicitly include a whole set of gauge-dependent diagrams, which do not have their match in the tree-level diagrams of the direct processes. The gauge dependence of resummed propagators is well documented in the literature (see e.g. [159]), and can be alleviated by using e.g. higher order effective actions. Note that in R ξ gauges the mass of the goldstones is given as ξ × m gb , where m gb is the mass of the corresponding gauge boson. Therefore, if we choose ξ = 1, i.e. the 't Hooft-Feynman gauge, we avoid introducing new scales into the problem.
The resummed lepton propagators are given by (see e.g. [160]): where the superscripts H and A stand for the hermitian (dispersive) and antihermitian (dissipative) parts of the lepton self-energy. These expressions significantly simplify if we use light-cone coordinates, where the denominator factorizes, which effectively gives us two propagators, one for the particles, and the other for the holes in the plasma. For the retarded propagator this gives us: where the indices ± indicate the light-cone coordinates a ± = (a 0 ± a i k i /|k|). The indirect FNV and FNC rates are then given by [44,45]: where Σ are the self-energies of the neutrinos. Interestingly, for large HNL masses, these self-energies have to be evaluated on the HNL mass shell, which leads to further uncertainties, as the active neutrino interaction rates have only been calculated for relativistic neutrinos.
To estimate this rate we use the same approach as for the rates in the symmetric phase, namely, we add the 1 → 2 contribution to the relativistic estimate when the HNL mass exceeds the decay thresholds. Fortunately, any such uncertainties are not important for larger HNL masses, since the off-shell active neutrino width does not grow faster than the HNL mass in the denominator itself, and the indirect contribution becomes negligible for M T , as shown in Fig. 8.
By combining the estimates of the rates we have a consistent set of rates in both the broken and symmetric phases for any HNL mass. In spite of the uncertainties related to our extrapolation, our rates remain consistent with the spin-averaged rates from the literature as shown in Fig. 5. For completeness we show both the self-energies (upper panels) and the interaction rates (lower panels).

A. Baryogenesis through freeze-in and freeze-out
As we have seen in the previous chapters, the same equations can be used to describe both leptogenesis through neutrino oscillations, and resonant leptogenesis. Because there is no clear cut between the two mechanisms, we may instead look at whether the majority of the asymmetry is produced during freeze-in or freeze-out of the heavy neutrinos.
For a clear separation between the two regimes, we perform three parameter scans: 1. Parameter scan with vanishing initial abundance of the heavy neutrinosboth freezein and freeze-out contribute to the BAU generation.
2. Parameter scan with thermal initial conditionsonly freeze-out contributes. Indeed, by freeze-in we understand the period during which HNLs reach the equilibrium. By artificially setting the thermal initial conditions we eliminate this period.
3. Parameter scan with vanishing initial abundance of the heavy neutrinos and without the deviation from equilibrium caused by the expansion of the Universeonly freezein contributes. Technically this is achieved by putting to zero the source term dρ eq N /dt from Eqs. (21).
Because of the approximate linearity of the evolution equations, the three parameter scans are not independent, and the baryon asymmetry from both freeze-in and freeze-out can be obtained by summing the other two parameter scans. However, to avoid numerical errors due to cancellations of large numbers, we perform all three parameter scans independently.

B. The range of allowed masses and mixing angles
The parameter space of leptogenesis is quite large, as we have discussed in section III A, in total we have six unknown parameters. These parameters are: the average mass M ; the mass splitting ∆M ; one Dirac and one Majorana phases of the PMNS matrix; the real and imaginary parts of the complex angle ω. Instead of constraining these parameters directly, it is instructive to see how they are related to the experimentally observable quantities. Among these the most relevant are the masses of the heavy neutrinos and their mixing angles. This is particularly important for direct searches as the size of the mixing angle determines the number of heavy neutrinos that can be produced.
The condition of reproducing both the baryon asymmetry of the Universe and the light neutrino masses imposes constraints on the allowed mixing angles of the heavy neutrinos.
The lower bound on the total mixing angle U 2 is mainly coming from the requirement of reproducing the light neutrino masses, while the upper bound comes from leptogenesis.
The mixing angle of the heavy neutrinos is tied to the size of their Yukawa couplings, see Eq. (10). For large values of the Yukawa couplings, the mixing angle itself will be large.
The main issue for leptogenesis is that a large value of Yukawa couplings at the same time leads to a large washout strength.
To perform the study of the parameter space we choose several specific benchmark points which extremize the ratios of the Yukawa couplings It is interesting that these points correspond to special values of the PMNS phases, with 13 δ = nπ/2 , η = mπ/2 , Re ω = lπ/4 with m, n, l ∈ Z .
The same phases have been show to maximize the total mixing U 2 in Ref. [42,49]. These choices of parameters are particularly favorable for leptogenesis as they allow for a hierarchy in the washout strengths. This means that even in the presence of a large overall washout parameter, the asymmetry may survive hidden in a particular lepton flavor.
Once the phases δ, η, and Re ω are fixed, we scan over the remaining parameters. The details of the numerical procedure are specified in A. We limit the HNL mass to the range

C. Constraints on the heavy neutrino mass splitting
The mass splitting between the heavy neutrinos is one of the most important parameters for both leptogenesis scenarios. The main reason why leptogenesis is so sensitive to the mass splitting ∆M between the heavy neutrinos is that this parameter sets the scale for the oscillations that violate CP and lead to a lepton asymmetry. The temperature corresponding to the onset of oscillations depends on the Hubble rate and is given as [21] T where M 0 = 90/ (8π 3 g * )M Pl and g * is the effective number of relativistic degrees of freedom. For heavier neutrinos, it is possible that the oscillations begin when they are already non-relativistic, which gives us a different temperature since the typical HNL energy is M This oscillation temperature scale can be lowered once we take the Yukawa-induced thermal masses into account, as the thermal masses are aligned with the HNL interaction basis, and do not lead to oscillations themselves. To understand the parameter space it is instructive to consider "slices" where we vary the mass splitting and the magnitude of the Yukawa couplings (i.e. the parameter Im ω), and keep the remaining parameters fixed.
The results of these parameter scans are shown in figure 10. For each mass we show three lines, one corresponding to leptogenesis via freeze-in, one to leptogenesis via freeze-out, and one containing both contributions. This helps us better understand the transition between the two regimes. produced during the freeze-in and oscillations of the heavy neutrinos [99]. However, the lepton number washout does not remain large indefinitely, as it also depends on the equilibrium density of the heavy neutrinos. This gives us the usual Boltzmann suppression factor [98] for the lepton number washout The lepton number washout rate therefore reaches its maximum when T ∼ M , and the  Re ω = π/4, δ = 0, η = π/2. The different contours correspond to different initial conditions for the heavy neutrinos. The area inside the regions corresponds to a BAU greater than the observed asymmetry. The (dark blue, full) curve includes both the contributions from freeze-in and freezeout. The (light blue, dashed) curve corresponds to freeze-out only and the (red, dotted) curve corresponds to baryogenesis via freeze-in. It is interesting that the largest mass splitting is realized exactly during freeze-in, as the HNL oscillations happen at high temperatures, and therefore before the HNLs begin to decay.
which we average over the two (almost degenerate) heavy neutrinosK = (K 1 + K 2 )/2. The value ofK indicates whether the washout is strong 1 or weak 1. It is interesting that the washout rate exceeds the Hubble rate by at least a factor O (30), which implies that the washout is already strong for small values of Im ω, and becomes even larger when At the same time, the expansion of the universe drives the heavy neutrinos out of equilibrium. When the universe cools down to T ∼ M , the heavy neutrinos are no longer relativistic, and they become over-abundant compared to their Boltzmann-suppressed equilibrium density. This late-time deviation from equilibrium is exactly what leads to freeze-out leptogenesis. It can typically be neglected in studies of leptogenesis with GeV-scale HNLs, since it is suppressed by a factor (M/T ) 2 for M < ∼ T , and only reaches its maximum when The results from the parameter scan in Section VI B imply that freeze-in and freeze-out leptogenesis parameter spaces have significant overlap: freeze-in leptognesis remains viable for masses much higher than what one could have expected, and freeze-out leptogenesis is already possible for masses as light as M ∼ 5 GeV. We find that the main reason why freeze-in leptogenesis remains possible is the difference of the flavored washout compared to the average washout strengthK, which we further discuss in section VI D 1. We show how freeze-out leptogenesis can be realized in a scenario with GeV-scale heavy neutrinos section VI D 2. Finally, we discuss why large mixing angles do not necessarily imply a large lepton number washout in section C.

Freeze-in leptogenesis for TeV-scale HNLs
For GeV-scale heavy neutrinos the washout of the baryon number stops when the sphaleron processes freeze-out, i.e. at T sph ≈ 130 GeV. The washout of the asymmetry in a particular lepton flavor does not only depend on the washout parameterK, but also on how strongly a particular flavor couples to the heavy neutrinos. In a realistic computation, we also need to take into account the so-called spectator effects (see, e.g [35]), which redistribute the asymmetries among the remaining SM species. These effects are typically encoded in the susceptibility matrices ω αβ , as included in equation (21a).
The SM lepton number washout is determined by the smallest flavored washout W α from Eq. (17). To determine the size of the lepton number washout, we have to take the spectator effects into account via the susceptibility matrix ω αβ . The susceptibility matrix is close to diagonal, and we can estimate the SM lepton number washout strength as The key quantity governing the size of the lepton number washout compared to the heavy neutrino washout rate is therefore the minimal branching ratio This means that the effective lepton number washout strength is given by For normal ordering this implies that the washout can be weak to moderate when Im ω ≈ 0, whereas for inverted ordering it can also be strong -depending on the values of the CPphases η and δ. In figure 11 we show an example of the parameter space in the IH case with a choice of PMNS phases that leads to strong washout, where the freeze-in parameter space completely closes around 500 GeV, and the BAU for larger masses becomes independent of the initial conditions.

Freeze-out leptogenesis for GeV-scale HNLs
Another reason why the parameter space of resonant leptogenesis remains connected with leptogenesis through neutrino oscillations is that HNL freeze-out provides a sufficient deviation from equilibrium even if their masses are as low as a few GeV. To obtain a qualitative estimate of the BAU, we can use the decay asymmetries (31) from section IV C. The deviation from equilibrium caused by the expansion of the Universe is suppressed by a factor M 2 /T 2 : where s(T ) is the entropy density of the Universe. Re ω = π/4 is the same as in figure 9, as it does not affect the washout of the lepton number.
If we neglect the lepton number washout, in equation (17) only the LNV part of the decay asymmetry survives, which gives us an estimate of the BAU where we used the mass splitting that saturates the resonant enhancement ∆M ∼ Γ I , and only included the leading term in γ − .
This approximation gives us an estimate similar to the results from [36,132], but the origin of this M 4 suppression is in reality quite different. In [36], the additional M 2 /T 2 factor arises through the interplay of the thermal HNL masses, and the thermal decay rates, while here it is a direct result of the helicity-dependent rates γ ± , where the FNV rate is suppressed by a factor γ − ∼ M 2 /T 2 .
We should however note that this approximation has very limited applicability, and an accurate lower bound can only be obtained by using the full density matrix equations for the following reasons: • the mass splitting that is needed to saturate the decay asymmetry ∆M ∼ Γ I corresponds to an oscillation temperature very close to the sphaleron freeze-out temperature, which means that the approximation of fast oscillations does not hold , where Y denotes either Y + ot Y − from Eq. (26) and we used Y ∼ F 2 ∝ M for a fixed value of m ν , as can be seen from the parametrization (5). We note that the equilibration rate remains the same, while the bare Hamiltonian term scales inversely with the mass of the heavy neutrinos. Fortunately, the mass splitting ∆M in resonant leptogenesis is typically much smaller than the mass M . This effectively allows us to scale ∆M independently of M . In order to cancel the ζ scaling of H 0 we replace ∆M → ∆M ζ 2 . This only induces a small correction to the Yukawa couplings F → F + O(ζ 2 ∆M 2 /M 2 ). Finally, this means that we do not need to repeat the parameter scan for every mass above M > 2 TeV, as the results are related by an appropriate re-scaling of the mass splitting ∆M . This conclusion has several caveats: • For mass splittings of ∆M ∼ M , the scaling relation does not hold and the correction to F becomes significant.
• The up quark only reaches equilibrium at temperatures below 10 6 GeV, therefore if a significant portion of the lepton asymmetry is produced above these temperatures, the simple scaling solution will not be sufficient, as the susceptibility matrices change, and the distribution of the lepton asymmetries among the flavors may be modified.
However, we note that this is a small change in the number of degrees of freedom entering the equations, and we expect it to lead to no more than an O(1) change in the baryon asymmetry.
• For temperatures above 10 9 GeV, the "flavored" approximation breaks down, and we have to keep track of the lepton doublet correlation matrices.
• The typical energy scale of the interactions changes with temperature, which means that the running of the couplings has to be taken into account, which leads to a modification of the rates entering the kinetic equations.
• Finally, when solving the equations for the heavy neutrino oscillations we implicitly neglect the vertex diagrams which are important for thermal leptogenesis. However, the contribution of this diagram is subleading compared to the wave-function contribution below the Davidson-Ibarra bound.
We may conclude that the biggest differences occur when the majority of the BAU is produced at temperatures above 10 9 GeV. • T. Hambye and D. Teresi [36,37] The authors study resonant leptogenesis in the GeV-regime. They apply the language of decay asymmetries to GeV-scale heavy neutrinos and indicate that leptogenesis in both freeze-in and freeze-out are possible for GeV-scale HNLs. One should note that although the decay asymmetry obtained in this work qualitatively differs from the expressions in (31), the authors find the same parametric suppression of the LNV decay asymmetry of order O(M 2 /T 2 ). In [37], the authors used the density matrix approach, and found semi-quantitative agreement with the approach using Boltzmann equations. The authors also find that the LNV asymmetry scales with the product of the helicity-dependent rates γ + γ − , similar to the LNV source term in Eq. (31).
This study of parameter space is limited to M < 10 GeV. As, we summarize in section IV C, the Boltzmann equations can be understood as an approximate limit of the density matrix approach, however in a limited range of validity. A conclusive study of the leptogenesis parameter space in this regime requires a unified treatment that can reproduce both baryogenesis via oscillations and resonant leptogenesis in the appropriate limits (as we present in this work).
• A. Granelli, K. Moffat and S. Petcov [132] The authors consider resonant and flavored leptogenesis in the GeV-regime. They find that the flavored decay asymmetry can lead to a BAU several orders of magnitude larger than what the LNV contribution alone would imply. Due to the reliance on the Boltzmann equations, this approach also has limited applicability (in the same sense as [36]). One should also note that the flavored decay asymmetry studied here can be understood as a limit of fast oscillations from (32). If interpreted this way, these results can be understood as an alternative approach to leptogenesis via neutrino oscillations using the language of Boltzmann equations.
One should note that [36,37,132] often discuss two separate decay asymmetries-one coming from mixing, and the other one coming from the oscillations [63,66,129,130] of the HNLs.
We do not make such a separation, as no separate sources of the asymmetry appear in derivation of the equations using either the canonical Raffelt-Sigl formalism of neutrino oscillations (see e.g. [20,21,45,52]), nor in the CTP approach in Wigner space [40,47,62].
Because of this, we do not introduce any such additional terms by hand as this can lead to double-counting, and therefore violate the symmetries that are present in the equations in the relativistic limit (for a detailed account of the conserved charges see section C).

VIII. DISCUSSION AND CONCLUSIONS
In this work we investigated the similarities and differences between resonant leptogenesis and baryogenesis via neutrino oscillations in the minimal extension of the standard model by two heavy neutrinos. We found that the two mechanisms are closely related, and that the equations used to describe the two mechanisms are virtually the same. Since the defining feature of resonant leptogenesis, namely the resonant production of the baryon asymmetry is also present in baryogenesis via neutrino oscillations, we focus on the major difference between the two mechanisms, namely the question whether the majority of the BAU is produced during the freeze-in, or freeze-out of the heavy neutrinos.
We found significant overlap between the two regimes, namely, freeze-in leptogenesis turns out to play a major role in generating the BAU even for TeV and heavier Majorana neutrinos. This regime mainly coincides with relatively large ∆M/M 2 ∼ 10 −7 GeV −1 mass splitting, compared to the one optimal for a resonant enhancement ∆M/M 2 ∼ 10 −13 GeV −1 .
Furthermore, this also implies a strong dependence on the initial condition which is absent in freeze-out leptogenesis.
On the other hand, we also find viable freeze-out leptogenesis with masses as low as The equations governing the evolution of the heavy neutrino densities, as well as the lepton asymmetries depend on multiple scales, and can in principle be stiff and therefore numerically expensive. Since the parameter space spans many orders of magnitude, there is no single analytical approximation that can be applied everywhere in the parameter space.
Nonetheless, there are several simplifications of the equations that can lead to improvents in accuracy and speed of the numerical calculations. In the following we summarize the main approximations and simplifications that we use.
Momentum averaging. Equations (21) are integro-differential equations. In principle one has to solve a differential equation for each of the HNL momentum modes, and integrate over them at each time step to accurately describe the evolution of the lepton asymmetries.
Such a calculation was performed in [26,48,50], however, the numerical intensity of such calculations makes it unfeasible for large scale parameter scans. In this work we use a momentum-averaging procedure, whereby we only consider how the time evolution of the integrated HNL matrix of densities: We then assume that the density matrices remain proportional to the equilibrium distribution, but differ by a total normalization factor: which allows us to replace the momentum-dependent matrices ρ N in Eq. (21) by the HNL number n N . This also gives us a procedure to calculate the momentum-average of the HNL rate coefficients: for the coefficient in Eqs. (22,23,24): γ ± , h ± and 1/E N .
This approximation gives rise to some artifacts, like the oscillatory behavior of the lepton asymmetries which are not present in the full system, where the oscillations of the different momentum modes do not happen simultaneously. However, in [26,48,50] it was shown that the momentum-dependent treatment leads to O(1) corrections, which justifies the approach taken in this work.
Transformation of the equations into vector form. The variables in the system of Eqs. (21) are a combination of different objects-numerically, the lepton asymmetries form a real vector, while the HNL densities δρ N , and δρ N are hermitian matrices with complex elements. In general, matrices with complex elements will have n 2 s × 2 degrees of freedom, while a hermitian matrix only has n 2 s degrees of freedom (where n s is the number of HNLs). Any numerical ODE solver will treat these matrices as general complex matrices and introduce more degrees of freedom than are strictly necessary to solve the system. To avoid this, we may choose parametrize all hermitian matrices as real four-component vectors: Once we transform all the equations accordingly, the full system (21) can be written in matrix form:q where all degrees of freedom can be encoded in an 11-dimensional vector: Approximate integration of fast modes. In general, the system of equations (21) contains many different time scales. The presence of significantly different time scales can makes the system stiff, which limits the speed of the numerical computations.
Several strategies of reducing the stiffness of the system have been laid out in [40], where the fast equilibration modes can be integrated-out in the case of large mixing angles and small mass splittings, and the equations can even be solved semi-analytically.
Let us briefly overview the main idea of this calculation. Let us assume that we can separate the fast and slow degrees of freedom as q = (q F , q S ) T , witḣ where A AB are block matrices, and the eigenvalues of A F F are much bigger than the eigenvalues of A SS . Due to the large eigenvalues of A F F , the fast modes reach the quasi-static equilibrium: We can re-insert this solution into the equation for the slow modes to track the evolution of the rest of the system, which gives us a differential equation for the slow modes: In practice, it is however impractical to keep track of the fast and slow modes, as these change with temperature (see, e.g. [56]). What we can do instead is to replace the matrix A with one that approaches the same quasi-static limit for the fast modes, but leaves the slow modes intact.
One such ansatz is: where we choose the cutoff scale Λ between the fast and the slow scales This choice ensures that the fast modes approach the same quasi-static limit as in Eq. (A9), but at the time-scale Λ, i.e. for A F F Λ we have: whereas for the slow modes we recover Eq. (A10).
What is perhaps most important is that we do not need to block-diagonalize the system into fast and slow modes at each point in time as the matrix A evolves. Instead, the fast and slow modes are automatically separated based on how they compare to the scale Λ.
Another useful feature of this approach is that it equally applies when the HNL oscillations become fast, since it does not depend on A F F being real. This ensures that the quasi-static limit of the off-diagonal source term is the same as in (29), but also prevents the system from becoming too stiff.
Appendix B: Pseudo-Dirac limit In this work we are interested in resonant leptogenesis and baryogenesis via oscillations which both require that the two HNLs have nearly equal masses. In this limit there is an approxiamte U (1) symmetry in the theory (two exactly degenerate Majorana particles form a Dirac one with the associated U (1)). This symmetry allows us to assign a lepton number to different flavors of the heavy neutrinos, and to effectively combine the two Majorana neutrinos into a single Dirac neutrino.
and the Lagrangian takes the form This limit is also realized in the Yukawa couplings. For large values of | Im ω|, they are much bigger than the naive estimate In such a scenario it is convenient to introduce the matrix of Yukawa couplings h αI related to the matrix F αI defined in (1) as follows where the couplings h αI are manifestly hierarchical with α |h α2 | 2 α |h α3 | 2 = e 4 Im ω + O(M 3 − M 2 ) .

(B6)
Appendix C: Conserved lepton numbers One of the main effects that lead to large values of the mixing angle even for U 2 m ν /M is the approximate conservation of lepton number.
There are several ways of assigning lepton number to heavy neutrinos so that it is approximately conserved. In the following we will discuss three that are important for leptogenesis.
a. Helicity as a lepton number The first, and most commonly discussed in low-scale leptogenesis mechanisms is the association of a lepton number with the helicity of the heavy neutrinos. This regime effectively corresponds to heavy neutrinos being relativistic so that their masses can be neglected. It was discussed already in the work [20], and is the main reason why three heavy neutrinos were required for successful baryogenesis via oscillations.
When this was further developed in [21], it was noted that this conservation is not necessarily a barrier for a production of a total lepton asymmetry as a part of the total (conserved) lepton asymmetry is hidden in the heavy neutrino sector. Furthermore, it is interesting that this lepton number remains conserved regardless of the number of the heavy neutrinos, as long as they are all relativistic.
This lepton number conservation is a direct consequence of one of the helicity-dependent rates dominating, in the case when γ + γ + , we can verify that αΓ α ≈ −Γ. Summing over the lepton asymmetries and the HNL helicity asymmetry, we find the conserved combination: which can be verified by inserting the expressions from (21) into the derivatives of the components. On the other hand, in the broken phase, the FNV rate γ − can be orders of magnitude larger than γ + , giving us αΓ α ≈ Γ. This gives us a different assignment of lepton number [44]: b. Approximate Dirac lepton number There is another type of lepton number that is typically considered in low-scale seesaw mechanisms. Instead of assigning a lepton number to all right-handed neutrinos, we may assign different lepton numbers to the superpositions of the heavy neutrino flavor eigenstates. We effectively combine the two Majorana spinors into a single Dirac spinor. There are however terms that violate this symmetry, namely the mass splitting between the two heavy neutrinos, and particular structure of the Yukawa couplings.
The conserved combination is most evident in the non-relativistic limit (when γ + = γ − , as well as h + = h − ). In the pseudo-Dirac limit we can neglect the mass difference between the HNLs, as well as the terms suppressed by e −4| Im ω| in the Yukawa couplings. This makes the equilibration matrix proportional to the identity (Γ ∼ 1 2×2 ), and all of theΓ α ∼Γ are proportional to each other. We can introduce a matrix κ ≡Γ/Γ, which selects the combination of HNL flavors that carries the lepton number.
We can then verify that the combination: as long as the superpositions of the HNL flavors remain coherent.
This effect was taken into account in [102,103], and leads to a suppression of the flavored washout parameter (which enter the washout term in the Boltzmann Equation (17)): It is important to note that this effect is already included in the density-matrix equations used here, which do not need to be modified further. and (D2). The allowed range of mass splittings and mixing angles for a benchmark point with normal hierarchy and parameters Re ω = π/4,δ = π, η = 3π/2 for NH, and Re ω = π/4,δ = 0, η = π/2 for IH. The (yellow, full) line shows the upper limit on the mass splitting (D1), and the (green, dashed) line shows the lower bound on the mass splitting (D2). For comparison we also show the magnitude of the mass splitting that is induced by the Higgs vev corrections to the HNL mass matrix. The Lagrangian mass splittings smaller than this value are completely obscured by this contribution, and because of this cannot be measured in any experiment.