GeV-scale hot sterile neutrino oscillations: a numerical solution

The scenario of baryogenesis through GeV-scale sterile neutrino oscillations is governed by non-linear differential equations for the time evolution of a sterile neutrino density matrix and Standard Model lepton and baryon asymmetries. By employing up-to-date rate coefficients and a non-perturbatively estimated Chern-Simons diffusion rate, we present a numerical solution of this system, incorporating the full momentum and helicity dependences of the density matrix. The density matrix deviates significantly from kinetic equilibrium, with the IR modes equilibrating much faster than the UV modes. For equivalent input parameters, our final results differ moderately (~50%) from recent benchmarks in the literature. The possibility of producing an observable baryon asymmetry is nevertheless confirmed. We illustrate the dependence of the baryon asymmetry on the sterile neutrino mass splitting and on the CP-violating phase measurable in active neutrino oscillation experiments.


Introduction
Explaining the matter-antimatter asymmetry of the Universe through experimentally verifiable laws of nature remains one of the most important open issues for particle physics and cosmology. The scenario of baryogenesis through GeV-scale sterile neutrino oscillations has established itself as a nice framework in which concrete progress can be made on all aspects of this problem. The original idea was put forward in ref. [1], and a significant reformulation, constituting the current understanding of various parametric dependences, was provided by ref. [2]. Representative examples of recent refinements can be found in refs. [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18]. Among these, the present investigation can most easily be contrasted with ref. [12], whose benchmark point we adopt as a central test case for our numerical solution.
The present paper is a follow-up to ref. [19], in which rates and rate equations were derived for the behaviour of baryon and lepton asymmetries and the sterile neutrino density matrix at complete leading order in Standard Model couplings. The derivation generalized and built up on techniques developed in several previous works [20][21][22][23][24][25][26]. In particular it required a resummation of infrared sensitive 1 + n ↔ 2 + n scatterings as well as a computation of all 2 ↔ 2 contributions to sterile neutrino production rates and chemical and kinetic equilibration coefficients. These coefficients display a non-trivial momentum dependence, which in combination with the general structure of the rate equations leads to non-trivial momentum dependences of different components of the density matrix as well.
The parameter space of the (type-I seesaw) model in question has been nicely delineated in ref. [6]. In a so-called "scenario I", two sterile neutrinos are responsible for generating active neutrino mass differences, the observed baryon asymmetry, and a large lepton asymmetry. A third sterile neutrino constitutes keV scale dark matter, whose production is resonantly boosted by the above-mentioned large lepton asymmetry. In a broader "scenario II", the production of a large lepton asymmetry is not considered, but the focus is otherwise on the same two-flavour problem for active neutrino mass differences and baryon asymmetry. In the parametrically most relaxed "scenario III", three flavours of sterile neutrinos participate in the production of active neutrino mass differences and the baryon asymmetry. In a technical sense, our study corresponds to scenario II, which is minimal in the dimension of its parameter space. However, the same methods would also permit to address the more restrictive scenario I if the solutions for the lepton asymmetries were followed deep into the Higgs phase, and the more relaxed scenario III if a larger-dimensional density matrix were considered. We postpone these numerically more demanding investigations into future.
The structure of this paper is as follows. The basic equations from ref. [19], transcribed into an expanding cosmological background, are reviewed in 2. The most important terms, helpful for analytic understanding and numerical estimates, are identified in sec. 3. The main numerical challenge of the problem, namely that both "fast" and "slow" processes play a role, is tackled in sec. 4. Numerical solutions are presented in sec. 5, and we conclude in sec. 6. Appendix A reviews the definitions and some relevant properties of the rate coefficients Q, R, S from ref. [19], appendix B explains our parametrization of neutrino Yukawa couplings, and appendix C summarizes our treatment of the so-called sphaleron rate.

Review of basic equations
We start by rewriting and completing the set of rate equations derived in ref. [19], transcribing them from a flat to an expanding background. The expansion is characterized by a Hubble rate H = √ 8πe/( √ 3m Pl ), where e is the energy density and m Pl = 1.22 × 10 19 GeV is the Planck mass. The entropy density s and the speed of sound squared c 2 s = ∂p/∂e also appear, where p is the pressure. Yield parameters are defined as where the n i stand for various particle number asymmetries ("particles minus antiparticles"). The coefficient functions A, B, ... introduced in ref. [19] are rescaled as Denoting furthermore where T max is a maximal temperature, T min is a minimal temperature, k T is a co-moving momentum, and k is the momentum at T = T min , the evolution equation for lepton asymmetry of generation a minus one third of baryon asymmetry reads coupling a sterile neutrino of flavour I to an active lepton of generation a;μ a ≡ µ a /T andμ Y ≡ µ Y /T are rescaled lepton and hypercharge chemical potentials; 2 and unexplained notation is identical to that in ref. [19]. A way to fix the values of h Ia in terms of observable quantities is reviewed in appendix B. The evolution equations of the density matrices, integrated along co-moving momenta, read The coefficients describing real processes (particle creations and annihilations) are whereas the unitary part of the evolution is determined by a Hermitean Hamiltonian with , (2.13) (2.14) We have here chosen H 0 to be traceless (the trace part drops out in eq. (2.8)). A further rate equation concerns the time evolution of the baryon asymmetry, and requires a careful discussion. Let us denote the right-hand side of eq. (2.4) as a "force", F a . If we were to write equations separately for Y a and Y B , they would have the forms where F diff is the anomalous baryon plus lepton number violating rate. Going over to the usual variables Y a − Y B /3 and Y B+L ≡ a Y a + Y B , the rate equations become At high temperatures, where F diff ≫ a F a , the first term is sometimes omitted from eq. (2.18) (cf. e.g. ref. [19]). However, we want to solve the equations down to low temperatures, where F diff ≪ a F a , and then this term must be kept. It guarantees that the baryon yield stops evolving below the electroweak crossover: (2.19) Following the notation of ref. [19], the anomalous force term here reads (n G ≡ 3) whereμ B+L is a chemical potential associated with the baryon plus lepton asymmetry, and Γ diff is the Chern-Simons diffusion coefficient, whose T -dependence is reviewed in appendix C. The equations above depend on the chemical potentialsμ a ,μ Y andμ B+L . The first two can be obtained by going through chemical potentials associated with lepton minus baryon asymmetries,μ a , and throughμ B+L , via [19] µ a =μ a +μ B+L T , Here, up to corrections of O(α 1/2 w , α s ) [22,27],  This closes the set of rate equations. (The matrix appearing in eq. (2.23) is modified in the Higgs phase [28], but for our considerations at T > ∼ 130 GeV where the Higgs expectation value is parametrically v < ∼ gT , this amounts to a higher-order effect.)

Identification of the most important terms
In order to solve the equations of sec. 2 numerically, it is convenient to go over into the interaction picture. Moreover, in order to understand the structure of the solution, it is helpful to identify which of the many terms on the right-hand sides of the equations are the most important ones. The latter maneuver is not necessary for a numerical solution at early times, however it facilitates finding a simplified solution valid at late times (cf. sec. 4.3).
As a first step, focussing for concreteness on two generations, we rename the upper diagonal component of H 0 , i.e. ( H 0 ) 11 , as The essential term here is the vacuum mass difference M 2 1 − M 2 2 . Let U be a rapidly varying phase factor satisfying and denote Apart from depending on time, the off-diagonal components play another important role: they contain the complex phases responsible for CP violation. Therefore, they act as sources for lepton asymmetry. This suggests a way to simplify the rate equations. Indeed we can define diagonal and off-diagonal components not only for the coefficients, but also for the density matrix. In particular, the rate equation for the lepton asymmetry, eq. (2.4), obtains a form in which the contributions from the diagonal and off-diagonal components of the density matrix are nicely separated (we denote n F ≡ n F (k T )): Here the coefficients have been evaluated up to leading order in small chemical potentials. The terms proportional to chemical potentials are washout terms, the others are source terms. At early times, the solution is dominated by the source terms on the second row. Consider then the source terms for the density matrix. The key point is that the helicity asymmetry, parametrized by ρ − II , is odd in parity (P). Given that sterile neutrinos are their own antiparticles, it is even in charge conjugation (C). Therefore it is odd in CP, just like lepton asymmetries. Consequently ρ − II is as small as lepton asymmetries, and in general much smaller than the other components of the density matrix. Moreover, the off-diagonal components ρ ± 12 are much smaller than the diagonal components ρ + II , because both their initial values and their equilibrium values vanish. 3 To summarize, we can assume that the solution satisfies In order to write the evolution equations in this limit, it is helpful to compactify the notation somewhat, denoting Now, for the diagonal helicity-symmetric ρ + II , all terms on the right-hand side involvingμ a , ρ − II , or ρ ± 12 , are small. Therefore the evolution equation reads In the numerics, other terms are trivially included, and in general they do affect the final results on a few percent level, however eq. (3.9) is sufficient for a qualitative understanding.
As far as the rate equations for ρ ± 12 are concerned, we need to include washout contributions from ρ ± 12 itself, as well as source terms from the large ρ + II . The latter can contribute both through a unitary oscillation part (parametrized by Q 0 ) as well as through a decay/production part (parametrized by Q ± 12 ): The other components follow from ρ + 21 = ( ρ + 12 ) * and ρ − 21 = ( ρ − 12 ) * .
Finally, with the same notation, the diagonal helicity-antisymmetric ρ − II obeys (3.13) More terms are needed than before because there are many effects of similar (small) magnitude. In fact, there is a substantial cancellation in the two terms proportional to Q 0 , which has to be properly tracked in the numerical solution.

Outline
There is a specific challenge with the solution of the rate equations of secs. 2 and 3, namely that certain modes evolve much faster than others. Normally, fast evolutions should be "integrated out", so that in the actual dynamics only slow modes appear. However, a fast rate can be important if it leads to a new effect, absent from the purely slow evolution. This is the case with sterile neutrino oscillations, leading to CP violation, and with anomalous baryon plus lepton number violation, converting a part of the total lepton asymmetry into a baryon asymmetry. More precisely, both of these rates cross the Hubble rate during the period under consideration [2], and therefore play a crucial role. At very high temperatures, the sterile neutrino oscillation rate is much smaller than the Hubble rate. Then there is no time for CP violation to take place, and no lepton asymmetries get generated. Around a certain temperature, referred to as the oscillation temperature T osc (numerically T osc ∼ 10 4 . . . 10 5 GeV for the benchmarks considered here), the oscillation rate is similar to the Hubble rate, and individual lepton asymmetries get generated. Later on the oscillation rate is much faster than the Hubble rate: fast oscillations can be averaged over, and the evolution becomes "decoherent".
In contrast, the baryon plus lepton number violation rate starts by being much faster than the Hubble rate. Later on it rapidly switches off, at a temperature referred to as the sphaleron temperature T sph (numerically T sph ∼ 130 GeV). For T ≪ T sph , this rate is exponentially small, and baryon number becomes a conserved quantity.
In the remainder of this section we show how the fast modes, whose direct numerical integration is challenging, can be handled. The basic idea for their treatment is that we solve their equations of motion within a "static" background of the slow modes, which appear effectively as parameters in the solution. This solution is then inserted into the equation of motion of the slow modes. Thereby the rate equations of the slow modes get modified through "virtual" fast corrections. This is similar in spirit to the usual effective theory approach.

Anomalous baryon number violation
Consider first the anomalous baryon number violation rate, discussed in eqs. (2.15)-(2.20). Let us define Y eq B+L as the value of Y B+L at whichμ B+L from eq. (2.23), and consequently F diff from eq. (2.20), vanishes [29]: Then the evolution equation for Y B+L (cf. eq. (2.18)) can be rewritten as Assuming that a F a , Y eq B+L and γ vary slowly, the fast evolution determined by γ can be integrated exactly in a short time interval

Fast sterile neutrino oscillations
The second fast term originates from sterile neutrino oscillations, described by H fast , cf. eq. (3.1). For our benchmark parameter values, H fast ∼ 10 8 for k ∼ 3T at T ∼ T sph , and tracking the corresponding oscillations on par with the slow evolution poses a challenge. Let us, however, look at the fast evolution on its own, in a given background of the slow modes. Consider the form of U from eq. (3.2), viz.
This is not integrable because of the non-trivial x-dependence of H fast . Suppose, however, that we integrate only over a short period of time (i.e. small interval of x, so that (x−x 0 )∂ x H fast ≪ H fast ). Then we can expand H fast in slow variations, and to leading order use a constant H fast , With this form, the fast oscillatory dynamics of eqs. (3.10)-(3.13) can be integrated. Concretely, denoting the solution for ρ ± 12 from eqs. (3.10) and (3.11), multiplied by U 2 as is needed in eqs. (3.5), (3.12) and (3.13), reads . Let now . . . denote an average of the solution over one oscillation period centered around x =x. Then, to leading order in 1/ H fast , The constant part iF ± /(2 H fast ) emerges because the phase factor (U * ) 2 in eqs. (3.10) and (3.11) is compensated for by U 2 in eqs. (3.5), (3.12) and (3.13). This yields a source term for lepton asymmetries as we now show.
In the evolution equation for the lepton asymmetries, eq. (3.5), the integration over the spatial momenta eliminates the second term from eq. (4.9), 5 up to corrections of order 1/ H 2 fast . Therefore we can replace Inserting F ± from eqs. (4.6) and (4.7) yields In this equation all terms on the right-hand side evolve slowly, i.e. without U 2 or (U * ) 2 .
We note in passing that if we sum over a, and subsequently undo the helicity symmetrization/antisymmetrization of Q ± (cf. eq. (A.3)), then the numerator on the second row of eq. (4.11) becomes and similarly for the third row. This is suppressed by the helicity-conserving coefficients . Nevertheless a total lepton asymmetry is generated even in the massless limit, because individual lepton asymmetries are generated through the source terms in eq. (4.11), and the washout terms (proportional toμ a in eq. (4.11)) depend on a. "Decoherent" evolution equations, such as eq. (4.11), can also be obtained for the density matrix. If we carry out an average like in eq. (4.9) but for ρ ± 12 , a simple exercise shows that Therefore the average value of ρ ± 12 evolves slowly towards equilibrium, where ... ′ ≡ ∂x ... . Consider finally the source terms for ρ − II , from the second rows of eqs. (3.12) and (3.13). Given that in the end we only need the integrals over momenta of these components, the oscillatory terms from eq. (4.9) lead to corrections suppressed by 1/ H 2 fast and can again be omitted. Inserting the non-oscillatory parts from eq. (4.9) yields Making use of the definitions of Q ± from eq. (A.3) shows that This structure is proportional to the helicity-conserving coefficients Q (−) and therefore suppressed by M 1 M 2 /(g 2 T 2 ) [19]. In addition, eq. (4.17) vanishes if the dependence on flavour indices is symmetric; this is violated only by soft corrections of order (M 2 1 − M 2 2 )/(g 2 T 2 ) [19]. In total, we thus find a suppression ∼ M 1 M 2 (M 2 1 − M 2 2 )/(g 4 T 4 H fast ) in the source terms. Obviously, the method presented above can only be used for H fast ≫ 1. Empirically, we find that it works well if H fast > ∼ 10 3 . Note that H fast depends strongly on k T , cf. eq. (3.1), so smaller values of k T decohere earlier than large values. Therefore Y ′ a − Y ′ B /3 should in general get a contribution both from a decoherent small-k T domain according to (4.11) and from a coherent large-k T domain according to eq. (3.5). We have verified that after the implementation of this setup, our results are independent of the precise value of H fast at which we switch from the coherent to the decoherent evolution.

Outline
For the numerical solution, we start from initial conditions at which both the sterile neutrino density matrix and all lepton asymmetries vanish, at some temperature T max . For theoretical consistency, this temperature has to be so high that sterile neutrino oscillations have had no time to take place [2]. It can therefore be determined from eq. (3.1), by requiring H fast ≪ 1. The solution depends on the co-moving momentum k T . By writing k T = T osc ≡ κ T osc , evaluating thermodynamic functions at leading order in Standard Model couplings, and omitting the thermal mass corrections from eq. (3.1), we get where M ≡ (M 1 + M 2 )/2 and ∆M ≡ M 2 − M 1 . We choose in practice T max = 10 7 GeV as the initial temperature, and keep track of momenta κ > ∼ 0.01. 6 An important aspect of the problem concerns the dependence of the solution on ∆M . With increasing |∆M |, the value of H fast at the low temperature T sph increases, and therefore the efficiency of baryon asymmetry generation, which is suppressed by 1/ H fast at low T (cf. sec. 4.3), decreases. At the same time T osc increases according to eq. (5.1), so that there is a longer period between T osc and T sph for the process to take place [2].
Another important dependence originates from the momentum k T . According to eq. (3.1), the oscillations start earlier for the smallest values of k T , and at a given temperature are fastest for the small-k T modes. At the same time, the damping coefficients Q, R, S grow rapidly with decreasing k T (cf. appendix A). This implies that the small-k T modes both oscillate and equilibrate much faster than the large-k T modes.

Parameter choices
As a main benchmark point we consider a case marked with ⋆ in fig. 4 of ref. [12], which lies in the middle of the viable domain of "scenario II" (cf. sec. 1) according to the parameter scans performed in ref. [12]. In the notation of appendix B, the input parameters read Here δ is a Dirac-like CP-violating phase, and Im z and φ 1 are other complex phases which are not observable in active neutrino oscillation experiments but enter when sterile neutrinos are considered. In order to consider the same physical situation as in ref. [12], we have inverted the signs of the complex phases, for reasons explained below eq. (B.9 Here, to a good approximation, h 2a ≈ ih 1a . This leads to cancellations in neutrino mass formulae whereby active neutrino mass differences can be kept at their physical values despite largish neutrino Yukawa couplings. In the domain M I ≪ gT that we are interested in, and restricting to temperatures T > 100 GeV so that processes relevant for the "symmetric phase" dominate [26], the coefficients Q, R, S only display a powerlike mass dependence: helicity-flipping coefficients are mass-independent, whereas helicity-conserving coefficients are quadratic in masses. We have evaluated the coefficients according to ref. [19], inserting M = 1 GeV as an IR regulator where needed. As an example, for T = 4 × 10 4 GeV (cf. fig. 3) and k T = 3T , the coefficients read The equation of state is taken from ref. [30] (cf. also ref. [31]). 7 The evolution equations are 7 The non-trivial feature of this equation of state is that the heat capacity has a noticeable peak at around the electroweak crossover temperature T ≈ 160 GeV. As a result, the Universe spends more time in this temperature range, diluting extra energy density into expansion. Therefore there is relatively speaking more time for various production and equilibration processes to take place at around T ≈ 160 GeV. We note that in principle the effect of sterile neutrinos should also be included in the equation of state, however this results in corrections on the percent level and is furthermore very difficult to implement correctly, as it requires solving the Einstein equations simultaneously with the other ones. integrated from T max = 10 7 GeV down to T min = 100 GeV, where the Chern-Simons diffusion rate has switched off and no more baryon asymmetry is being produced.

Results
In fig. 1, the separate lepton minus baryon asymmetries Y a − Y B /3 are shown for the benchmark point of eqs. (5.2), (5.3), together with the corresponding full baryon and lepton asymmetries. In fig. 2, the integrals over components of the density matrix are illustrated, normalised to the entropy density. In fig. 3, the momentum dependence of the density matrix is shown at T ≈ 4 × 10 4 GeV, where the lepton asymmetries are being most efficiently produced (cf. fig. 1). All shapes differ significantly from the Fermi distribution, with in particular the IR modes of ρ + II having already reached equilibrium. 8 Remarkably, the total baryon asymmetry that we obtain with the parameter values of eqs. (5.2), (5.3) is Y B ≈ 1.3 × 10 −10 , i.e. ∼ 50% larger than the value 0.86 × 10 −10 in ref. [12]. In other words, the parameter scans carried out in ref. [12] could be somewhat conservative in their predictions for the viable domain.
It can be noted from fig. 1(right)  against each other to a good approximation at T > ∼ 120 GeV. This is because in the symmetric phase Y L-B + 2 I Y − II , sometimes called a fermion number, remains zero up to corrections suppressed by M 1 M 2 /(g 2 T 2 ) [19]. At lower temperatures the coefficient Q (−) kicks in (cf. appendix A and ref. [13]) and fermion number violation becomes rapidly visible.
Finally, in fig. 4, we illustrate the dependence of the final baryon asymmetry on the sterile neutrino mass splitting and on the CP-phase δ. The parameters have been fixed according to eqs. (5.2), (5.3), except that we now consider the less favourable normal hierarchy. It is seen how the value of δ is important for obtaining the correct sign of the baryon asymmetry, and how the magnitude of the baryon asymmetry is strongly affected by ∆M .

Conclusions
The purpose of this study has been to numerically integrate the evolution equations derived in ref. [19], in order to determine how the sterile neutrino density matrix and the lepton and baryon asymmetries evolved in the Early Universe. We find that the momentum dependence of the density matrix plays an important role in the solution, with the IR modes oscillating and equilibrating much faster than the UV modes. Therefore the shape of the density matrix differs substantially from the Fermi distribution at the time when leptogenesis is most efficient, cf. fig. 3. This effect was not included in an extensive recent parameter scan which otherwise employed similar rates and rate equations as our study [12]. 9 As a drastic illustration for the importance of the momentum dependence, we find that even very soft modes 0.01T < k < 0.1T can give a surprisingly large ∼ 5% contribution to the final baryon asymmetry. For our benchmark point the soft modes add up to the contribution from the hard modes. Understanding more precisely the physics of the IR modes may merit further study. 10 As another observation from tracking the momentum dependence, we note that even though single-k T modes experience oscillations (cf. fig. 3 for a snapshot of spectra), the lepton asymmetries, which contain an integral over all momenta, oscillate much less (cf. fig. 1), because different momentum modes add up incoherently.
As a benchmark point, taken from ref. [12], we focussed on two sterile neutrinos which are somewhat but not extremely degenerate in mass, cf. eqs. (5.2), (5.3). Then the production of lepton asymmetries is fastest at T osc ∼ 4 × 10 4 GeV, much above the temperature T sph ∼ 130 GeV at which sphaleron processes cease to be active. This parameter choice represents a typical case for the so-called "scenario II" outlined in sec. 1. For this benchmark point we find a baryon asymmetry ∼ 50% larger than the observed value (i.e. the result of ref. [12]). However it would be easy to re-adjust the baryon asymmetry to the observed value, by modestly changing the sterile neutrino mass splitting or CP-violating phases, cf. fig. 4.
In the more restrictive "scenario I", which aims to generate not only a baryon asymmetry but subsequently also much larger lepton asymmetries, it is natural to choose parameters leading to T osc ∼ T sph . This case has recently been studied in refs. [13,17], and we hope to apply our methods to that situation in the future. Another case meriting further scrutiny is 10 For M, m 2 φ /(4T ) ≪ k T ≪ T , where m φ is the thermal Higgs mass, the coefficient Q (+) grows as ∼ m 2 φ T /k 2 T , cf. eq. (A.7). Inserting into eq. (4.11), the contribution from small k T is ∼ dk T /k T (n F − ρ + II ). Therefore there is a logarithmic IR sensitivity, dynamically lifted if the small-k T part of ρ + II has equilibrated. The IR sensitivity is even stronger in the terms influenced by the helicity-conserving coefficients Q (−) , cf. eq. (A.8). the so-called symmetry protected scenario, ∆M/M → 0 and |Imz| → ∞ in the language of eq. (5.3), which leads to large neutrino Yukawa couplings and therefore to the best prospects for experimentally detecting sterile neutrinos (cf. ref. [18] for an overview). We end by remarking that our main results, including the rate coefficients Q, R, S used, are publicly available from the web site http://www.laine.itp.unibe.ch/leptogenesis/. We have also tabulated final results for many more benchmark points than discussed in this presentation. Examples of additional points are those included in the parameter scan illustrated in fig. 4, showing the dependence of the results on the sterile neutrino mass splitting and on the CP violating parameter measurable in active neutrino oscillation experiments. We would be happy to add further results on the web site, should readers provide us with their desired input parameters in the format of eqs. Appendix A. On the rate coefficients Q, R and S The coefficients Q, R and S that parametrize the rate equations of sec. 2 (cf. eqs. (2.5)-(2.7) and (2.9)-(2.12)) capture the processes relevant for sterile neutrino production, their kinetic and chemical equilibration, as well as lepton number washout. They can be defined by considering the Euclidean correlator whereφ = iσ 2 φ * is a Higgs doublet, ℓ a = (ν e) T a is a left-handed lepton doublet of generation a, and k n is a fermionic Matsubara frequency. The analytic continuation k n −iµ a → −i[k 0 +i0 + ] gives the retarded correlator Π R (K), whose imaginary part equals the spectral function ρ a (K). Taking matrix elements of ρ a (K) with on-shell spinors leads to the desired rate coefficients: a R are chiral projectors; and u kτ I is an on-shell spinor for sterile flavour I in the helicity state τ = ±. The lepton and hypercharge chemical potentials have been scaled with the temperature,μ a ≡ µ a /T andμ Y ≡ µ Y /T . The specific combinations playing a role in the main body of the text are obtained by symmetrizing or anti-symmetrizing the original coefficients with respect to helicity, and in some cases by symmetrizing them with respect to flavour indices: As an example, consider very high temperatures and Born level processes. As discussed in ref. [23], one has to omit the lepton thermal mass m ℓ from the Born computation since it is not a mass in the usual sense but a modification of the dispersion relation at large momenta. The dominant processes are Higgs decays and inverse decays [20], and we may write Here n B and n F are Bose and Fermi distributions, ǫ φ ≡ (p + k) 2 + m 2 φ , µ H ≡ µ Y /2, µ La ≡ µ a − µ Y /2, P ≡ (p, p) is the lepton momentum, and k ± ≡ (k 0 ± k)/2. It is straightforward to carry out the integral over p, leading to logarithms and dilogarithms.
Taking matrix elements according to eq. (A.2) and expanding in chemical potentials yields the coefficients Q, R and S. For transparent expressions, let us restrict to M ≪ k. Then, employing the functions the coefficients read We observe that the coefficients grow rapidly at small k but are then cut off at k ∼ m 2 φ /(4T ). At the same time, they overestimate the correct values at k > ∼ T , because they do not contain the lepton thermal mass m ℓ that restricts the phase space in that region.
In ref. [19], not only the 1 ↔ 2 processes but also 1 + n ↔ 2 + n and 2 ↔ 2 contributions to ρ a (K) were included. However, in order to complete this task, use was made of the "collinear" kinematic simplification m 2 φ /T, M, m ℓ , m φ ≪ k. Given that we observe the domain k < ∼ m φ to give a significant numerical contribution to lepton asymmetries, we need to extrapolate the coefficients to that domain. In order not to grossly overestimate their values, we replace the collinear 1 ↔ 2 contributions by eqs. (A.7)-(A.12) at small k for m φ > m ℓ . We furthermore apply an overall scaling factor to the small-k corrections, in order not to inadvertently change the sign of the resulting coefficients in a region where their determination is not trustworthy.
In the domain m ℓ < m φ , i.e. close to the electroweak crossover, the small-k region cannot be corrected as above. Particularly at 120 GeV < ∼ T < ∼ 140 GeV, there is a lot of structure but also some numerical uncertainty in the determination of Q (−) , R (−) and S (−) . At the same time, "indirect" contributions, i.e. oscillation from active neutrinos, become important at these temperatures. Adopting the notation and results of ref. [26], we have included them as δQ (−) = Im Π R | indirect /k, which indeed dominates over the direct contributions in the broken phase. 11 A similar correction is expected for the chemical potential dependence, parametrized in the symmetric phase by R (−) , S (−) , however this would require a comprehensive re-organization of the framework, because in the broken phase the dependence on chemical potentials is non-linear and because the gauge potential A 3 0 develops an expectation value in addition to the hypercharge gauge potential. We do not dwell on these issues further here, apart from noting that we have checked that in practice the broken phase values of R (−) and S (−) play very little role for our benchmark point.
Finally we remark that one of the 2 ↔ 2 contributions, namely scattering off soft Higgs bosons, was also observed to give an IR-sensitive contribution in ref. [19]. Its eq. (3.34) needs to be refined at k < ∼ m φ , as the energy conservation constraint δ(q 0 − k + (k − q) 2 + m 2 φ ) can only be satisfied for k > m φ in the range 0 < q 0 < k. Concretely, we now evaluate eq. (3.34) of ref. [19] as 12 rather than approximating the square brackets through ln(2k/m φ ) − 1 for all k.

Appendix B. Parametrization of neutrino Yukawa couplings
We provide here a self-contained exposition of the parametrization of neutrino Yukawa couplings, in order to be clear about our sign and phase conventions.

B.1. General discussion
Let us consider the leptonic sector of a Lagrangian including right-handed neutrinos. In order to be transparent about minus-signs, we employ Euclidean conventions here: Here ℓ L ≡ (ν L e L ) T ;φ ≡ iσ 2 φ * is a Higgs doublet; ν c R ≡ Cν T R denotes a charge-conjugated spinor; and M , h e and h ν are complex matrices with generation indices.
Given thatν c R M ν R =ν c R M T ν R , the mass matrix M is symmetric, M T = M . Through the so-called Takagi factorization (a special case of singular value decomposition), it can be written as M = V ∆ M V T , where V is unitary and ∆ M is a diagonal matrix with real non-negative entries. The matrices V and ∆ 2 M can be found by diagonalizing the Hermitean matrix M M † . Subsequently V can be eliminated through a unitary rotation of ν R . In the following we assume that this field redefinition has been carried out, and that therefore M = diag(M 1 , M 2 , M 3 ), where M I ≥ 0 are referred to as the Majorana masses.
The Yukawa matrix h e can also be assumed to be real and diagonal. Indeed, a biunitary transformation permits us to write it as h e = W † R ∆ h e W L , where W R,L are unitary matrices. There is no unique choice for W R,L , but possibilities can be found by diagonalizing the Hermitean matrices h e h † e and h † e h e , respectively. In the following, we assume that ℓ L has subsequently been rotated as ℓ L → W † L ℓ L and e R as e R → W † R e R , so that h e is diagonal, with real positive entries proportional to charged lepton masses.
After the field redefinitions of ν R and ℓ L , the matrix h ν is in general complex and nondiagonal. There are three free phases in W L which can be used to remove redundancies. Therefore, the total number of parameters introduced by N = 3 right-handed neutrinos is 18 (N from M I and 2N 2 − N from the complex matrix h ν with three unphysical phases projected away). Of these, 5 are currently known (two active neutrino mass differences and three mixing angles) and 2 are frequently considered accessible (absolute mass scale of active neutrino masses and "Dirac-like" CP-violating phase in the active neutrino mixing matrix). The remaining 11 can be chosen as the three 3 Majorana masses M I , 2 "Majorana-like" phases in the active neutrino mixing matrix (see below), and 3 complex angles related to the so-called R matrix of the Casas-Ibarra parametrization [34] (see below). Combinations of these can possibly be constrained by 0νββ and B-factory-type experiments.
As a next step, let us go to the Higgs vacuum, settingφ ≃ (v/ √ 2, 0) T where v ≃ 246 GeV. We denote where Y corresponds to the notation of ref. [12]. Then, from eq. (B.1) and recalling the transformation carried out with M , the mass terms in the neutrino sector read Inserting −1 = CC and noting that ν T R C =ν c R , we can writē Thereby Here M ν corresponds to the notation of ref. [35], representing a matrix multiplying (ν L ν c R ) T . The matrix M ν is symmetric and can again be represented via the Takagi factorization: where m ν and M h are real matrices containing the active and sterile neutrino masses, respectively. According to eq. (2.17) of ref. [35], in the seesaw limit we can write where U PMNS is the Pontecorvo-Maki-Nakagawa-Sakata matrix, and R is orthogonal.
As the final step, U PMNS can be rotated away from active neutrino masses through ν L → U PMNS ν L . It is then re-introduced into the non-diagonal parts of the weak interaction term, The relation in eq. (B.7) underlies the so-called Casas-Ibarra parametrization [34] and its generalization beyond the seesaw limit [35]. Specifically, combining eqs. (B.5)-(B.7), inspecting the upper right block, and expanding to leading order in 1/M , we obtain We note that eq. (2.5) of ref. [12] cites the left version for M D , so in comparisons with ref. [12] we need to flip the signs of complex phases, if we want to study the same physical situation.

B.2. Parametrization of U PMNS
We proceed to the parametrization of U PMNS , which fixes the neutrino Yukawa couplings according to eqs. (B.2) and (B.9). As mentioned above, 6 real parameters are needed: 4 "Dirac-like" parameters like for the Cabibbo-Kobayashi-Maskawa matrix, and two additional parameters, which can be chosen as "Majorana-like" phases. Ref. [12] writes where c ij ≡ cos θ ij and s ij ≡ sin θ ij . For the mass differences, we denote ∆m 2 ij ≡ m 2 i − m 2 j . Two cases are considered, normal hierarchy (NH) and inverted hierarchy (IH). According to ref. [

B.3. Specialization to two sterile generations
After the general discussion above, we now focus on a special case. In the so-called νMSM parameter corner (scenarios I and II in the language of sec. 1), one Majorana mass is very small (∼ keV), and the corresponding Yukawa couplings are tiny, so that the contribution from these Yukawas to active neutrino masses is vanishing. Following ref. [12], the small Majorana mass is denoted by M 3 , and we set (h ν ) 3a → 0. 13 Consequently, the smallest of the active neutrino masses necessarily vanishes. Therefore active neutrino masses are now fixed: where the projection operator is P NH = 0 1 0 0 0 1 , P IH = 1 0 0 0 1 0 . (B.20) 13 In the line of work reviewed in ref. [6], it is rather the Majorana generation I = 1 that is decoupled.