Parameter space of baryogenesis in the νMSM

The Standard Model accompanied with two right-handed neutrinos with masses below the weak scale can explain the observed baryon asymmetry of the Universe. Moreover, this model is at least partially testable in the forthcoming experiments such as NA62, SHiP, and MATHUSLA. The remarkable progress in understanding of various rates entering the kinetic equations describing the asymmetry generation along with considerable improvements of the numerical procedures allow us to perform a comprehensive analysis of the parameter space of the model. We find that the region of parameters leading to successful baryogenesis is notably larger than it was previously obtained for light HNLs. Our results are presented in a way that they can be readily used for studies of sensitivity of various experiments searching for the right-handed neutrinos responsible for the baryon asymmetry of the Universe. We also present a detailed comparison with the studies by other groups.


JHEP07(2019)077
obtained. The results are presented in a way that they can be used for a detailed study of sensitivity of different experiments.
The paper is organized as follows. First, we introduce the νMSM and the parameters the model in section 2. Then we describe the experimentally relevant quantities in section 3 and present the cosmologically favourable values of these quantities in section 4. These values are determined by imposing the requirement of the successful baryogensis. In section 5 we provide all necessary information on the open-access datasets. All theoretical and technical details are presented in the subsequent sections. In section 6 we overview the kinetic equations derived in ref. [36]. We discuss our approach for the numerical solution of these equations and describe the impact of the improvements in section 7. The study of the parameter space is performed in section 8. Section 9 contains a detailed comparison with the works [21, 23, 32-34, 37, 40]. We summarise in section 10. Appendix A describes the mixings of active neutrinos and HNLs in our parametrization of Yukawas. Appendix B contains the derivation of the kinetic equations. Finally, in appendix D we list several sets of the model parameters along with the corresponding values of the BAU. These sets can be used by other groups as benchmarks to compare numerical results.

The νMSM
In this section we fix our notations by introducing the Lagrangian of the νMSM [7,8] and the parametrization of Yukawa couplings [55,56]. Even though these expressions are well known and have been presented many times, we list them to make the paper self-consistent.
The Lagrangian of the νMSM is the usual see-saw one [1][2][3][4][5][6] L = L SM + iν R I γ µ ∂ µ ν R I − F αILαΦ ν R I − M IJ 2ν c R I ν R J + h.c., (2.1) where L SM is the Lagrangian of the SM, ν R I are right-handed neutrinos labelled with the generation indices I, J = 1, 2, 3, F αI is the matrix of Yukawa couplings, L α are the lefthanded lepton doublets labelled with the flavour index α = e, µ, τ andΦ = iσ 2 Φ * , Φ is the Higgs doublet. We work in a basis where charged lepton Yukawa couplings and the Majorana mass term for the right-handed neutrinos M IJ are diagonal.
In the broken phase, the Higgs field acquires a temperature dependent vacuum expectation value Φ(T ) , which is 174.1 GeV at zero temperature. The Yukawa couplings in the Lagrangian (2.1) lead to the Dirac mass terms [M D ] αI = F αI Φ . The 6 × 6 symmetric mass matrix of neutrinos can be diagonalized by a complex orthogonal transformation. We will restrict ourselves to the see-saw limit |[M D ] αI | M I . In this limit the active neutrino flavour states are given by where U PMNS is the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix [57,58], ν i are mass eigenstates of active neutrinos, N I are the mass eigenstate of HNLs. The activesterile mixing matrix in the leading order of the see-saw mechanism is

JHEP07(2019)077
The parameters of the theory (2.1) are restricted by the see-saw mechanism since one has to reproduce the observed values of the mass differences and mixing angles for the active neutrinos [59]. 2 A convenient parametrization of the Yukawa couplings which automatically accounts for these observables was proposed by Casas and Ibarra in ref. [55]. The application of the Casas-Ibarra parametrization to the νMSM has been studied in ref. [56]. In the matrix form, the Yukawa couplings entering the Lagrangian (2.1) read (in the notations of refs. [36,38]) where m ν and m N are the diagonal mass matrices of the three active neutrinos and HNLs correspondingly. The matrix Ω is an arbitrary complex orthogonal N ν × N N matrix, where N ν is the number of left-handed neutrinos and N N is the number of right-handed neutrinos.
In the νMSM, the lightest HNL N 1 is the dark matter particle. A combination of Lyman-α and X-ray constraints puts strong bounds on the magnitude of its Yukawa couplings, see [60] and references therein. As a result, N 1 is almost decoupled and does not contribute to the see-saw masses of active neutrinos. Therefore, the masses and mixings of active neutrinos correspond to the case of two HNLs. In this case the matrix Ω can be chosen in the form with a complex mixing angle ω. The sign parameter is ξ = ±1. We fix it to be ξ = +1, since the change of the sign of ξ can be compensated by ω → −ω along with N 3 → −N 3 [61]. Throughout this work, we will use the abbreviations NH and IH to refer to the normal and inverted hierarchy of neutrino masses. In what follows it is convenient to introduce X ω = exp(Im ω). (2.7) In the case of two right-handed neutrinos, the PMNS matrix contains only two CPviolating phases, one Dirac δ and one Majorana η, see appendix A for the details of our parametrization of the PMNS matrix. Two Majorana masses in (2.1) can be parametrized by the common mass M and the Majorana mass difference ∆M . Note that the physical mass difference controlling the oscillations of two HNLs is a sum of ∆M and a term proportional to the product of Yukawa couplings with Φ . The expression for this mass difference can be found in ref. [15].  Table 1. Parameters of the theory: common mass; mass difference; Im ω; Re ω; Dirac and Majorana phases. In the second line we indicate the ranges of these parameters which were considered in this work.
We end up with six free parameters of the theory which are listed in table 1 along with their ranges considered in this work. The common mass M of HNLs is restricted to the [0.1−10] GeV interval. The smaller masses are in tension with the Big Bang nucleosynthesis (BBN) [62]. Heavier HNLs, which we do not consider in this work, deserve a separate study. The ranges of the Majorana mass difference ∆M and Im ω are determined by the a posteriori requirement of generating enough BAU. The real part of the complex angle ω plays a role of a phase, therefore it is enough to restrict it, along with the single Majorana phase, to the interval [0, 2π]. The range of the Dirac phase δ is somewhat more subtle since it was restricted in the recent global analysis of neutrino oscillation measurements. We will comment on this in section 8.
Note that the relation (2.4) is not an isomorphism, i.e. more than one set of the parameters lead to the same Yukawas F αI and therefore are physically equivalent. Still, the parametrization (2.4) spans all possible values of Yukawas compatible with the oscillation data. Dependence of the resulting BAU on the Yukawas (and parameters in table 1) is in general very complicated. Therefore a thorough study of the parameter space of the model is required in order to determine the boundary of the region where successful baryogenesis is possible in terms of the experimentally interesting quantities.

Experimentally observable quantities
The parameters listed in table 1 are useful for the theoretical understanding of the model, but the last four of them cannot be directly measured. In this section, we discuss the experimentally observable quantities and their relations to the parameters of the model.

The total mixing
The formula (2.2) establishes the basis for experimental searches of the HNLs. It shows that an amplitude of a process involving HNL N I is equal to the analogous amplitude involving active neutrino ν α multiplied by Θ αI .
In order to understand, how weakly HNLs are coupled to the SM in general, it is helpful to sum |Θ αI | 2 over flavours of active neutrinos and over I = 2, 3. This defines the total mixing where m 2,3 are masses of active neutrinos and the normal hierarchy (NH) of the active neutrino masses is assumed (the inverted hierarchy (IH) case can be obtained by replacing m 2 → m 1 , m 3 → m 2 in (3.1)). The total mixing (3.1) controls the amount of HNLs produced in an experiment and the lifetime of these HNLs.

Individual mixings
The total mixing (3.1) is useful to quantify interactions of the HNLs with the SM particles, however, it is not sufficient for determining sensitivity of experiments to the HNLs. Therefore we also consider individual, or flavoured, mixings.
To clarify the role of flavoured mixings, let us consider, e.g. the SHiP experiment [49,63]. This is the beam damp-type experiment. An intense proton beam from the SPS accelerator hits the target. The main detector consists of a large empty decay volume with calorimeters and trackers at the end. In the SHiP set-up, HNLs are supposed to be produced mostly in decays of heavy mesons and the observational signatures consist of boosted charged particles originating from a vertex in the empty volume. The production is proportional to the partial decay width of a heavy meson into an HNL Γ(H → N I α ), which is in turn proportional to the |Θ Iα | 2 . It is important to note that the HNL production channels with different accompanying leptons α are in principle distinguishable. In the SHiP, they could be discriminated if the mass of HNLs is close to upper bounds of kinematically allowed regions. Let us illustrate this in an example. Suppose that one observes a decay of an HNL with the mass exceeding m Bc − m µ to a muon. This means that the HNL was produced along with an electron in the process B c → N I e since the process B c → N I µ is kinematically forbidden.
The decay widths of HNLs are in turn proportional to |Θ Jβ | 2 . The channels with charged leptons β in the final state are also distinguishable. In the example above the product of |Θ Iα | 2 and |Θ Iβ | 2 is important.
Therefore, individual mixings are phenomenologically relevant. Notice also that for the mass difference range which we are studying here the characteristic oscillation length is several orders of magnitude smaller than the length SHiP shielding and fiducial volume. 3 Namely, the oscillation length 100 m coincides to ∼ 10 −7 eV physical mass difference. Therefore for the lepton α in the target and β in the detector one has to sum incoherently over I, J. So the total dependence on mixings is The individual mixings as functions of the parameters are given by where the combinations C ± α for the NH case are and for the IH case where U PMNS αi with i = 1, 2, 3 are the elements of the PMNS matrix (should not be confused with |U α | 2 ). JHEP07(2019)077 Figure 1. Within the white regions it is possible to reproduced the observed value of the BAU (black solid curves). The dashed and dotted curves demonstrate how large the possible theoretical uncertainties could be. Namely, the dashed curves correspond to the condition Y B ≥ 2 · Y obs B , whereas the dotted lines correspond to Y B ≥ Y obs B /2 accounting for the factor of 2 uncertainty in the computation of the BAU. The thin grey lines show the see-saw limit, i.e. it is impossible to obtain the correct masses of active neutrinos below these lines. The blue line shows the projected sensitivity of the SHiP experiment ref. [65] as presented in ref. [66]. Left panel : normal hierarchy, right panel : inverted hierarchy.

Cosmologically motivated values of mixings
In this section, we present our main results, namely, the values of the total and individual mixings of HNLs with active neutrinos for which the observed BAU can be explained by the νMSM. We describe how these results were obtained in a separate section 8.
The value of the BAU can be characterised in different ways. Throughout this work, we use the variable Y B = n B /s, where n B is the baryon number density (particles minus antiparticles) and s is the entropy density. The observed value is Y obs B = (8.81 ± 0.28) · 10 −11 [64]. For each set of the model parameters, we numerically find the value of Y B . We are interested in the regions of the parameter space where one can reproduce the observed value Y obs B . The regions of successful baryogenesis are shown in figure 1. In order to indicate how large can be the effect of the theoretical uncertainties in BAU computation, discussed in section 6, we show the borders of the regions where one can generate 2 · Y obs B , and Y obs B /2. The cosmologically favoured region of the parameter space is larger for light HNLs than it was previously recognized (see also the discussion in section 9, in particular, figure 7). The fact that successful baryogenesis is possible for quite large values of the mixings rises the question about the upper bounds of sensitivity of the direct detection experiments. To illustrate this point, we estimate the lifetime of an HNL using expressions for the decay rates of HNLs from ref. [67]. 4 For instance, let us consider an HNL with the mass M = 5 GeV and mixings close to the upper boundary in figure 1. For such an HNL the lifetime is of the order of 5 · 10 −9 s. Estimating the gamma factor to be 10 we see that the decay JHEP07(2019)077 length in the lab frame is less than 15 m. This implies that, e.g., in the SHiP experiment this HNL will decay well before the detector. Therefore it might be interesting to revisit the current experimental bounds on HNLs. Let us note in passing that there also exist bounds from the Big Bang nucleosynthesis [23,60,67]. The question of the derivation of such bounds has been addressed in details for HNLs with the mass below 140 MeV in ref. [62]. For heavier Majorana HNLs an accurate derivation is still missing.
Results for the individual mixings |U α | 2 and products |U α | · |U β | are presented in figures 2 and 3.

Open access datasets
In the previous section we have presented the boundaries of the regions where all observed BAU can be addresses within the νMSM in terms of various combinations of the mixings of HNLs and active neutrinos. However, the parameter space of the νMSM is not completely JHEP07(2019)077 determined by the plots presented above, and there are more hidden parameters. These parameters can be essential for experimental searches for different signatures, and, e.g. it is interesting to know branching ratios, such as N → π α determined by U e : U µ : U τ . For instance, what ratios of U e : U µ : U τ are possible for some point in the allowed region of figure 2 or 3? This information is crucial to determine the decay length and branching ratios of various detection channels. In order to fill in this gap, we publish several datasets [68].
• Upper and lower boundaries of the M − |U | 2 region where successful baryogenesis is possible as functions of the common mass of HNLs. Note that these lines correspond to the region where one can obtain Y B ≥ Y obs B .
• The dataset of models with successful baryogenesis. The value of the BAU for every model (a parameter set) in this list lies in the range [Y obs B /2, 2 · Y obs B ]. Note that the value of the Y B is recorded for each parameter set so one can easily perform another cuts. The models of this dataset can be used to perform detailed Monte
• The dataset of models leading to various values of the BAU. Even though not all of these models provide a correct value of the BAU, they can be used to compare different theoretical approaches.
• Selected benchmark points are gathered in appendix D.
6 Generation of the baryon asymmetry

Kinetic equations
In this section, we discuss the machinery of baryogenesis in the νMSM. We present the kinetic equations which form the basis of the numerical analysis of this paper. These equation possess the same generic structure as those used in refs. [8,15,23,[32][33][34]. However, several important improvements are incorporated. These are: (i) splitting of the rates to fermion number conserving and fermion number violating ones [36,37]; (ii) accounting for neutrality of the electroweak plasma (this requirement was added to kinetic equations in [24]) and (iii) non-instantaneous freeze-out of sphalerons studied in [38,40]. The rates entering the kinetic equations are updated using the recent results of ref. [69]. The only remaining source of possible uncertainties is the averaging procedure described below. The detailed derivation of the equations is presented in appendix B. Here we start from the system of kinetic equations, introduce an ansatz which allows us to integrate these equations over the momentum and show how the gradual freeze-out of the sphaleron processes can be accounted for. The subscripts 2 and 3 are inherited from the νMSM and used to distinguish two HNLs participating in the generation of the BAU.
We are interested in coherent oscillations of HNLs and their interactions with leptons. The HNLs N 2 and N 3 are Majorana fermions with two helicity states each. Helicities are used to distinguish particles from anti-particles. We assign positive fermion number to HNLs with positive helicity and vice versa. Distribution functions and correlations of two HNLs are combined into matrices of density ρ N (ρN for antiparticles). The kinetic equations for leptons are presented in terms of the densities of the ∆ α = L α − B/3, where L α are the lepton numbers and B is the total baryon number. These combinations are not affected by the fast sphaleron processes and change only due to interactions with HNLs, therefore their derivatives are equal to the derivatives of the lepton number densities n Lα . Here we present the equations determining the generation of asymmetries in terms of ∆ α . In the next subsection, these asymmetries are related to the BAU. The system of kinetic equations reads

JHEP07(2019)077
In (6.1) f ν = 1/ e k/T + 1 is the Fermi-Dirac distribution function of a massless neutrino. The effective Hamiltonian describing the coherent oscillations of HNLs is where E N = k 2 N + M 2 and σ 1 is the first Pauli matrix. 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. This region includes both symmetric and Higgs phases.
In the Higgs phase, the rates can be split -in terms of ref. [28] -into "direct" and "indirect" contributions. 5 The direct contributions correspond to the processes where the Higgs field actually propagates. These reactions are also present in the symmetric phase. The processes with the Higgs field replaced by its temperature dependent expectation value give rise to the indirect contributions. These indirect contributions are crucial at low temperatures. As we show below, they are also important for the baryogenesis.
In eqs. (6.5) the first and second terms represent the direct and indirect contributions respectively. Below we discuss the direct contributions, the neutrino dumping rates γ ν,(±) , and the neutrino potential in medium b entering eq. (6.5) through E ν = k − b. The derivation of the indirect contributions is presented in appendix B.

JHEP07(2019)077
The direct contributions to the effective Hamiltonians h ± come from the real part of HNL's self energy / Σ N (p) = / pα + / uβ. Namely, one has 6 At hight-temperature limit the function h + reproduces the standard Weldon correction T 2 /(8k) [70]. If one neglects α, which is numerically insignificant, the h − is suppressed compared to h + by a factor M 2 N /p 2 . Thus h − is very small in the symmetric phase and, at the same time, the indirect contribution dominates in the Higgs phase. In our numerical computations we use the real part of the HNL safe energy calculated in ref. [28].
We now move to the direct contributions γ direct ± . They originate from 1 ↔ 2, 2 ↔ 2, and 1 + n ↔ 2 + n processes. The latter two require proper resummations [18,37]. Fermion number conserving rate comes mainly from the 2 ↔ 2 scatterings . Another contribution to γ direct + comes from 1 + n ↔ 2 + n. Fermion number violating rate comes from the 1 + n ↔ 2 + n processes (note that 2 ↔ 2 scatterings do not contribute to γ direct − ). For our numerical analysis we use γ direct ± kindly provided by the authors of ref. [69]. It is important to stress that in ref. [69] the temperature of the electroweak crossover T c has been extracted from the one loop correction to the Higgs potential. Computed this way, T c 150 GeV [28], whereas the non-perturbative result is T NP c 160 GeV [71,72]. Since we are using the rates from ref. [69], the crossover temperature is set to T c 150 GeV. We also implement one-loop running of the couplings following the approach of ref. [73].
The last two ingredients entering (6.5) are the neutrino dumping rates γ ν,(±) and the neutrino potential in the medium b. The function b can be calculated following, e.g. refs. [74,75]. The neutrino damping rates are related to its self energy / In the temperature region of interest, the neutrino damping rates are dominated by 2 ↔ 2 process mediated by soft gauge bosons. The calculation of these rates requires proper resummations [28]. We use analytical results presented in ref. [69]. The dependence on the Yukawa coupling constants factorises out from the rates (6.2)-(6.4). It is convenient to introduce the matrix of Yukawa couplings h αI related to the matrix F αI defined in (2.1) as follows In terms of these couplings we have (6.9) 6 Note that the term proportional to α has been omitted in ref. [69] since it is subleading. Here we keep it for completeness. Note also that according to the formal power counting of ref. [69] the contribution h direct − has to be omitted as well.

JHEP07(2019)077
The relation between the number densities and the chemical potentials to leptons in eq. (6.1) 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. [28,38]. In the symmetric phase its diagonal elements are ω αα = 514/(237 T 2 ), while the off-diagonal ω αβ = 40/(237 T 2 ), α = β. Note that relation (6.10) should be modified for temperatures below the T sph 131.7 GeV at which the sphalerons decouple [76]. Full expressions of the susceptibility matrices can be found in ref. [38].
The set (6.1) is a system of coupled integro-differential equations for each momentum mode of HNLs. Numerical solution of this system is a very complicated task [19,40]. However, a certain ansatz could be made to simplify the system. Namely, let us assume that the momentum dependence of the distribution functions is the equilibrium one, is the Fermi-Dirac distribution of the massive HNLs. Then it is possible to integrate the kinetic equations over the momentum and obtain a set of ordinary differential equations. This procedure is the main source of the theoretical uncertainty. The error can be estimated by comparing solutions of the averaged equations with solutions of the full set (6.1). This has been done first in ref. [19]. Results of this work indicate that the error in the value of the BAU doesn't exceed 40%. Authors of the recent study [40] have also solved the full system. They have found that the accurate result differs by a factor of 1.5 from the benchmark of ref. [33]. However, note that the equations of ref. [40] include effects that haven't been accounted for in ref. [33]. We have also tested the benchmark points listed in [77] using the same neutrino oscillation data as in [40] and found a surprisingly good agreement for the most of the benchmark points. A file with the values of the BAU for all benchmark points could be found in [68]. In what follows we will be conservative and assume that the averaging can lead to a factor of two error.
Let us finally present the system of equations that we actually solve numerically. It is convenient to introduce the CP-even and CP-odd combinations ρ with the integrated rates defined as

JHEP07(2019)077
Equations (6.11) are formulated in a static universe. The expansion of the Universe can be accounted for by rewriting eqs. (6.11) in terms of the so-called yields Y X = n X /s, where n X is the number density of species X and s is the entropy density which is conserved in the co-moving volume. For the numerical computations we use s(T ) calculated in refs. [71,78].

Gradual freeze-out of sphalerons
The asymmetry generated in the lepton sector is communicated to the baryon sector by the sphaleron processes. As long as these processes are fast compared to the rate of the asymmetry generation, the following equilibrium relation holds [79,80] where Φ is the Higgs vacuum expectation value, which is equal to 174.1 GeV at zero temperature. However, as was demonstrated in ref. [38], a deviation from the equilibrium with respect to sphalerons happens at temperatures around 140 GeV, i.e., before the freezeout. It was also shown that the errors stemming from the usage of the equilibrium formula can exceed an order of magnitude. To overcome this problem, one can implement the method suggested in ref. [38]. Namely, one solves the kinetic equation for the baryon number [80,81] where for the three SM generations were Y B eq is given by eq. (6.13) and α Y ∆α is calculated from the main system (6.11). It is enough to solve eq (6.14) starting from T = 150 GeV. This finishes the presentation of the kinetic equations. To summarize, our equations incorporate all physical effects that are relevant for the range of HNL masses considered here. The only source of errors is the momentum averaging of the kinetic equations. Based on the previous studies [19,40], we conservatively estimate that these errors do not exceed a factor of two.

Physics of asymmetry generation
Before solving eqs. (6.11) numerically, we briefly discuss the physics of the asymmetry generation. There are several important temperature scales. One of them is the already mentioned sphaleron freeze-out temperature T sph 131.7 GeV [76]. The lepton asymmetry generated at temperatures lower than T sph does not affect the final value of Y B .

JHEP07(2019)077
Second scale is a temperature at which the HNLs enter thermal equilibrium, T in . This temperature can be determined from the condition Γ N (T in )/H(T in ) 1, where H(T ) is the Hubble rate and we take the largest eigenvalue of the matrix Γ N (T in ). The Hubble rate during radiation-dominated epoch is Another important scale is the so-called oscillation temperature, T osc . Coherent oscillations of the HNLs play the crucial role in the generation of the individual lepton asymmetries. In fact, they provide a CP-even phase, which, being combined with the CPodd phase from Yukawas, leads to the generation of the individual asymmetries, see, e.g. the discussion in ref. [15]. This mechanism becomes efficient around the first oscillation. The temperature at which the first oscillation takes place can be estimated as [8,13] T osc δM M M * where δM is the physical mass difference. This mass difference is the splitting between two eigenvalues of the effective Hamiltonian (6.2). Since the asymmetry generation is efficient at T around T osc , one can roughly estimate the lower bound on the δM by requiring T osc > T sph . Let us consider two cases: This regime is sometimes referred to as oscillatory. In this case the kinetic equations could be solved perturbatively [8,13] if also T osc > T sph . 7 Late thermalization implies that Yukawa couplings are relatively small. In terms of the parameters listed in table 1, this regime is realized if the value of |Im ω| is small.
In this case, the dampings of the generated asymmetries are efficient, so this scenario is referred to as strong wash-out regime. Yukawa couplings must be relatively large, this can be in agreement with the oscillation data if |Im ω| is large as well. Sizeable damping causes a wash-out of asymmetries before the freeze-out of sphaleron processes. However, the production of HNLs and interactions with left-handed leptons at high temperatures are also enhanced, so the asymmetry generation is more efficient.
In order to account for all relevant processes, the kinetic equations have to be solved numerically.

Numerical analysis of the kinetic equations
Equations (6.11) together with (6.13) or (6.14) allow one to determine the value of the BAU for each parameter set in table 1. For most values of the parameters, the equations (6.11) have to be solved numerically. In this section, we discuss the procedure of solving these equations and demonstrate how the improvements in the equations influence the results.

The procedure for numerical solution of the equations
First of all, it is necessary to determine initial conditions. Right after the inflation the baryon and lepton numbers of the Universe as well as number densities of HNLs are equal to zero. 8 Thus we start from the vanishing Y − (T 0 ) ≡ n − (T 0 )/s(T 0 ) and Y ∆α (T 0 ). According to the definition of ρ + , at initial stage Y + (T 0 ) = n eq (T 0 )/s(T 0 ), where n eq (T ) is an equilibrium number density of a fermion with mass M .
The appropriate initial temperature can be specified on physical grounds. As we have discussed, the asymmetry generation starts around the temperature of the first oscillation, T osc . We have checked numerically that no significant asymmetry is generated at the temperature 10 · T osc , i.e. much before the onset of oscillations. Therefore we take T 0 = 10 · T osc if 10 · T osc < 10 3 GeV, or T 0 = 10 3 GeV otherwise. Now, having set up the initial conditions, we can solve equations (6.11). It is convenient to implement them using z = log(M/T ) as a variable. Even though the problem is reduced to the solution of the set of 11 ordinary differential equations (ODE) by means of averaging (6.12), it still remains challenging. The reason is that significantly different time scales are present in the system (6.11). Appropriate stiff ODE solvers, such as LSODA [82], can handle our equations quite efficiently. However, the integration time can be reduced further. Notice that the effective Hamiltonian entering equations (6.1) can be decomposed as Therefore, we can move to the 'interaction picture' with respect to H 0 . After this transformation, the equations can be solved using a non-stiff method.
The final value of the BAU is founded by solving eq. (6.14). This ensures that the value of the BAU is not affected by the assumption of an instantaneous freeze-out of sphalerons.
In order to find an appropriate method we have implemented equations (6.11) in the Python programming language. The SciPy library [83] allows one to use several different ODE solvers. We have found that the most efficient (in terms of the number of calls of the r.h.s.) one for our purposes is the LSODE [82]. The equations were then coded in the Fortran 77/95 along with the native Fortran implementation of the LSODE [84]. Note that for the successful integration it is important to carefully tune the parameters of the solver (such as absolute and relative tolerances for each variable).
We have also implemented the same system of kinetic equations in Mathematica [85]. This allowed us to validate the results obtained using the Fortran code. However, since the computation of r.h.s. takes much longer, the overall computation time is also very large. The Fortran realization outperforms the Mathematica one by more than four orders of magnitude.
The whole computation of the BAU for a single set of the parameters takes from 0.05 sec in the oscillatory regime (approximately |Im ω| < 2) to 1.0 sec in the strong wash-out regime (approximately |Re ω| > 5.5). The very efficient numerical procedure allows performing a comprehensive study of the parameter space.  figure 4 indicates the zero-temperature Higgs contribution to the physical mass difference δM in the limit ∆M/M → 0. Below this line the Higgs contributions to the physical mass difference dominates, whereas above the line, the physical mass difference is mostly determined by the Majorana mass difference. This means that δM cannot be much smaller than ∆M in the region above the line and δM cannot be much smaller than the Higgs contribution below the line. Smaller values of δM -which are interesting, e.g. for studies of resolvable HNL oscillations at the SHiP experimentare only possible if there is a cancellation between ∆M and the Higgs contribution. This cancellation can happen only close to the blue line.
It is also important to understand the role of the improvements that we consider. We want to address the questions:  sphalerons. In order to answer the first question, we effectively switch off fermion number violating processes in our kinetic equations. It is possible because the fermion number conserving and fermion number violating processes are neatly separated in eqs. (6.2), (6.3) and (6.4). So we can set γ − = h − = 0 for the whole range of temperatures. In order to model the absence of the Higgs phase at temperatures down to 130 GeV we put Φ(T ) = 0 and also put all direct rates to zero at T < 150 GeV . We consider the NH case and two different values of the common mass, 1 GeV and 10 GeV. The phases are fixed to the values (7.1a). We present the results in figure 5. In order to see how accounting for the charge neutrality of plasma modifies the results we replace the susceptibility matrix in (6.10) by a diagonal one. In figure 6 we compare results with and without susceptibilities. One can see that the effect is quite sizeable. In the same figure we demonstrate the results with and without careful treatment of sphalerons.
Inspecting figures 5 and 6 one can arrive at the following conclusions.
• Fermion number violating rates. Accounting for the fermion number violation increases the Y B . See figure 5, green dashed lines.
• Broken phase. Equations without the fermion number violation solved under the assumption that the Higgs vacuum expectation value is zero at all temperatures above the T sph lead to larger amount of the Y B for heavy HNLs at large |Im ω|. See figure 5, blue dotted lines.
• Neutrality of plasma. Accounting for the neutrality of plasma by means of susceptibilities is important. The effect is stronger for lighter HNLs. See figure 6.
• Freeze-out of sphalerons. The boundary of the allowed region in figure 6 is not sensitive to the method of calculation of the BAU from the lepton asymmetry. In fact, if one is interested in the upper bounds on the mixings (large |Im ω|) the instantaneous freeze-out of sphalerons can be assumed. See figure 6. JHEP07(2019)077

Study of the parameter space
In this section, we describe how the study of the parameter space of the model have been performed. Our strategy is a direct sampling of the parameters defining the theory. In subsection 8.1 we fix specific values of the phases δ, η, Re ω which maximize the generated asymmetry. In subsection 8.2 we sample the whole 6 dimensional parameter space.

Total mixing
In order to set the bound on the value of the total mixing (3.1) we need to find, for each value of the common mass, the largest value of |U | 2 .
Since the value of |U | 2 for a given mass depends only on Im ω, one can marginalize over phases δ, η, Re ω and mass difference ∆M and solve an optimization problem for Re ω. The optimization problem consists of maximizing (or minimizing for negative values) Im ω subject to Y B Y obs B . Several comments are in order. The value Y B can be both positive or negative. If it is possible to obtain some value Y 1 B for some X ω and mass difference, then it is also possible to obtain −Y 1 B for the same X ω and ∆M provided that the phase parameters in the model can vary freely. In what follows we would always take the absolute value |Y B | of the computed BAU.
Next, it is important to clarify what does Y B Y obs B actually mean. The kinetic equations that we solve contain an inherent error stemming from the assumption of equilibrium momentum dependence of density matrices. In order to account for this theoretical uncertainty we impose the following condition Y obs B /2 < |Y B | < 2Y obs B . In practice, it is easier to maximize |Y B | for given values of Im ω and M . If the maximal value of |Y B | exceeds e.g. Y obs B , then it is also possible to generate a smaller value of asymmetry. One can iterate this procedure on a grid in Im ω and M space. Then by

JHEP07(2019)077
interpolating |Y B | as a function of Im ω for a given M and finding roots of the equation |Y B | = κY obs B , κ = 0.5, 1, 2, one can find the upper and lower bounds on |U | 2 . The case of κ = 2 corresponds to the conservative assumption that the averaging procedure amounts to a twice larger asymmetry compared to the accurate treatment. Authors of ref. [40] have solved the full system of equations for several parameter points. Their results indicate that the averaged equations rather tend to underestimate the value of BAU. This case is indicated by κ = 0.5.
Maximizing |Y B | with respect to ∆M, Re ω, δ, η is a resource demanding task. It can be significantly simplified in the strong wash-out regime, i.e. for large values |Im ω|. In this regime the value of the total asymmetry strongly depends on the difference in the damping rates of active neutrinos. For a given mass, X ω and mass difference the damping rates are controlled by Dirac and Majorana phases together with the real part of ω. Note that a set of phases that maximizes the difference among these damping rates (and, the total lepton asymmetry consequently) also minimizes (maximizes in the case of IH) |U e | 2 . We have used the following values 9 NH: δ = π, η = 3π/2, Re ω = π/4, (8.1a) These choices of phases maximize one of the individual mixings U α (see appendix A). For the NH case the phases (8.1a) maximize U µ , while for the IH case the phases (8.1b) maximize U e . Since the phases are fixed, we need to find only the value of ∆M which maximizes BAU at each point of the M − Im ω grid. The upper bounds in figure 1 were obtained using the method described above. Let us stress that the same upper bounds can be obtained by the random sampling described in the next subsection. We have checked that two methods agree with each other. The lower bounds on the |U | 2 obtained with the fixed phases (8.1) are not optimal, since the asymmetry is generated in the oscillatory regime. Therefore, the lower bounds in figure 1 are obtained by the direct sampling.

Individual mixings
The mixings |U α | 2 depend on δ and η through the elements of the PMNS matrix entering eqs. (3.3) or (3.4). Therefore it is no longer possible to solve a simple optimization problem, as was the case for |U | 2 . However, our numerical routines solve kinetic equations for different values of parameters very efficiently. This allows us to perform a scan of the full parameter space. The parameter space is sampled as follows. As was already mentioned, we are restricted to the discrete grid in the common mass M in the interval [0.1, 10.0] GeV. The rest of parameters we sample randomly, so that log 10 ∆M, Im ω, Re ω, δ, η are distributed uniformly in the intervals specified in table 1. 10 Note that according to eq. (A.21), the 9 The value δ = 0 for the IH case is incompatible with the recent 3σ bounds of the NuFit 3.2 analysis.
However, setting δ = 354 • doesn't change the results presented here. 10 Note that for the Dirac phase we have actually used the 3σ interval from the NuFit 3. uniform distribution in Im ω approximately coincides to a uniform distribution of |U α | 2 in log-scale. However, in order to obtain the upper bounds more accurately, we also perform a flat sampling in the X ω . After computing the value of BAU for each point we select the points according to the criterion |Y B | > Y obs B . In order to plot figures 2 and 3 we have generated 2 800 000 points for each hierarchy type and selected only those points for which |Y B | > Y obs B . We also have utilized these datasets to obtain the lower bounds and to cross check the upper bounds on the total mixing |U | 2 .

Comparison with other works
Baryogenesis in the νMSM has attracted a lot of attention of the community in recent years. The first scan of the parameter space was performed in refs [21,23]. Authors of ref. [24] have accounted for the neutrality of the electroweak plasma which leads to O(1) corrections to the final asymmetry.
More recently, scans of the parameter space were performed by two groups, see refs. [32,33]. The role of fermion number violating processes was clarified in refs [29,36,37,86]. Implications of a non-instantaneous freeze-out of sphalerons were addressed in refs. [38,40].
In what follows we list corresponding works.

L. Canetti, M. Drewes, T. Frossard and M. Shaposhnikov [21, 23]:
The first detailed study of the parameter space. Only the symmetric phase has been considered. Asymmetries in the leptonic sector were described by means of the chemical potentials, i.e. neutrality of the plasma has not been accounted for. The rates were underestimated by a factor of two (see table 2 below). In the scan of the parameter space the values of phases were fixed to non-optimal values. As a result, the allowed region of the parameter space is much smaller compared to what we have obtained in this work.

JHEP07(2019)077
M. Drewes, B. Garbrecht, D. Gueter and J. Klaric [32,34]: In ref. [32] only the symmetric phase has been considered. The kinetic equations were generalized to the broken phase in ref. [34]. The rescaling of the parameters that simplified computations has been suggested in ref. [32]. The relation between leptonic chemical potentials and number densities accounting for the neutrality of the plasma has been implemented. This relation is analogous to eq (6.10), however it is valid only at large temperatures. In high-temperature limit this relation agrees with eq (6.10).
Note the persistent disagreement between the damping rate of the active neutrinos in ref. [32] and in our work (see discussion below).

P. Hernández et al. [33]:
Only the symmetric phase has been considered. The neutrality of the plasma has been accounted for, however, apparently, the susceptibilities disagree with those in ref. [32] and with ours at high-temperature limit. 11 The approach to the study of the parameter space is different from what we use in this work. The parameter space has been sampled by means of the Markov Chain Monte Carlo (MCMC) with certain priors. The cosmologically allowed regions of the parameter space were presented as contours 90% of all generated points. This method resulted in regions which are much smaller compared to what we have obtained in this work.

S. Antusch et al. [39]:
The scan of the parameter space of heavy HNLs (M > 5 GeV). Fermion number violating processes have been accounted for in the symmetric phase. The parameter space has been sampled by means of the Markov Chain Monte Carlo (MCMC). However, the selection criteria is different from [33]. Namely, the models leading to were selected. This approximately corresponds to 0.68 · Y obs B < |Y B | < 1.32 · Y obs B . Let us emphasize that the uncertainties in the value of |Y B | are theoretical, whereas the experimental uncertainty, characterized by σ Y obs B is much smaller. This is the reason why throughout this work we consider a larger interval for |Y B |.
J. Ghiglieri and M. Laine [28,37,40,69]: There were no scans of the parameter space. However, a thorough derivation of all rates has been performed. The susceptibilities have been calculated accounting for 11 It is important to clarify that the matrix C αβ , entering eq. (2.20) from ref. [33] agrees with our matrix ω αβ from eq. (6.10) provided that in ref. [33] the symbol µα denotes the chemical potential to left-handed leptons. Note that our µα are the chemical potentials to all leptons of flavour α. Therefore, µ (ref. [33]) α = µ (our) α − µY /2, where µY is the chemical potential to the hypercharge. We thank Jacopo Ghiglieri for pointing this out. However, once the chemical potentials are eliminated from the kinetic equations by means of eqs. (6.10) and (2.20) from ref. [33], the equations are actually different. Namely, in the r.h.s. of the equations the terms proportional to β ω α,β Y∆ β will appear. The matrices multiplying Y∆ β are different in ref. [33] and in this work.

JHEP07(2019)077
the non-zero masses of the fermions in ref. [28]. The full non-averaged system has been solved for several benchmark points in ref. [40].
After the ArXiv version of this paper had been released, a new study [69] appeared. This work contains the most up-to-date determination of both fermion number conserving and violating rates in the whole temperature region relevant for baryogenesis. It was pointed out in ref. [69] that the 2kΓ k part of the γ ν(+) was missing in the ArXiv version of the present paper. In the current version of the paper we correct this point. We have also updated all rates entering the kinetic equations using the results of ref. [69].
T. Hambye and D. Teresi [29,86]: There were no scans of the parameter space. A role of fermion number violating Higgs decays has been discussed. The considerations of ref. [86] are limited to the Higgs decays and inverse decays. The rate of the fermion number conserving processes has been underestimated compared to the ones including 2 ↔ 2 scatterings. This can be seen, e.g. from figure 4 of ref. [37]. Therefore, a direct comparison between our study and ref. [86] is not straightforward.
It is important to note that the generic structure of kinetic equations is the same in all studies of the low-scale leptogenesis. Therefore it is possible to compare the rates in the kinetic equations independently of their derivation. In order to be able to compare refs [23,32,33,40], we compute the corresponding rates at temperature T ref = 10 3 GeV. At this temperature the rates are dominated by lepton number conserving processes.
The production rate of HNLs, the communication term of HNLs, the damping term of the lepton asymmetries and their communication term can be described as where h 2 is a symbolic representation of an appropriate product of Yukawa coupling constants for each term. The values of the coefficients C i in considered works are summarized in table 2. Note that since authors of ref. [40] treat the momentum dependence exactly in their numerical computations, we cannot compare their rates, however the hierarchy among the rates and their values at k = 3T are the same as ours. Leaving aside ref. [23], one can see that the values of the coefficient C 1 do agree with a reasonable precision. However, the values of the other coefficients differ roughly by a factor of two from work to work. In order to understand this difference, we numerically solve our kinetic equations with the rates multiplied by a constant coefficients κ a (a = 1, 2, 3, 4) as follows.

JHEP07(2019)077
Article  Table 2. The coefficients of the rates in considered works. For the cases 2, 3, and 4 the values above do not reproduce the kinetic equations in each works exactly, but allow us to understand the qualitative behaviour in each case. We demonstrate the time-evolution of asymmetries up to T = 160 GeV in figure 8. The qualitative picture of figure 8 agrees with the results presented in figure 7. 12 After the preprint of this paper had been released, we received a comment from the authors of refs. [32,34]. They found a missing factor of two in their calculations. Once this factor is corrected, the relative sizes of the coefficients Ci in the corresponding row of table 2 will approximately agree with these of ref [33]. Namely, the Case 4 will be realized.

JHEP07(2019)077
There is also an important comment regarding studies of the parameter space. In fact, each point in this space defines a theory. It is not clear at all what could be a prior probability in the space of theories. The problem is not entirely philosophical. This can be most easily seen comparing the first columns of sub-plots in figures 4 and 5 from ref. [33] with the diagonal sub-plots in our figures 3 and 2. Our allowed regions of the parameters space a much larger than the contours shown in ref. [33]. The reason for this difference is that the study of ref. [33] relied on a Bayesian analysis of the Markov Chain Monte Carlo (MCMC). This analysis assumes certain prior probabilities in the space of theories and depends strongly on the chosen priors [33]. We advocate the point of view that each parameter point leading to the correct values of the observables (such as neutrino mixing angles and the value of the BAU) should be accounted for.

Conclusions and outlook
In this work we have performed the thorough study of the parameter space of baryogenesis in the νMSM. All important effects have been accounted for in our kinetic equations. Our study improves that of previous works in several respects.
(i) The rates entering kinetic equations are calculated from the parameters of the theory.
In the symmetric phase, as one can see from the table 2, in ref. [23] the values of the rates were consistently underestimated. Moreover, apart from a factor of two difference in the damping rates, there is an agreement among all studies. Note also that all considered rates are practically the same in our work and in ref. [40].
(ii) In the broken phase the effects of the fermion number violation were systematically taken into account for the first time. These effects are important for the baryogenesis even though the temperature interval between the electroweak crossover and the sphaleron freeze-out is rather small.
(iii) We have accurately accounted for the sphaleron freeze-out utilizing the 'improved approach' of ref. [38].
(iv) Last but not the least improvement is related to the performance of the ODE solver which was used to solve the kinetic equations numerically. Impressive increase of efficiency of the numerical routine allowed us to perform a comprehensive sampling of the parameter space.
Our main results are upper and lower bounds of the region where successful baryogenesis in the νMSM is possible. We list them and stress significant points.
• Bounds in the |U | 2 − M plane, figure 1. The allowed region is significantly larger for light HNLs compared to the previous studies. Let us emphasize that the position of the upper bound is actually important for the direct detection. Even though this region seems to be the easiest for the direct detection owing to the most efficient production of HNLs, it might be actually not the case, because the life time of HNLs

JHEP07(2019)077
is short. HNLs can decay before they reach the detector. See the line of the SHiP experiment in figure 1. Also, it might be interesting to update the study of the neutrinoless double beta decay in the νMSM, refs. [30,31,33].
• Bounds on individual mixings |U α | · |U β | as functions of M . Note that we present also the off-diagonal elements. These are important for thorough simulations of the experimental sensitivity.
• The dataset of different choices of the parameters of the νMSM. This dataset can be used to compare our approach with other groups. As we have already stressed, we use momentum averaged kinetic equations. Computation of the BAU in the full system is highly non-trivial and a scan of the parameter space is very demanding. Therefore our parameter sets can be used as benchmark points to test different regimes of the BAU production with the accurate non-averaged equations. Models from the dataset could also be used by experimental collaborations for Monte Carlo simulations.

A Mixings of HNLs and active neutrinos
In this appendix we collect the formulae for the mixings of HNLs and active neutrinos. We considered the two-HNL case here. All formulae presented here are obtained for the normal hierarchy (NH) of the neutrino masses. The case of the inverted hierarchy (IH) can be obtained by the following replacement So, for example, m 2 + m 3 becomes m 1 + m 2 in the IH case. The Yukawa coupling constants entering the Lagrangian (2.1) can be decomposed using the Casas-Ibarra parametrization (2.4). Formula (2.4) can be rewritten as X ω = X ω e −iRe ω , (A.6) the Higgs vev at zero temperature is Φ(0) = 174.1 GeV. In the case of two HNLs the PMNS matrix contains two phases: Up to the leading order of the seesaw mechanism the mixing elements of HNLs are It is possible to show that the flavour components of the mixings are given by Using the unitarity of the PMNS matrix one can derive the following expression for the total mixing Corresponding formulae in the IH case can be obtained by means of the replacement (A.1).

B Derivation of the kinetic equations in Higgs phase
In this appendix we present the derivation of the kinetic equations in the Higgs phase. We consider only processes where the Higgs phase is substituted by its vev, i.e. only indirect processes (see, e.g. eq. (6.5)) which give the dominant contribution.
The kinetic equations accounting for the fermion number violating processes in the Higgs phase were derived in ref. [36]. A certain ansatz about the modification of the neutrino energies by the SM plasma has been made there. Here we extent the method of ref. [36] so that the interactions with the SM particles are consistently accounted for. This consideration is motivated by recent ref. [69] where an extra active neutrino rate, missing in the equations of ref. [36], is involved. Our derivation here confirms the results of ref. [69].
Below we derive the kinetic equations (6.1). First, we overview the idea behind the derivation and then present the actual calculations.

B.1 Overview of the procedure
• The lepton asymmetry in the νMSM is generated due to interactions of active neutrinos with HNLs and coherent oscillations of the latter. Therefore we need to derive the equations describing the evolution of number densities of both active neutrinos, ρ να and HNLs ρ N I as well as correlations of HNLs.
• We work in the Heisenberg picture and introduce a time-independent density matrix of the complete system ρ. The distribution function of a particle created by an operator a † is given then by Tr[a † a ρ].
• The time evolution of an operator O is governed by the Heisenberg equation where H is the Hamiltonian of the system. We are interested in the operators of the type O = a † a. So we need to derive the Hamiltonian in terms of creation and annihilation operators.
• The evolution equations for the operators describing the number densities of neutrinos and HNLs and correlations of HNLs involve some new operators. We write down the evolution equations for these new operators. These equations, in turn, involve new operators. We keep going until the system of the equations closes (note that this is very different from the Bogolyubov-Born-Green-Kirkwood-Yvon hierarchy which should be truncation at some level).

JHEP07(2019)077
• We obtain a set of a large number of first-order ordinary differential equations. Noticing that distinct time scales are present in these equations, one can integrate out fast oscillations and obtain a system describing the slow evolution of the quantities of interest (number densities and correlations).

B.2 Lagrangian and Hamiltonian
The Lagrangian in the mass basis (2.1) is useful for a study of the phenomenology of the νMSM. For the derivation of the kinetic equations it is more convenient to change the basis [15] so the Lagrangian reads (we use tildeÑ I to indicate the different basis) where ∆M is the Majorana mass difference so that the mass matrix of two heavier HNLs is M I = diag(M − ∆M, M + ∆M ). The matrix of Yukawa couplings h αI can be related to the matrix F αI defined in (2.1) as follows It is convenient to further rewrite the Lagrangian (B.2) by unifying two Majorana fields into one Dirac field Ψ = N c 2 + N 3 . After that, the Lagrangian in the Higgs phase reads We treat the mass difference of HNLs and their interactions with left-handed neutrinos as small perturbations. It is important that all the SM interactions -including these of active neutrinos -occur with much larger rates compared to those originating from the (B.6). It is therefore reasonable to formulate a perturbation theory in small parameters of (B.6). In what follows we realize this program.
The momentum expansion of the HNLs field is given by We orient p along the z axis. choose plane wave solutions u(p, s) and v(p, s) the helicity states of Ψ are s = ± and the operators a s and b s obey the usual anticommutation relations
In the HNL sector, we assign a positive fermion number to a particle with a positive helicity and to an antiparticle with a negative helicity. For instance, one HNL is created by a † + (±p) and another one is created by b † − (±p), see table 3. Attributed in this way, the fermion number is conserved in the limit M → 0, ∆M → 0.
We will work in the matrix of densities formalism inspired by ref. [87]. At the end of the day, we want to describe the distribution functions of HNLs and their coherent oscillations. Therefore we are interested in time evolution of bilinears of the type O = a † N (p, +1)a N (p, +1). The time evolution of such operators is governed by the Heisenberg equation where H is the total Hamiltonian of the system. This Hamiltonian can be decomposed as In the last expression, H 0 is the Hamiltonian of the free Dirac field Ψ and the free Weyl fields ν Lα ; H int describes quadratic interactions of the HNLs and left-handed neutrinos and H SM int is the Hamiltonian describing all SM interactions of ν Lα . In order to be able to use eq. (B.9) we express the Hamiltonian (B.10) in terms of creation and annihilation operators. For the first and the second terms in eq (B.10) it is a tedious but straightforward task. The SM part has to be treated separately.
Using (B.7) and analogous decomposition for the neutrino fields, one can find that the free Hamiltonian H 0 is given by where E N (p) and E να (p) are the energies of HNLs and active neutrinos. As the SM model effects are accounted for, the energy of the active neutrino must be replaced by a temperature dependent dispersion relation, see e.g. refs [74,75]. We will use the same symbol E να for both vacuum and thermal energies of neutrinos. Note, however, that E να in medium can deviate significantly from the vacuum value which is just |p| (the small neutrino masses can be safely neglected).

JHEP07(2019)077
The interaction Hamiltonian H int can be further decomposed into the Majorana and Yukawa parts coming from the first and the second parenthesis in (B.6) respectively The Majorana part of the interaction Hamiltonian is where p ≡ |p|.
The part of the interaction Hamiltonian coming from the Yukawa interactions reads The interactions of active neutrinos with the SM plasma are described by the following Hamiltonian.
where J is the SM current coupled to neutrino. It is important that J anticommutes with the creation and annihilation operators from table 3.

B.3 Commutators and averaging
Now we can use the Hamiltonian (B.10) to derive the evolution equations. However, first we need to identify the operators of interest. Note that the interaction Hamiltonians (B.14) and (B.16) contain both positive and negative p. The distribution functions of HNLs and active neutrino can be constructed from the operators 20) or from the operators with inverted sign of the spatial momentum In what follows we will call "slow" the operators whose time derivatives are proportional to the small parameters, such as Yukawas. All other operators are "fast".
where x i are the slow operators (B.22) (their derivatives are proportional to small parameter ), while variables y k are fast. All coefficients a, b, c, d are time dependent functions of order of unity and E k are combinations of energies of HNLs and active neutrinos of type E N +E να , E να −E ν β , etc. Note that there is no summation over k in (B.26b). For the sake of clarity we have assumed the following power counting, both Majorana mass differences and Yukawas times the Higgs vacuum expectation value are proportional to a small dimensional parameter ∆M, h · Φ ∝ . We will show now, that the fast oscillations can be averaged or -in the language of effective theories -integrated out in the way that the final system describes the slow evolution of operators.
Let us first consider the system at moment t.
With the fixed x i (t) = x i eqs. (B.26b) read iẏ k = −E k y k + l ( c kl x l + d kl y l ) .

(B.27)
We choose a time interval t ∈ [t, t + T ], such that 1/E T 1/ , where E is of order of E N or E να . The first inequality means that x i (t) does not change significantly on this time interval. We solve the system (B.27) on this interval and denoted the solution asỹ k (t, x).

JHEP07(2019)077
The second inequality allows us to exclude the fast oscillations fromỹ k (t, x) by means of the averaging y k (t, x) = 1 T In practice, we first integrate out the interactions with medium. This allows us to express the operators O medium m in the r.h.s. of (B.24) in terms of O fast k . We are working in the leading order in small parameters , as a result, it is only the first term in the r.h.s. of (B.24) which should be modified. This modification can be summarized as where Σ(k) is the active neutrino self energy and we have suppressed the flavour indices of the active neutrinos. The imaginary part of the self energy can be parametrized as [70] Im / Σ(k) = / k Γ k /2 + / u Γ u /2, (B.31) where u µ = (1, 0, 0, 0) is the four-velocity of the medium. Once the fast operators related to the medium are integrated out, we also get rid of the fast operators of the type O fast k . Eventually one gets the closed set of equations in terms of the slow variables only. These equations have the following generic form where we have separated the coefficients c j,l into two groups, c and c . Eventually, the terms with c and c will lead to processes with and without fermion number violation correspondingly, see eq. (6.5).

B.4 The final form of the equations
In order to obtain the system of kinetic equations in the matrix form we introduce the convenient notations of ref. [36], Using these notations we arrive at the following kinetic equations Expressions for the rates and the effective Hamiltonian are present in section 6 so we do not repeat them here. Notice that ρ N , ρN , ρ να and ρν α in eqs. (B.36) depend on momentum so there is fact a set of equations for each momentum mode. During the whole period of the asymmetry generation, the leptons are in thermal equilibrium and different momentum modes communicate to each other. Therefore, the appropriate variables for the r.h.s. of equations (B.36) are the chemical potentials. Subtracting (B.36b) from (B.36a) and introducing the Fermi-Dirac distribution function for massless neutrino f ν = 1/ e k/T + 1 we can rewrite (in the limit of small chemical potentials) equations (B.36) in the following form [38] i dn ∆α dt where ρ eq N is the equilibrium distribution function of HNLs. 13 Note that the r.h.s. of eq. (B.37a) is written in terms of the density of the ∆ α = L α − B/3, where L α are the lepton numbers and B is the total baryon number. These combinations are not affected by the fast sphaleron processes and changes only due to interactions with HNLs, therefore their derivatives are equal to the derivatives of the lepton number densities n Lα .

C Rates in the symmetric phase
In this section we describe how the rates entering eqs. (B.37) are defined. For completeness we also include the fermion number violating Higgs decays and inverse decays in the symmetric phase. 13 The equilibrium distribution function appears as a result of application of the detailed balance principle.

JHEP07(2019)077 D Benchmark points
In this section we present several parameter sets along with the corresponding values of the Y B . These sets can be used to compare the numerical results among different groups. The full datasets could be downloaded from [68]. Note that we use the latest global fit to neutrino data [88].