Leptogenesis from Oscillations of Heavy Neutrinos with Large Mixing Angles

The extension of the Standard Model by heavy right-handed neutrinos can simultaneously explain the observed neutrino masses via the seesaw mechanism and the baryon asymmetry of the Universe via leptogenesis. If the mass of the heavy neutrinos is below the electroweak scale, they may be found at LHCb, BELLE II, the proposed SHiP experiment or a future high-energy collider. In this mass range, the baryon asymmetry is generated via $CP$-violating oscillations of the heavy neutrinos during their production. We study the generation of the baryon asymmetry of the Universe in this scenario from first principles of non-equilibrium quantum field theory, including spectator processes and feedback effects. We eliminate several uncertainties from previous calculations and find that the baryon asymmetry of the Universe can be explained with larger heavy neutrino mixing angles, increasing the chance for an experimental discovery. For the limiting cases of fast and strongly overdamped oscillations of right-handed neutrinos, the generation of the baryon asymmetry can be calculated analytically up to corrections of order one.


Motivation
Over the past decades the Standard Model of particle physics (SM) has been established as a powerful theory explaining almost all phenomena that are observed in particle physics. Its full particle content has been discovered eventually, and its predictions to this end pass all precision tests [1]. Nevertheless, it is clear that the SM cannot be a complete theory of Nature. Any attempt to explain the observed neutrino flavour oscillations with the SM field content relies on nonrenormalisable interactions mediated by operators of mass dimension larger than four, which are generally associated with the existence of new heavy degrees of freedom that have been integrated out. Moreover, the SM fails to explain several problems in cosmology. These include the origin of the matter-antimatter asymmetry in the Universe that can be quantified by the baryon-to-photon ratio [2,3,4] η B = 6.1 × 10 −10 .
The addition of n s ≥ 2 right-handed (RH) gauge-singlet (sterile) neutrinos N i (i = 1 . . . n s ) can simultaneously explain the observed light neutrino masses via the seesaw mechanism [5,6,7,8,9,10] and the baryon asymmetry of the Universe (BAU) via leptogenesis [11]. 1 . The sterile neutrinos are connected with the SM solely through their Yukawa interactions Y with the SM lepton doublets a (a = e, µ, τ ) and the Higgs field φ. The Lagrangian of this model is given by where L SM is the SM Lagrangian. The spinors N i observe the Majorana condition N c i = N i , where the superscript c denotes charge conjugation. Besides, ε is the antisymmetric SU(2)-invariant tensor with ε 12 = 1. 2 The eigenvalues M i of M , which in good approximation equal the physical masses of the N i particles, introduce new mass scales in nature. The requirement to explain neutrino oscillation data does not fix the magnitudes of the masses M i , as oscillation experiments only constrain the combination An overview of the implications of different choices of M i for particle physics and cosmology is e.g. provided by Ref. [13]. The relation between the parameters in the Lagrangian (2) and neutrino oscillation data is given in Appendix A. The magnitude of the M i is often assumed to be much larger than the electroweak scale. However, values below the electroweak scale are phenomenologically very interesting because they may allow for an experimental discovery of the N i particles and to potentially unveil the mechanism of neutrino mass generation. Various experimental constraints on this low-scale seesaw scenario are summarised in Ref. [14] and references therein. In the present work, we focus on masses M i in the GeV range. Apart from some theoretical arguments [15,16,17,18], the study of this mass range is motivated by the possibility to test it experimentally. Heavy neutrinos with M i < 5 GeV can be searched for in meson decays at B and K factories [19,20,21] or fixed target experiments [22], including NA62 [23], the SHiP experiment proposed at CERN [24,25,26] or a similar setup at the DUNE beam at FNAL [27,28]. Larger masses are accessible at the LHC [29,30,31,32,33,34,35,36,37,38,39,40,41,42,43], either via vector boson fusion (M i > 500 GeV), s-channel exchange of W bosons (500GeV > M i > 80 GeV) or in real gauge boson decays (M i < 80 GeV), but the perspectives would be best at a high energy lepton collider ILC [29,44,19,45,35], FCC-ee [46,19,47,48] or the CEPC [19,49].
Since the N i are gauge singlets, they can interact with ordinary matter only via their quantum mechanical mixing with left-handed (LH) neutrinos that arises as a result of the Higgs mechanism and is the reason why the SM neutrinos become massive. This mixing can be quantified by the elements of the matrix Event rates in experiments are proportional to combinations of the parameters which determine the interaction strength of the heavy neutrino N i with leptons of flavour a. Here U N is a unitary matrix that diagonalises the heavy neutrino mass matrix. For convenience, we also introduce the parameter 1 The possibility that sterile neutrinos compose dark matter is discussed in detail in the review [12]. 2 Note that SU(2) group indices remain suppressed throughout this paper.
that quantifies the total mixing of a given heavy neutrino of flavour i as well as the quantity It is of interest to determine for which range of values of U 2 ai heavy neutrinos can simultaneously explain neutrino oscillation data and the BAU. In the present work, we improve on the network of equations that describes the generation of the BAU from GeV-scale sterile neutrinos and develop analytic approximations to the solutions for phenomenologically relevant limiting cases.

Leptogenesis Scenarios
Any mechanism that generates a non-zero BAU has to fulfil the three Sakharov conditions [50] of i) baryon number violation, ii) C and CP violation and iii) a deviation from thermal equilibrium. 3 Parity and baryon number are already sufficiently violated in the SM, the latter by weak sphalerons [51] at temperatures larger than T ws 130 GeV [52]. In the Lagrangian (2), CP is (in addition to the CP violation in the SM) violated by complex phases in the Yukawa coupling matrix Y and the mass matrix M . The non-equilibrium condition can be addressed by the heavy neutrinos N i in different ways. These can qualitatively be distinguished by the relative magnitude of different time scales, which we express through the variable z = T ref /T . Here T is the temperature of the primordial plasma and T ref some arbitrarily chosen reference temperature, which we set to T ref = T ws for convenience, such that sphalerons freeze out at z = 1. We assume that the radiation dominated era starts with a vanishing abundance of N i , which appears reasonable due to the smallness of their couplings Y [53]. The heavy neutrinos are produced in a flavour state that corresponds to an eigenvector of Y (interaction basis). Since Y and M are in general not diagonal in the same flavour basis, they start to undergo flavour oscillations at z = z osc . Their abundance reaches thermal equilibrium at z = z eq . They decouple (freeze out) from the plasma and subsequently decay at z = z dec . While this picture qualitatively holds for all parameter choices in the Lagrangian, the values of z osc , z eq and z dec greatly vary.
In the original leptogenesis scenario [11], the N i are superheavy (M i T ws ). In this case, their production, freezeout and decay all happen long before sphaleron freezeout (z osc < z eq < z dec 1). The final lepton asymmetry is produced in the CP -violating decay of N i particles and partly converted into a BAU by the sphalerons. The non-equilibrium condition is satisfied by the deviation of the N i distribution functions from their equilibrium values at z > z dec . This scenario and various modifications have been studied in the literature in great detail and are reviewed in Refs. [54,55,56]. For M i in the GeV range under consideration here, however, the smallness of the light neutrino masses (3) implies that the Yukawa couplings Y ia must be very small. In this case the N i production proceeds much more slowly, and the non-equilibrium condition is satisfied by the initial approach of their distribution functions to equilibrium prior to sphaleron freezeout at z = 1. This scenario is often referred to as leptogenesis from neutrino oscillations [57] because coherent oscillations of the heavy neutrinos during their production lead to CP -violating correlations between their mass eigenstates at z ∼ z osc . These are then transferred into matterantimatter asymmetries ∆ a = B − L a /3 in the individual SM flavours a = e, µ, τ when scatterings convert some of the N i back into SM leptons. Here L a are flavoured lepton asymmetries, which overdamped oscillatory M = 1 GeV Reω = 3π/4 ∆M 2 = 10 −6 M 2 ∆M 2 = 2 × 10 −5 M 2 δ = 3π/2 α 1 = 0 Imω = 4.71 Imω = 2.16 α 2 = −2π U 2 = 3.6 × 10 −7 U 2 = 2.2 × 10 −9 are kept in equilibrium with the baryon asymmetry B by sphaleron processes. Since the violation of total lepton number due to the Majorana masses is suppressed at T > T ws M i , the total lepton number remains small initially: |∆ a | | a ∆ a | 0. A total asymmetry a ∆ a = 0 is, however, generated because part of the asymmetries ∆ a are converted into helicity asymmetries in the Majorana fields N i by washout processes with an efficiency that depends on the different flavours a. 4 If the washout is completed before sphaleron freezeout, all asymmetries are erased. If the washout is incomplete at z = 1, then a baryon asymmetry B survives, as B is conserved for z > 1.
Based on the relation among the time scales z osc and z eq , which is controlled by the Yukawa couplings of the sterile neutrinos and their Majorana masses, we can distinguish between two regimes: • In the oscillatory regime oscillations occur much earlier than the equilibration (z osc z eq ) such that the charges ∆ a are mainly generated at early times during the first few oscillations. This requires weak damping rates and hence small Yukawa couplings in order to prevent the charges from being washed out too early. In turn, this setup allows for a perturbative analysis in the Yukawa couplings.
• In the overdamped regime the equilibration of at least one heavy neutrino happens before any full oscillation among the heavy neutrinos can be completed (z osc z eq ). This requires either some degree of mass degeneracy amongst the M i because the mass differences govern the oscillation time or anomalously large Yukawa couplings Y . Yet, for a successful generation of the BAU, we must have at least one sterile neutrino that does not fully equilibrate. This setup allows for an analytic approximation in terms of quasi-static solutions that are driven by the slow approach of one of the sterile flavour eigenstates toward equilibrium.
Within this work, we illustrate our analytic and numerical results through two parametric example points that are specified in Table 1.
We shall introduce two theoretical benchmark scenarios that roughly correspond to the two regimes. The naive seesaw corresponds to a situation in which the Yukawa couplings are of the order where m 2 atm is the larger of the two observed light neutrino mass splittings and m lightest is the unknown mass of the lightest neutrino. In this scenario, there are no cancellations in the seesaw relation (3). This leads to rather small mixing angles U 2 ai and makes it very difficult to find the heavy neutrinos in experimental searches. Larger mixing angles can be made consistent with the observed neutrino masses if there are cancellations in the seesaw relation (3). One way to motivate this is to promote B−L, which is accidentally conserved in the SM, to a fundamental symmetry that is slightly broken. This possibility is usually referred to as approximate lepton number conservation, as it implies that the violation of the total L at low energies is suppressed compared to the violation of individual lepton numbers L a . In this limit one finds that heavy neutrinos with Yukawa couplings much larger than suggested by the relation (8) must be organised in pairs of mass eigenstates N i and N j which in the limit of exact B − L conservation form a Dirac-spinor Ψ N = 2 −1/2 (N i + iN j ). This implies Moreover, the heavy neutrino mass basis (where M is diagonal) and interaction basis (where Y † Y is diagonal) are maximally misaligned in the flavours i and j. One of the interaction eigenstates does not couple to the SM at all, corresponding to a zero eigenvalue in Y , while the other one can have arbitrarily large Yukawa couplings without generating large neutrino masses or a rate of neutrinoless double β decay that is in conflict with present observational bounds.

Goals of this Work
The seesaw Lagrangian (2) contains 7n s − 3 new parameters, where n s is the number of sterile neutrinos. For five of these (two mass splittings and three light neutrino mixing angles) best fit values can be obtained from neutrino oscillation data [58], see Appendix A. In view of upcoming experimental searches, it is highly desirable to identify the range of the remaining parameters that allow to explain the BAU via leptogenesis from neutrino oscillations. This question has been addressed by a number of authors in the past [57,59,60,61,62,63,64,65,17,66,67,68,69,70,71]. The viable parameter space in the minimal model with n s = 2 has first been mapped in Refs. [61,62,63]. 5 The results of this analysis have been used to examine the physics case for the SHiP experiment [25] and the discovery potential of a future lepton collider [46]. More recent studies [69,70] suggest that the viable parameter region is smaller. In particular, the maximal values of U 2 i that are for given M i compatible with successful leptogenesis are smaller than claimed in Refs. [62,63], making an experimental discovery more difficult. With the present paper, we aim to clarify this question. For this purpose, we derive approximate analytic solutions for the time evolution of the asymmetries in the oscillatory and overdamped regimes. This is in contrast to the initial study in Refs. [62,63], which was entirely numerical. Analytic solutions for the oscillatory regime have previously been found in Refs. [59,64,69,70], but cannot be used to identify the maximal U 2 i compatible with leptogenesis because the N i oscillations tend to be overdamped when some of the U 2 ai are comparably large. We confirm numerically that our analytic solutions are accurate up to factors of order one in the regimes where they are applicable. We make use of the analytic understanding to identify the parameter region that leads to the largest possible U 2 ai that is consistent with successful leptogenesis. Within this region, we search for the maximal value of U 2 numerically. Compared to the previous numerical scan in Refs. [62,63], we apply the results of improved calculations of the thermal production and washout rates in the plasma [72,73,74,75], include spectator processes, use a more realistic model for the temperature dependence of the Higgs expectation value v(T ) and an updated result for the value of T ws .
The parameter space in the model with n s = 3 is considerably larger and has been studied only partially in the context of leptogenesis from neutrino oscillations to date [64,65,70,67,68]. In Ref. [64] it has been pointed out that in this scenario the generation of the BAU does not necessarily rely on a mass degeneracy amongst the M i , which is required in the case with n s = 2 [59] as well as for resonant leptogenesis from N i decays [76,77,78,79,80]. This results have been confirmed in Refs. [65,70,67,68]. It has also been pointed out that leptogenesis can be achieved for larger values of U 2 i for n s > 2 [65,70]. A complete parameter scan for the model with n s = 3 would be highly desirable, but is numerically challenging. Our analytic understanding in specific corners of the parameter space will be helpful in this context, as it allows to identify the relevant physical effects and time scales.
This paper is structured as follows: In Section 2 we present the evolution equations for both the sterile neutrinos and the SM asymmetries, and we discuss the qualitative behaviour of the solutions. In Sections 3 and 4 we derive analytic approximations to the solutions in the oscillatory and the overdamped regimes, respectively. Constraints on the active-sterile mixing are derived in Section 5. We discuss the implications of our results and conclude in Section 6. Technical details can be found in a number of appendices. In Appendix A, we summarise the parametrisation of the masses and couplings in the seesaw Lagrangian (2) that is employed in this paper. We also explain the phenomenological interesting case of scenarios with an approximate lepton number conservation that can lead to a large active-sterile mixing. Appendix B contains an extensive derivation of the kinetic equations for the sterile neutrinos based on first principles of non-equilibrium field theory, while in Appendix C the kinetic equations for the SM particles, that also include spectator effects, are reviewed more briefly. Finally, Appendix D contains some details on the oscillations of the sterile neutrinos that are omitted in the main text.

Evolution Equations
We need to describe the real-time evolution of the fields appearing in the seesaw Lagrangian (2) as well as of the spectator fields these couple to in the early Universe from the hot big bang down to T = T ws (or z = 1). Since quantum correlations of the different mass eigenstates of the heavy neutrinos are of crucial importance, there is an immediate need to go beyond a formulation in terms of Boltzmann equations for classical distribution functions. The evolution of sterile neutrinos in the early Universe is often described by density matrix equations [57,59,60,61,62,63,67,69,70] that can be motivated in analogy to the more detailed derivation for systems of SM neutrinos [81].
An alternative way to derive quantum kinetic equations and systematically include all quantum and thermodynamic effects from first principles is offered by the closed-time-path (CTP) formalism of non-equilibrium quantum field theory [82,83,84]. We describe this approach in Appendix A. The main advantage is that it allows to derive effective kinetic equations that hold at the desired level of accuracy from first principles in a series of controlled approximations. More specifically, overcounting issues as well as ambiguities related to the definition of asymptotic states in a dense plasma can be avoided, and necessary resummations of infrared enhanced rates at finite temperature are straightforward.

Charge and Number Densities
We can safely assume that the charged fields are maintained in kinetic equilibrium by gauge interactions such that we can describe these by chemical potentials, which are in linear approximation proportional to the comoving charge densities, We use a parametrisation where corresponds to a comoving temperature in an expanding Universe with Hubble parameter H. Here, m Pl = 1.22 × 10 19 GeV is the Planck mass and g = 106.75 the effective number of relativistic degrees of freedom. The physical temperature is given by T = a R /a, where a is the scale factor. The main quantity of interest is the baryon asymmetry of the Universe or, more precisely, the comoving density B of baryon number as a function of time. It is violated by sphaleron processes that are fast compared to the expansion rate for z < 1 and connect B to the comoving lepton number density L = a=e,µ,τ L a . The slowly evolving quantities relevant for leptogenesis are which are conserved by all SM interactions (including weak sphalerons). Here where q a and q Ra are the comoving lepton charge densities of flavour a stored within left and right chiral SM leptons, respectively, and g w = 2 accounts for the SU(2) doublet multiplicity. Among the SM degrees of freedom, only and φ directly interact with the sterile neutrinos. Nonetheless, the remaining degrees of freedom can also carry asymmetries and participate in chemical equilibration. They are referred to as spectator fields [85,86,87]. The main effect of the spectators is to hide a fraction of the asymmetries from the washout, which only acts on the L a . Taking account of these, one arrives at relations where the coefficients are derived in Appendix C.2. The Majorana fields N i strictly speaking cannot carry any lepton charges. However, at temperatures T M i , their helicity states effectively act as particles and antiparticles. We describe the N i by the deviation δn h of their number density from equilibrium, that is formally defined in Eq. (B.46). Here, h = ± denotes the sign of the helicity ± 1 2 , and δn h is matrix valued. In the flavour basis where M is diagonal, the diagonal elements are the number densities and the off-diagonal entries correspond to quantum correlations. This allows to define sterile charges in terms of the helicity-odd deviations of the occupation numbers from their equilibrium values, which is introduced more precisely in Appendix B. The Yukawa interactions Y violate individual lepton flavour numbers L a at order O[Y 2 ] (e.g. by light neutrino oscillations). The Majorana mass M also violates the total lepton number However, at temperatures T M i most particles are relativistic and spin flips are suppressed, such that the quantityL is approximately conserved (up to terms of order M 2 i /T 2 ). Since the N i start from initial conditions that are far from equilibrium, the assumption of kinetic equilibrium is not justified for them in principle. We briefly discuss the error introduced by the use of momentum averaged equations in Appendix D.2, see also Ref. [88].
In terms of these charge densities, we next write down the set of quantum kinetic equations used in our analysis. A detailed derivation for the evolution of the sterile neutrinos within the CTP framework is given in Appendix B, while a sketch of the derivation for the equations of SM charges is presented in Appendix C.

Evolution of Sterile Neutrino Densities
In terms of the variable z the time evolution of the number densities and flavour correlations of the sterile neutrinos is governed by the equation The flavour matrix H vac N can be interpreted as an effective Hamiltonian in vacuum, and H th N is the Hermitian part of the finite temperature correction. The contributions involving the matrix Γ N and the vectorΓ N are collision terms. Explicit expressions for these are derived in Appendix B, 23 and v(z) being the z-dependent Higgs field expectation value. As pointed out in the previous section, we make use of the freedom of choice of the reference temperature scale T ref to fix it as the temperature T ws of weak sphaleron freezeout. However, for the sake of generality we keep the notation T ref throughout this paper.
Before explicitly solving Eq. (19), we discuss the basic properties of the solutions. For this purpose we neglect the backreaction term withΓ N . The qualitative behaviour of the system is governed by the eigenvalues of H vac N and Γ N , which determine the time scales on which the sterile neutrinos oscillate and come into equilibrium. While H vac N is diagonal in the flavour basis where M is diagonal (mass basis), Γ N is diagonal in the same basis as Y Y † (interaction basis). The misalignment between the two leads to sterile neutrino oscillations. That means, particles are produced in the interaction basis and then oscillate due to the commutator involving H vac N . At sufficiently high temperatures the correction H th N due to thermal masses is larger than H vac N , but by itself cannot initiate oscillations because it is diagonal in the same basis as Γ N . For n s flavours of heavy neutrinos, there are of course n s relaxation times z eq and n s (n s − 1)/2 oscillation times z osc , all of which in general can be different. For a qualitative classification of the oscillatory and overdamped regimes it is useful to consider the largest eigenvalues of the matrices H vac N and Γ N . We use the norm || · || of a Hermitian matrix as the modulus of its largest eigenvalue. In case of Y * Y t it is, for instance, associated with the interaction eigenstate with the strongest coupling to the primordial plasma. The first oscillation involves the sterile neutrino mass states N i and N j with the largest mass splitting and occurs at a time such that z 3 osc ||H vac N || = O(1). The relaxation time scale at which a sterile neutrino interaction state comes into thermal equilibrium is given by such that z eq Γ = O(1) with γ av being the averaged relaxation rate (over temperature).
If the slowest oscillation time scale is shorter than the fastest relaxation time scale, then leptogenesis occurs in the oscillatory regime. In this case the heavy neutrinos undergo a large number of coherent oscillations before coming into equilibrium, which in terms of the variable z become increasingly rapid. The baryon asymmetry is most efficiently generated during the first few oscillations, before the oscillations become fast (compared to the rate of Hubble expansion), cf. Figure 1. There is a clear separation between the time z osc when the asymmetry gets generated and the time z eq when the N i come into equilibrium and the washout becomes efficient. This allows to treat these two processes independently. We discuss this regime in Section 3.
If, on the other hand, at least one heavy neutrino flavour eigenstate comes into equilibrium before a neutrino that is produced in this state has performed a complete flavour oscillation, then the oscillations are overdamped, cf. Figure 6. As we illustrate in Section 4, this allows for baryogenesis with larger Yukawa couplings and consequently also larger active-sterile mixing angles U 2 ai . In the scenario with n s = 2, the largest possible values of U 2 ai can be realised when the first oscillation happens rather late (z osc ∼ 1), as otherwise the washout tends to erase all asymmetries before sphaleron freezeout. As a result of the integration over a long time, the power counting in Y that allows to estimate the magnitude of the asymmetries in the oscillatory regime may not be applied, and the backreaction term involvingΓ N may not be neglected. Eqs. (21,22) allow to relate the mass difference to the Yukawa couplings in order to determine which regime a given parameter choice corresponds to: Figure. 2 schematically illustrates where the oscillatory and the overdamped regime are located in the M i − U 2 plane for various mass splittings. We also indicate the points from Table 1 that we use in our examples in order to illustrate the two parametric regimes. For n s > 2 the situation becomes more complicated because there are more oscillation and equilibration time scales, which can be ordered in various different ways. Moreover, the constraints on the relative size of the individual U 2 ai from neutrino oscillation data are much weaker and allow for a flavour asymmetric washout (while for n s = 2 there is not enough freedom in the unconstrained parameters in Eq. (A.1) to realise vastly different values of individual U 2 ai [89,90]).

Evolution of SM Charge Densities
The time evolution of the asymmetries ∆ a is governed by the equation A sketch of its derivation is presented in Appendix C. The first term on the right-hand side is the washout that is complementary to the damping rate for the sterile charges, while the second term is referred to as the source term It describes the generation of SM asymmetries in the presence of off-diagonal correlations of sterile neutrinos.

Oscillatory Regime
We now study the oscillatory regime, where the first oscillations of the off-diagonal correlations of the sterile neutrinos happen much earlier than their relaxation toward equilibrium, i.e. z osc z eq . The separation of scales z osc z eq allows to treat the generation of flavoured asymmetries from N i oscillations and their washout (which leads to B = 0) independently. At early times when z ∼ z osc , we can expand the solution to the coupled system of Eqs. (19,24) in the Yukawa couplings |Y * Y t |, as we specify within Section 3.1 in detail. At late times, when z ∼ z eq , the off-diagonal correlations have either decayed or their effect averages out due to the rapidity of their oscillations. Therefore, we can neglect the commutator term in Eq. (19) as well as the source term in Eq. (24) (i.e. the contributions explicitly depending on δn ij for i = j). This is done in Section 3.2. Our solutions hold for arbitrary n s as long as the slowest oscillation time scale is faster than the fastest equilibration time scale. Throughout this section, we work in the mass basis (where M is diagonal). In Figure 1, we present a characteristic example for the evolution of the particular charge densities for n s = 2.

Early Time Oscillations
We now identify in more detail the truncations that may be applied to Eqs. (19) and (24) when z ∼ z osc and solve the problem thus simplified analytically.
Oscillations of Sterile Neutrinos First, consider the thermal correction to the oscillation frequency of the sterile neutrinos due to thermal masses. While in the parametrisation of Eq. (19), the oscillation frequency induced by the vacuum term H vac N is growing with z 2 , the thermal contributions given by h th remain constant. As a result, at very early times, the thermal effects exceed the contributions from the vacuum masses. However, because H th N is generated by forward scatterings mediated by the Yukawa interactions, it is diagonal in the same flavour basis as Γ N , i.e. These act as a source for the generation of flavoured lepton asymmetries. We cut off the oscillations at the point when they become too rapid to make a significant contribution to the source term, as indicated in the plot. The middle panel shows the individual asymmetries generated in the three SM flavours. It is clearly visible that the total lepton asymmetry is only generated when the washout begins, and that its modulus remains smaller than that of the asymmetries in individual flavours at all times. The lowest panel shows the generated baryon asymmetry, where the green band indicates the error bars of the observed value.  The regions above/below the blue/red lines correspond to the overdamped/oscillatory regimes for the mass splittings indicated in the plot. The blue and red dots correspond to the two example parameter sets specified in Table 1. We can see that the blue point lies in the oscillatory and the red point in the overdamped regime.
the interaction basis in which heavy neutrinos are produced. H th N therefore commutes with δn h at early times (before H vac N becomes sizeable) and does not lead to oscillations. 6 For this reason, the thermal masses only lead to subdominant corrections in the oscillatory regime, and we neglect these in the following. In conjunction with the condition z eq 1, the same applies for the effect of the Higgs field through the term h EV as it only becomes important right before sphaleron freezeout, such that the relevant source term around z ∼ z osc remains unaffected. A more detailed discussion about these time scales is presented in Appendix D.1. The relation z osc z eq also leaves the backreaction mediated throughΓ in Eq. (19) as a higher order effect at early times z ∼ z osc , such that it only becomes important later, when the charges ∆ a have already been generated by the source term. In summary, for z ∼ z osc , and given the relation z osc z eq , Eq. (19) can be simplified to In order to compute q a as well as q N ii = 2n odd i we have to solve Eq. (26) both for helicity-even and helicity-odd distributions. The relation z osc z eq allows for a perturbative expansion in the coupling term |Y * Y t |. Solutions to order O(|Y * Y t | 0 ) are obtained when neglecting the right hand side of Eq. (26), what results in the diagonal terms 6 One may wonder whether the large thermal masses can lead to a big enhancement at z z osc by somehow amplifying a small population of the helicity-odd occupation numbers generated during the first fraction of an oscillation. However, it turns out that the main part of the charges ∆ a is produced well during the first full oscillation.
with the equilibrium solution (B.45), whereas the off-diagonal entries vanish. The first nonvanishing contribution to the charges ∆ a is O(|Y * Y t | 2 ), and it arises from the off-diagonal components of δn odd . These can be obtained by solving Eq. (26) with the replacement on the right hand side, such that we are left with solving with The general solutions to these equations are where C ij is an integration constant that in case of zero initial charge is determined to be and Sterile charges The helicity-odd off-diagonal elements δn odd ij are crucial for the generation of flavoured asymmetries q a . The diagonal elements (in the mass basis), on the other hand, can be interpreted as sterile charges q N , cf. Eq. (16). Within the present approximations, they vanish at z osc , when the flavoured asymmetries are generated. To show this, we solve Eq. (26) for diagonal, helicity-odd charge densities, where The solutions (31) lead to such that, when using the symmetry properties of the various tensors, F i (z) vanishes and consequently so does δn odd ij since we assume zero sterile charge as an initial condition. In total this results in which is valid at O(|Y * Y t | 2 ). In Appendix D.3 we show that for n s = 2 sterile neutrino flavours this even holds to all orders. However, in case of n s ≥ 3 flavours, already at O(|Y * Y t | 3 ) there appears a non-vanishing contribution that is however negligible in the oscillatory regime.
Asymmetries in Doublet Leptons and Sterile Neutrinos Likewise, in order to calculate the charge densities ∆ a in the oscillatory regime, we can neglect the washout term in Eq. (24) during the initial production process around z ∼ z osc . Since the generalised lepton number a q a + i q N i is conserved when T M i and we have previously shown that q N i 0 at z ∼ z osc in the oscillatory regime, we can conclude that B 0 and ∆ a −q a at z ∼ z osc . This immediately leads to the solution Now, when neglecting the washout that only becomes important at later times, we can obtain the flavoured lepton charge densities by substituting the source (25) into Eq. (38). To evaluate the resulting expression, we make use of the solutions (31) and integrate with the generalised hypergeometric function for p, q ∈ N 0 and w ∈ C, where Γ(x) is the Gamma function. Because soon after the first few oscillations the charges ∆ a saturate close to their maximal values ∆ sat a , cf. also Figure 1, we can use where the approximation holds for z moderately larger than z osc . On the other hand, as we have shown, the diagonal sterile charges q N i are negligible at early times [cf. Eq. (37)], so that the only asymmetries present in the plasma are flavoured asymmetries in the SM fields. To obtain these, we need the limit z → ∞ of Eq. (39) Putting these elements together and dividing by the comoving entropy density s = 2π 2 g a 3 R /45, we obtain In Figure 3 we compare the analytic results for δn odd 12 as well as for the late-time asymmetries (43) with the numerical solution. The discrepancies can be attributed to the fact that backreaction and washout effects are neglected so far. In a similar way as Figure 1, Figure 3 also illustrates the validity of the approximation in Eq. (38), where z is taken to infinity, because ∆ a indeed saturates after the first few oscillations.

Late Time Washout
At late times, when z ∼ z eq , we can neglect the oscillations of the sterile neutrinos because they have already decayed or they are so rapid that their effect averages out. In particular, there is no sizeable source for the flavoured asymmetries any more and also no other effects from off-diagonal correlations of the sterile neutrinos. This implies that the network of kinetic equations can be reduced to the following form where we use ∆ sat a and q N = 0 as initial conditions. Equation (44a) is easily obtained from Eq. (24) when dropping the source term. In order arrive at Eq. (44b), we keep the decay term Γ N as well as the backreaction termΓ N , while dropping the commutator in Eq. (19) and solve it for the helicity-odd, diagonal charges. This procedure is justified since the oscillations of the sterile charges around z = z eq are fast enough for their effect to average out. Note that the backreaction terms can be identified with the contributions involving q N i in Eq. (44a) as well as ∆ b and q φ in Eq. (44b). In Figure 4 the effect of the backreaction and spectator effects is presented where in particular the latter can have a substantial impact on the final result. The matrix A and the vector C appearing here specify the way how the spectator processes redistribute charges in the SM. Spectator processes have been neglected in most studies to date (except [67]), which corresponds to setting C = 0 and A = −1. The importance of including spectator effects is more pronounced than for conventional leptogenesis without flavour effects [87] because in the present scenario, the asymmetries are purely flavoured and the net result is due to an incomplete cancellation in the relation (C.15) that is rather sensitive to corrections in the individual terms.
Due to the hierarchy z osc z eq , we can use the charge densities generated through sterile oscillations around z ∼ z osc , cf. Eqs. (43) and (37), as initial conditions for solving the equations governing the washout process. For n s sterile flavours we can reduce Eqs. (44) to a linear first-order (upper panel), as well as the resulting time evolution of the three active charges (blue, solid), compared to their saturation limit given by Eq. (38) (red, dashed). The approximation (43) does not include washout effects since the washout time scale is assumed to be much later than the time scale of the oscillations. Furthermore, backreaction of the produced asymmetries on the N i evolution, as well as effects due to thermal masses and the Higgs expectation value are neglected. Note that the sum of the three charges ∆ a vanishes since lepton number violation only occurs at order |Y * Y t | 3 when washout effects are included.
where the components of the matrices K ∆∆ , K ∆N , K N ∆ and K N N read with i, j = 1, 2, . . . , n s sterile and a, b = 1, 2, 3 active flavours. Here A and C as defined in Eq. (15) account for the spectator processes. After diagonalising the Matrix K where T is a transformation matrix with the eigenvectors of K as column vectors, we are left with the solution with ∆ in = ∆ sat and q in N = 0 the asymmetries generated during the oscillation process at early times z ∼ z osc , cf. Eqs. (43) and (37). As the washout processes are suppressed during the initial creation of the asymmetries and because of relation z osc z eq , we can impose these initial conditions at A comparison of the evolution of the baryon asymmetry in the analytic treatment with the full numerical solution is shown in Figure 5.

Overdamped Regime
There are phenomenologically interesting parameter choices where the equilibration of one of the heavy neutrino interaction eigenstates happens before the first oscillation is completed, leading to an overdamped behaviour of the oscillations. This is particularly important in the case of mass-degenerate neutrinos, for which the first oscillation can happen at times as late as sphaleron freezeout, and in scenarios in which the Y ia are much larger than the naive seesaw expectation (8). Both of this can e.g. be motivated in scenarios with an approximate B − L conservation. In these scenarios one eigenvalue of Y † Y is always much smaller than the other, see Appendix A, so that one interaction eigenstate couples only very feebly to the plasma. Instead of being produced through direct scatterings, the feebly coupled state gets populated through oscillations with a sterile neutrino that has already equilibrated. Using the same perturbative approximation as in the oscillatory regime is no longer justified, because the larger decay rate cannot be treated as a small perturbation to the vacuum oscillation any more. Instead, we use a quasi-static approximation in a similar manner to applications to resonant leptogenesis from N i decay [91,92]. In the following we derive analytic expressions to treat the overdamped regime for n s = 2. Throughout this computation, we work in the interaction basis of the sterile neutrinos. An example plot for the generation of net baryon charge in the overdamped regime for two sterile flavours is shown in Figure 6.

Source of the Asymmetry
In the interaction basis, where Y Y † is diagonal, the fact that one interaction state decouples in the B − L conserving limit implies a |Y 1a | 2 b |Y 2b | 2 . We can therefore treat the smaller Yukawa coupling |Y 2a | as an expansion parameter throughout the following calculation. We will solve the equations for the positive helicity distribution δn +,ij , while all remaining distributions can be obtained through complex conjugation of the mass and Yukawa matrices.
The averaged sterile neutrino decay matrix Γ N inherits the flavour structure of the Yukawa matrices Y Y † . Therefore, in the interaction basis the decay rate Γ N as well as the thermal mass matrix H th N are both diagonal: We consider the regime where the equilibration of N 1 happens before the oscillations between the sterile flavours begin, which means that the rate at which δn 11 reaches its quasi-static value is much faster than the rate of the oscillations, Effects arising from the expectation value of the Higgs field can also be neglected up to a critical time z v ≈ 0.8 since h EV (z) = 0 for z 0.8. We separate the evolution equations into the directly damped equations, containing [Y Y † ] 11 , and the ones that are damped indirectly, through mixing with other sterile flavours, At this point we make the quasi-static approximation [91,92] to the solutions of Eqs. (53) by assuming that the interactions of the highly damped neutrino N 1 and its flavour correlations instantaneously reach values that are determined by the deviation of the feebly coupled state N 2 from equilibrium, i.e.
which allows us to express δn 11 , δn 12 , and δn 21 = δn * 12 in terms of δn 22 , where we have introduced Inserting these results into the equation for the weakly washed-out sterile neutrino N 2 yields the differential equation with the parameter Its absolute value introduces a new time scale The time scale |z c | indicates the instance when the vacuum part of the Hamiltonian z 2 H vac N becomes comparable to the thermal contribution H th N . The general solution to Eq. (57) is given by For times z |z c |, we can approximate this solution by which results in the equilibration time-scale for N 2 Therefore, unless |(H vac N ) 12 | 2 (H vac N ) 2 , N 2 will reach equilibrium before |z c |, justifying the usage of Eq. (61). Note that this situation naturally occurs in the pseudo-Dirac scenario, where the flavour and mass bases are maximally misaligned, such that (H vac N ) 11 = (H vac N ) 22 . The source of the lepton asymmetry is caused by the CP -odd correlation which yields the source term and is non-vanishing only at first order in the smaller eigenvalue |Y 2a |. The z dependence of the source term divided by T ref is shown in Figure 7. Note that the trace of the source a S a vanishes as we have (Y * Y ) 12 = 0 in the interaction basis.
Validity of the Approximations For times (Γ N ) −1 11 z |z c | 7 , Eq. (61) implies that dδn 22 /dz is small. Furthermore, we can approximate Hence, it is straightforward to see that the assumption made in Eq. (55) is justified in this regime, as the derivatives of δn 11 and δn 12 are much smaller than any of the individual terms on the right hand sides of Eq. (53),  (57,63). The quasi-static approximation still describes the solution correctly up to a rapidly oscillating transient correction. Taking the approximate form described in Appendix B, we can see that the contribution from the Higgs field expectation value quickly surpasses the contributions from the vacuum and thermal masses. Therefore it is sufficient to approximate the off-diagonal correlations by Where we assume δn 22 (z > z v ) to be approximately constant. For practical purposes, however, it is most convenient to completely neglect the CP -violating correlation, i.e. the source term (64), after the beginning of the electroweak phase transition.

Time Evolution of the SM Charges in the Overdamped Regime
At least one of the damping rates for the charges ∆ a is of the same order in |Y 1a | 2 as the larger of the sterile neutrino production rates. This implies that the washout of the active leptons typically happens at the same time as the overdamped oscillation of the sterile neutrinos. Neglecting the backreaction of the active flavours onto the sterile sector, as suitable for the oscillatory regime during the initial production of the asymmetries, is no longer an applicable approximation here. However, as all charges ∆ a are of first order in the smaller Yukawa coupling |Y 2a |, see Eq. (64), the calculation of the sterile charges at zeroth order in |Y 2a | remains unchanged. To correctly describe the evolution of the charge ∆ a , one has to solve the whole set of coupled differential equations at first order in |Y 2a |.
Suppression due to Backreaction To include effects coming from the backreaction of the active flavours onto the sterile sector, we consider once more the system of Eqs. (19,24). Among the CP -odd sterile distributions, the entry δn odd 11 receives the biggest correction due to backreaction. When neglecting the smaller Yukawa coupling |Y 2a |, the matricesΓ N take the form By applying the quasi-static approximation to the sterile neutrinos as in the previous section, we obtain the approximate densities of δn odd 11 = 2q N 1 , as well the off-diagonal correlations δn 12 . Inserting the quasi-static solutions back into the evolution equations of the SM leptons and the indirectly damped neutrino δn 22 gives with the effective washout matrix When we express the δn odd 22 dependence in Eq. (70) through the derivative dδn 22 /dz, the expression simplifies to To calculate the individual charges ∆ a , we can neglect the derivative dδn odd 22 /dz, as it is small for times z z c . The solution for ∆ a (z) can now be computed by integrating where w 1,2 are the two non-vanishing eigenvalues of the matrixW ab , and v bc the set of the corresponding eigenvectors. As a result of the conservation of the generalised lepton number (18), there is a vanishing eigenvalue. The generalised lepton number remains conserved when neglecting For practical purposes it is sufficient to completely neglect it for times before the equilibration of N 2 , z z eq N 2 , and to replace it by its quasi-static value for later times. By including corrections to δn odd 11 of order d∆ a /dz, and partially integrating the rate of change of the baryon asymmetry dB/dz, we can obtain the baryon asymmetry of the Universe up to an O(50%) error for z ≥ z eq N 2 . For the parametric example from Table 1, a comparison between this analytic approximation and the numerical result is shown in Figure 8.

Limits on the Heavy Neutrino Mixing
With Eqs. (49) and (75), we have found approximate analytic expressions for the BAU in the limiting cases that the oscillations of the sterile neutrinos occur either deeply in the oscillatory regime (the slowest oscillation time scale is much faster than the fastest equilibration time scale) or in the strongly overdamped regime (equilibration of one interaction eigenstate occurs long before the onset of the oscillations). From an experimental viewpoint it is interesting to identify the maximal values of U 2 i for which leptogenesis is possible. If one ignores constraints from direct searches for heavy neutrinos (see e.g. Ref. [14] and references therein for a recent summary), then these maximal values occur in the overdamped regime, which is characterised by a strong washout. There are two possibilities for preserving the baryon asymmetry at the electroweak scale from this At the original scale M

∆M 2
Im ω S(z) B(z = 1) Rescaled ξM η∆M 2 Im ω + log(η/ξ 3 )/6 S(η 1/3 z)ξη −1/3 B(η 1/3 )ξη −1/3 Table 2: Rescaling of the asymmetry washout. Either there is a strong hierarchy among the Yukawa couplings of heavy neutrinos to the different SM flavours e, µ, τ , causing one of the charges ∆ a to be approximately conserved, or the asymmetry is produced close to the electroweak scale, such that there is no time for a complete washout before sphalerons freeze out. In the case of n s = 2, a strong hierarchy among the doublet Yukawa couplings is not possible while being consistent with neutrino oscillation data. Therefore we need to resort to a strong mass degeneracy in order to delay the generation of the asymmetry until z 1. Yet, we are interested in maximising the mixing angles while keeping the washout rate of the SM flavours as small as possible, which constrains the parameters δ and α 2 . Minimizing the washout rate of the active flavours also introduces a difference between the normal and inverted hierarchies, as the minimal washout for the inverted hierarchy can be an order of magnitude smaller than the one for normal hierarchy given the same total mixing angle. Furthermore, maximising the analytic expression for the source also determines Re ω. Therefore it is only necessary to scan over the remaining three parameters: Im ω, M and ∆M .
When solving Eq. (19) for the helicity-even correlation function, we can use the fact that a solution with a rescaled time dependence δn even (zζ) corresponds to a solution of the same equation with the vacuum Hamiltonian replaced by H vac N → ζ 3 H vac N , the thermal mass by H th N → ζH th N , and the rate Γ N → ζΓ N . For parameter choices with large mixing angles, one of the eigenvalues of the decay rate of the sterile neutrinos is typically much larger than the other, (Γ N ) 11 (Γ N ) 22 , and the misalignment between the mass and flavour eigenstates is maximal, which implies that the only parameters playing a role in the evolution of the δn even correlation are the average Majorana mass M , the mass splitting ∆M 2 = M 2 1 − M 2 2 , and the imaginary angle in the Casas-Ibarra parametrisation Im ω. Any change of the mass scales M → ξM , or ∆M 2 → η∆M 2 , can therefore be compensated by a shift in Im ω → Im ω + log(η/ξ 3 )/6, as well as replacing δn even (z) → δn even (η 1/3 z). Note that although the oscillation and equilibration time scales change, their ratio remains the same.
To determine how the helicity-odd charges δn odd and ∆ a change under this parameter transformation, we first need to determine the change in the source term. In contrast to the decay rate where we can typically neglect the smaller Yukawa coupling |Y 2a | in the interaction basis, it is essential for the source term. By correctly applying the scaling transformation, the source term and with it the baryon asymmetry are rescaled according to Table 2. As a result, even if we do not achieve the observed BAU for some choice of parameters, by keeping the ratios ∆M 2 : |Y 1e | 2 : |Y 1µ | 2 : |Y 1τ | 2 constant, these transformation rules tell us how to find the parameters that lead to the desired result for the BAU just by changing the absolute mass and the mass splitting of the right handed neutrinos. Furthermore, by maximising B(η 1/3 )/η 1/3 , we can find the optimal mass splitting for producing the baryon asymmetry and then find the corresponding mass by determining ξ = B obs η 1/3 /B(η 1/3 ). For that mass these parameters give the maximal mixing consistent with leptogenesis. By using the scaling of the baryon asymmetry from Table 2, we find the maximal mixing angles consistent with baryogenesis for the mass range between 0.1 − 10 GeV as presented in Figure 9.
Note that including the continuous change of the expectation value of the Higgs field during the electroweak crossover reduces the final asymmetry as this effectively turns off the source at z = 0.8, while the washout continues until z = 1. We find numerically that this yields an error of Figure 9: The solid, dark blue lines show the largest and smallest value of U 2 we find to be consistent with neutrino oscillation data and the requirement to explain the observed BAU as a function of M = (M 1 +M 2 )/2. They are compared to the upper bound from direct search experiments summarised in Ref. [14] (solid black line), the lower bound from neutrino oscillation data (gray dashed "seesaw" line) and the lower bound from the requirement that the N i have a lifetime of less than 0.1s so that their decay does not modify primordial nucleosynthesis (dotted gray "BBN" line). The upper panel corresponds to normal neutrino mass hierarchy, the lower panel corresponds to inverted hierarchy. less than 50% for the parametric configurations with maximal mixing.

Discussion and Conclusion
In this work we study the production of lepton and baryon asymmetries from the oscillations of sterile neutrinos with GeV-scale masses in the minimal seesaw model. The main goal is to obtain an analytic understanding of the maximal heavy neutrino mixing angles U 2 ai consistent with the requirement to explain the observed BAU, while correctly accounting for backreaction and spectator effects. This is of crucial importance in order to assess the possibility of an experimental discovery of heavy neutrinos that may be responsible for the generation of light neutrino masses via the seesaw mechanism and for the BAU via low-scale leptogenesis. Baryogenesis via heavy neutrino oscillations can happen in different regimes, which can qualitatively be understood in terms of three time scales: the oscillation time z eq at which the first heavy neutrino flavour oscillation occurs, the equilibration time z eq at which the first heavy neutrino eigenstate comes into thermal equilibrium with the primordial plasma and the time z ws when weak sphalerons freeze out and baryon number becomes a conserved quantity, i.e., the BAU is frozen in. The generation of a baryon asymmetry can be understood analytically in the two extreme cases z osc z eq < z ws (oscillatory regime) and z ws > z osc z eq (overdamped regime). For heavy neutrino parameters that interpolate between these two regimes, we have to resort to solving the kinetic equations numerically.
In the oscillatory regime asymmetries in the individual lepton flavours are generated within the first few oscillations of the right-handed neutrinos at z z osc . At a much later time z z eq , the flavour asymmetric and lepton number violating washout generates a non-zero total lepton number from these, which is partly converted into a net baryon number by weak sphalerons. Once all heavy neutrinos come into equilibrium, all lepton asymmetries are washed out. However, if the washout is incomplete at z = z ws , a non-zero baryon asymmetry remains. The latter requirement implies that the Yukawa interactions of the sterile neutrinos must be sufficiently weak, and the analytic treatment is based on a perturbative expansion in the Yukawa couplings. To this end, our results agree with those previously found in the literature [59,64,69].
In the overdamped regime at least some of the heavy neutrino flavour eigenstates have Yukawa couplings that are much larger than the naive seesaw relation would suggest in absence of cancellations in the neutrino mass matrix. This can be made consistent with the smallness of the light neutrino masses if an approximate conservation of B − L is realised in Nature. This underlying symmetry implies that each strongly coupled heavy neutrino flavour eigenstate is accompanied by a feebly coupled eigenstate that completely decouples in the limit of exact B − L conservation. The two corresponding mass eigenstates form a Dirac spinor in that limit. We find an approximate analytic description in this regime by expanding in the tiny Yukawa coupling, and by employing a quasi-static approximation to the evolution of the strongly coupled flavour eigenstate, which comes into equilibrium before the flavour oscillations amongst the two can begin. In contrast to the oscillatory regime, the effect of the thermal masses and the backreaction of the produced lepton asymmetries on the heavy neutrino evolution cannot be neglected in this regime. Both of these tend to suppress the generated asymmetry. An additional suppression can be caused by the temperature-dependent Higgs expectation value if the asymmetry is generated near z = z ws . A complete washout of all asymmetries due to the large Yukawa couplings can be prevented in two different ways: Either one of the SM leptons couples to heavy neutrinos much more weakly than the two others (leading to a highly flavour asymmetric washout and a survival of the asymmetry stored in the weakly coupled SM flavour), or the heavy neutrinos have degenerate masses (in which case the oscillations and asymmetry generation occur very late at z ∼ z ws and there is no time for a complete washout before sphalerons freeze out). In the scenario with only two heavy neutrinos, a strong hierarchy amongst the couplings to different SM lepton flavours is ruled out by neutrino oscillation data, and leptogenesis can only be achieved with degenerate heavy neutrino masses. If there are more than two heavy neutrinos, then the extended parameter space allows to make a highly flavour asymmetric washout compatible with neutrino oscillation data, and baryogenesis is possible without a mass degeneracy [64,70].
The main new results of the present work are: • The equations of motion have been derived from first principles of quantum field theory in the CTP formalism. We have, for the first time, included the effects of thermal masses, backreaction from the generated asymmetries and spectator fields in this derivation.
• We have derived analytic approximations to the baryon asymmetry in case of both the oscillatory and the overdamped regime. While analytic solutions in the oscillatory regime have previously been found by several authors [59,64,69], the solutions in the overdamped regime are, to the best of our knowledge, presented here for the first time. Up to O(1) corrections they are consistent with numerical cross-checks.
• Based on these results, we have identified the largest possible heavy neutrino mixings consistent with leptogenesis in the scenario with two heavy neutrinos. Spectator effects, which account for the redistribution of SM charges due to fast SM interactions, and thermal masses have been included in both the analytic and the numerical treatment. While they have been neglected in recent studies so far, we have shown that they have a non-negligible impact on the final baryon asymmetry. Quantitatively we find that leptogenesis is possible for larger mixing angles than previously thought, which increases the chances of an experimental discovery.
In spite of this significant progress, several technical issues remain to be clarified in the future: • Our treatment relies on momentum-averaged kinetic equations. Since the assumption of kinetic equilibrium is not justified for the heavy neutrinos, this introduces an error of order one.
• Throughout this paper we have considered all SM Yukawa interaction to be in equilibrium, which is true for temperatures T 10 5 GeV, when the electron has finally equilibrated. However, the physical interesting regime, i.e. the time of the first oscillation, may already occur at higher temperatures.
• We assume the weak sphalerons to freeze out suddenly, which however is not completely true when electroweak symmetry is broken in a crossover, as it is the case for the SM. This could be phenomenologically important in the strongly overdamped regime, when the creation of the baryon asymmetry continues throughout the electroweak crossover.
• Our analytic solutions in the oscillatory regime are valid for an arbitrary number of heavy neutrinos. The treatment of the overdamped regime is, however, focused on the minimal realistic model, in which only two of these exist. A generalisation to the case with three or more heavy neutrinos, which includes a larger number of different oscillation and equilibration time scales, may be very helpful for efficient phenomenological studies.
To fully explore the discovery potential of present and future experiments, it would be highly desirable to perform a complete parameter scan of low-scale leptogenesis in the scenario with three heavy neutrinos. This should consistently include constraints from a wide range of past experiments that are sensitive to the existence of heavy neutrinos, in particular direct searches for these particles and indirect searches for lepton number or flavour violation. The present analytic results, in particular the new description of the overdamped regime, should also be applicable to assess possibilities of generating the BAU in extensions of the minimal seesaw model (2).

Acknowledgements
This research was supported by the DFG cluster of excellence 'Origin and Structure of the Universe' (www.universe-cluster.de).

A Parametrisation of the Seesaw Model and Neutrino Oscillation Data
The extension of the SM by n s sterile neutrinos introduces 7n s − 3 new physical parameters, i.e. 11 or 18 for the cases n s = 2 or n s = 3 considered in this paper. Various experimental constraints on these parameters are discussed in detail in Ref. [14]. The relation between the parameters in the Lagrangian (2) and constraints on the (presently incompletely determined [93]) light neutrino mixing matrix U ν , light neutrino mass matrix m ν can be expressed in term of the Casas-Ibarra parametrisation [94] The PMNS matrix can be factorised as with U ±δ = diag(e ∓iδ 1 /2 , 1, e ±iδ/2 ) and where the non-vanishing entries of the matrix V = V (23) V (13) V (12) are given by: The parameters θ ij are the mixing angles, δ is referred to as the Dirac phase and α 1,2 as Majorana phases. 8 The misalignment between sterile mass and interaction eigenstates is given by the complex orthogonal matrices R that fulfil RR T = 1. In case of three flavours it can be written as where the non-vanishing entries read with three complex angles ω ij , while for two flavours we have to deal with one complex angle ω and additionally a distinction between normal hierarchy (NO) and inverted hierarchy (IO) has to be applied: where ξ = ±1. In both cases Im(ω) determines the absolute size of the largest eigenvalue of the combination Y Y † . One can express the overall size of the mass eigenstates N 1 and N 2 defined in Eq. (7) as for normal hierarchy, for inverted hierarchy.
Finally, we shall make connection to the benchmark scenarios defined in Sec. 1.2. The naive seesaw is characterised by small values of Imω (or Imω ij ). In the approximately lepton number conserving scenario unitary transformations amongst the heavy neutrino fields can be used to bring Y and M into the form [95,96] Here a , a Y a and µ i M 3 ,M are lepton number violation (LNV) parameters, which must vanish if B − L is exactly conserved.M is the common mass of the two heavy neutrino mass eigenstates N 1 and N 2 that have comparable large mixing angles, the µ i quantify the mass splitting M 1 −M 2 . The deviation from maximal misalignment between the heavy neutrino mass basis (where M is diagonal) and interaction basis (where Y † Y is diagonal) in the flavours is quantified by the a . It is straightforward to see that U 2 a1 = U 2 a2 in the mass basis, i.e., both mass eigenstates couple with the same strength to SM leptons. The maximal misalignment implies that one interaction eigenstate has couplings of order Y a while the interactions of the other one are suppressed by the small parameters a , i.e., Y † Y has two eigenvalues of very different magnitude ∼ Y 2 a and ∼ 2 a . The analytic solution in Section 4 is effectively obtained in an expansion in a . In the parametrisation (A.1) the B − L conserving limit corresponds to large values of |Imω| 1. A third heavy neutrino (if it exists) must decouple in the B − L conserving limit, all its interactions are suppressed by a .

B Derivation of the Quantum Kinetic Equations
In this Appendix we provide a brief derivation of the quantum kinetic equation (19)based on first principles of non-equilibrium quantum field theory using the Schwinger-Keldysh CTP approach. For a pedagogical review of this topic see e.g. Refs. [97,98].

B.1 General Considerations and Definitions
We start our discussion assuming Minkowski background spacetime and generalise it to the radiation dominated Friedmann-Robertson-Walter metric in the subsequent Subsection.

Correlation Functions in a Medium
The use of S-matrix elements is not always suitable to describe non-equilibrium systems because there is no well-defined notion of asymptotic states, and the properties of quasiparticles in a medium may significantly differ from those of particles in vacuum. In contrast, observables can always be expressed in terms of correlation functions of the quantum fields, without reference to asymptotic states or free particles. There are two linearly independent two-point functions for each field. For a generic fermion Ψ these can be expressed in terms of the Wightman functions Here α and β are spinor indices, which we suppress in the following; flavour indices can be included equivalently. The . . . is to be understood in the sense of the usual quantum statistical average . . . = Tr( . . .) of a system characterised by a density operator . In the present context, we choose = eq SM ⊗ vac N , where eq SM is an equilibrium density operator for all SM fields and vac N is the vacuum density operator for sterile neutrinos. Physically this represents a situation in which the N i are absent initially and all SM fields have reached thermal equilibrium before the N i have been produced in significant amounts, which is justified by the smallness of the Yukawa coupling Y . The expressions (B.1) apply to both, Majorana fields (such as N i ) and Dirac fields (such as a ). The linear combinations have intuitive physical interpretations. The spectral function S A encodes the spectrum of quasiparticles in the plasma. The statistical propagator S + provides a measure for the occupation numbers. The correlators fulfil the symmetry relations If Ψ is a Majorana fermion, then there is an additional symmetry where C is the charge conjugation matrix and the transposition t acts on spinor as well as flavour indices.
It is often useful to introduce the retarded, advanced and Hermitian propagators, From this it follows that The usual Feynman propagator S F can be expressed as S F = S R + S < = S A + S > . Spectral, statistical, retarded, advanced and Hermitian self-energies / analogously, see e.g. [99,100] for a list of explicit definitions.

Equations of Motion
The correlation functions for quantum fields out of thermal equilibrium can be obtained from the Schwinger-Dyson equations which can be derived from two-particle irreducible effective action [101] in the CTP framework [84]. An explicit derivation is given in Ref. [98]. If the initial state at time t i is Gaussian (i.e. can entirely be specified by the initial conditions of the statistical propagators and one-point functions of all fields), then the above equations of motion are exact. Strictly speaking this is not true for (B.2) because eq SM is not Gaussian [102]. However, vac N is Gaussian, and we are primarily interested in the equation of motion for the heavy neutrinos.
The equations (B.8a) and (B.8b) can in principle be solved directly in position space [103,104,105,106,107,108,109], but it is often more practical to perform a Fourier transform in the relative coordinate x 1 − x 2 to Wigner space [110,111]. 9 This is the approach we take here. In order to perform the Wigner transformation, it is convenient to rewrite (B.8a) and (B.8b) with integration limits ±∞. For this purpose, we send t i → −∞ 10 and note that it can be seen that causality is maintained when substituting the retarded and advanced propagators and self energies by virtue 9 See also [112,113,114] for an alternative approach. 10 Boundary conditions at finite time can still be imposed by formally introducing singular external sources [109]. of the relations (B.6a) and (B.6b). By using Eqs. (B.6c) and (B.7) one finds S A,R = S H ± iS A . and From the kinetic equation (B.13b) it already becomes clear that H is the Hermitian part of an effective Hamiltonian that leads to oscillations of the sterile neutrinos, and G ≷ are dissipative gain and loss terms. N can be interpreted as a noise term that owes its existence to the fluctuationdissipation theorem. It is convenient to express whereH andḠ are H and G evaluated in local thermal equilibrium (with vanishing chemical potentials). The deviations δH and δG arise due to finite chemical potentials of the SM fields. 12 We now define the static solutionsS + = (S > +S < )/2 andS A = i(S > −S < )/2 as the solutions to the algebraic equations If the self energies / Σ are dominated by interactions with degrees of freedom that are in good approximation in equilibrium, then to leading order in the small couplings and gradients [103,108]. This yields The term ∂ tS + is due to the fact thatS + in the early Universe slowly changes due to Hubble expansion.

B.2 Quantum Kinetic Equations for Heavy Neutrinos
We now apply the general kinetic equation (B.19) to the case of heavy neutrinos. The associated correlators, self energies, effective Hamiltonians and rates will be attached with the subscript N .
It makes sense to split the self energies up into a partΣ computed in thermal equilibrium and for vanishing chemical potentials and a deviation due to non-zero chemical potentialsΣ, such that (B.20) In absence of chemical potentials, the self energies of the sterile neutrinos / Σ N and of the SM leptons / Σ can be factorised into a flavour dependent matrix of couplings and a reduced self energyΣ /, The absence of a superscript A , + , R , A or H indicates that the definition holds for either of these self energies. Here, g w = 2 accounts for the SU(2) multiplicity due to the SM doublets running in the loop. For non-zero chemical potentials, the reduced self energy also depends on the active flavour a. We can decompose Since gauge interactions keep all SM degrees of freedom in kinetic equilibrium, the deviations δH, δG ≷ and δN can be parametrised in terms of the chemical potentials, and the self energies thus fulfil the generalized Kubo-Martin-Schwinger (KMS) relationŝ where µ a and µ φ are the chemical potentials of the SM leptons and the Higgs. These chemical potentials are small at all times of interest, and we expand in µ/T to linear order. This yields a linearised KMS relation where in the second step we have suppressed the term that is quadratic in the chemical potentials. Note thatΣ From the definition of / Σ + N and the KMS relation (B.23) it is clear that δN N is quadratic in the chemical potentials, and we can neglect it. We also neglect the term δH N . In principle, it is of the same order as the term with δG N , but it only appears in a commutator, and δS N is approximately proportional to a unit matrix for T M i , Y 1 (δS N = −S N initially). This allows us to write The Fermi-Dirac distribution f F = f F (k 0 ) = 1/(e k 0 /T +1) arises from the KMS relation (B.25). We used the Leibniz rule for the term ∂ tS + N that can be expressed in terms of the stationary quantitȳ S A N and the distribution f F , where the derivative only acts on the latter one.H N andḠ N are the dispersive and dissipative part of an effective Hamiltonian. The termG N is responsible for the feedback or backreaction of the generated lepton asymmetry on the heavy neutrino dynamics.

Lorentz Decomposition and Off-Shell Kinetic Equation
We employ the decomposition of the non-equilibrium part δS N = iγ 0 δS N of the heavy neutrino propagator S N in Wigner space in Lorentz components used in [115,116], with the helicity projectors Note that we are using the Weyl (chiral) representation of the Dirac matrices. In the situation we consider here, they can all be expressed in terms of the functions g 0h , as shown explicitly in Refs. [115,116]. The relations can be found by taking traces over the products of the constraint equation (B.13a) with P h and different combinations of γ matrices. To linear order in M/k 0 and Y , they read 13 which is an approximation applicable in the regime M i T , Y ia 1. The equilibrium function iS + N can be decomposed in analogy to iδS N , where we replace g bh with the functionsḡ bh . In the self energies Σ N of heavy neutrinos and Σ of SM leptons we consider terms up to second order in Y . The kinetic equation for δg 0h can be obtained by taking the trace of Eq. (B.26) and inserting the relations (B.30a), On-Shell Kinetic Equation in an Expanding Universe The kinetic equation (B.31) for g 0h very much resembles an equation for density matrices commonly used in neutrino physics [81]. However, it is still an equation of motion for correlation functions (rather than particle numbers) because all quantities are defined for general k 0 that may also be off shell. The feeble strength of the Yukawa interactions implies that the narrow width approximation holds for the sterile neutrinoquasiparticles, and all phase space integrals are strongly dominated by the quasiparticle poles Ω i , which are defined by the poles of H −1 in the flavour basis where H is diagonal. In that basis we can approximateḡ In the ultrarelativistic regime T M i , we can further approximate δ(k 2 0 −Ω 2 i ) δ(k 2 ) in Eq. (B.33) because kinematically, the (thermal and vacuum) masses are negligible. We do, however, have to include these as a part of H N in the kinetic equation owing to their importance for flavour oscillations. In Eq. (B.33) we have used that iS + The above relations allow us to express the full equation (B.31) in terms of the equilibrium quantitiesḡ bh and the perturbative part δf bh as and there is no need to track particle and antiparticle numbers independently. For that reason we will focus on particles, hence we restrict to the case sign(k 0 ) = 1, and use Eq. (B.35) when needed. Using we eventually obtain an equation for on-shell distribution functions by integrating over k 0 , Note that we keep the same notation for Γ N ,Γ N and H N while these quantities are restricted to on-shell arguments k 0 = |k| in above equation. So far we have carried out our derivations in Minkowski spacetime. During the radiationdominated era, the expansion of the Universe can simply be included by using conformal time η instead of physical time [99], A prime denotes a derivative with respect to the conformal time η, and we additionally have made the flavour content explicit. During radiation domination, a = a R η, and when using the parametrisation introduced in Section 2, η = T /T ref . Note also that since a R /a = T can be interpreted as a comoving temperature, the equilibrium distribution for massless sterile neutrinos that appears in Eq. (B.38) is given by This distribution does not depend on conformal time, such that the term f eq = 0 in Eq. (B.38).
The effective Hamiltonian can be decomposed into a vacuum mass and a thermal mass term such that ( where M is the vacuum mass matrix [cf. Eq. (2)] such that Eq. (B.40a) is written in the most general form allowing for a complex symmetric mass matrix.Σ H N is the Hermitian part of the reduced self-energy of the sterile neutrino as defined in Eq. (B.21) [115] and v(z) is the temperature dependent expectation value of the Higgs field, 14 which we discuss in Section B. The anticommutator terms involve the reduced spectral self-energyΣ A N , and they are responsible for the decay of the deviations δf 0h toward equilibrium. Inserting the momentum and flavour structure one finds [115] The oscillations of sterile neutrinos induce flavour asymmetries in the active sector. The produced SM charges, i.e. those within the doublet leptons a of flavour a and the Higgs field φ, then lead to a backreaction effect that is described by the term Here, µ a and µ φ are chemical potentials for the doublet leptons and the Higgs boson, where we assume kinetic equilibrium for these species. We have linearised here in the chemical potentials, which is a valid approximation when µ a,φ T . It is convenient to define helicity-even and helicity-odd parts of the distribution functions In this work we assume that all lepton asymmetries remain small at all times. This allows to perform expansions to linear order in the µ a and δf odd . We cannot expand in δf even because the initial state (B.2) implies that at initial time, where the second equality holds on shell and at T M i , M j .
Rate Equations for Number Densities Though Eq. (B.38) has the same form as a density matrix equation for (quasi)particle occupation numbers, it is an equation of motion for the propagator (S N ) ij (which can be expressed in terms of the distribution functions δf 0hij ). In particular, it is valid for off-shell values of k 0 and holds for each momentum mode k individually. When accounting for backreaction effects, there will also be a coupling among the modes viaG. This is a considerable complication, and resolving the full momentum dependence would be a road block toward the goal of finding simple analytic approximations as well as fast numerical solutions. We therefore follow the common procedure [57,59,63] of reducing the problem to number densities in flavour of distribution functions by averaging the rates over the momentum. As we discuss below, this leads to order one uncertainties in the final result that should be resolved in future work. Some progress in this direction has been made in Ref. [88]. The developments that we present here may be helpful in order to address this issue. In order to cast Eq. (B.38) into a relation for the number densities of the sterile neutrinos, we perform an integration over momentum space. We are lead to introduce the equilibrium number density and the deviations The number densities δn even and δn odd are then defined analogously based on the distribution functions (B.43). We face the usual problem of approximating the momentum integral over products on the right-hand side of Eq. (B.38) by products of momentum integrals. Under the integral, the distribution functions are in general multiplied by different powers of the momentum. Inspection of the individual terms in Eq. (B.38) (that we discuss explicitly below) reveals that there are factors independent of k as well as factors of 1/|k|. In order to account for the latter, we replace 1/|k| by its average value .
For T M i , the spectral self-energy of the sterile neutrinos, that appears in G andG, is dominated by the t-channel exchange of a doublet lepton in association with the radiation of a gauge boson [73,117,68]. We follow Ref. [64], where the momentum averaging is applied through the replacement with the averaged relaxation rate γ av ≡ Γ av /T . This rate has been computed in different regimes by various authors [72,118,119,120,116,115,121,122,73,74,123,124,125,123,75,126].
Here we use γ av = 0.012, corresponding to the value from Ref. [68] based on Refs. [73,117]. In the backreaction term (B.42), it is useful to replace the chemical potentials with charge densities according to Eq. (10). In the effective Hamiltonian H N , we substitute the leading hard thermal loop contribution to the Hermitian self-energy given by k ·Σ H N = T 2 /8 [127] 15 Besides, assuming an electroweak crossover, there is a contribution from the background Higgs field that starts to continuously follow a temperature dependent expectation value , (B.49) eventually approaching the zero temperature limit v(z → 0) ≈ 174 GeV. An approximate form of v(z) for temperatures above the sphaleron freezeout can be obtained from a lattice calculation which can be found in e.g. [128]. For our purposes it is sufficient to approximate v 2 (z)z 2 by where the temperature dependence is shown in Figure 10. The contributions involving v(z) in Eq. (B.40b) may become important in the limit of degenerate eigenvalues in the Majorana mass matrix M , where it can have a large relative impact on the oscillation frequency, while we neglect the small admixture of doublet leptons to the sterile neutrinos for v(z) = 0. In total, we can decompose the Hermitian part as follows where using Eq. (B.47), the coefficients are given by The result (B.52) immediately leads to Eq. (19) when including spectator effects. Note that if all distribution functions appearing under the momentum integral were of the form f eq (k), the averaging procedure would not incur any inaccuracy. However, this form cannot be assumed for δf 0h (k) and it neither holds for the statistical factor in Eq. (B.42). Nonetheless, since all of these distributions take the form of a Boltzmann tail for |k| a R , the error incurred is only of order one. For comparison, along the same lines, we can see that momentum averaging for leptogenesis from non-relativistic sterile neutrinos does not lead to a leading order inaccuracy because all distributions are well approximated by the Maxwell-Boltzmann form.

C.1 Kinetic Equations
The evolution equations for the charge densities q a of doublet leptons are where S a is the source for the asymmetry (which is defined as the part of the collision term that is non-vanishing even if all µ a = 0) and W a is the rate for washout (which is defined as the remainder of the collision terms). Besides, we introduce the sterile charge as the helicity-odd part of the deviation of the number densities from equilibrium, This quantity is useful because in the limit M → 0, q N can be identified with a charge density contributing to a conserved (modulo weak sphalerons) generalised lepton number along with the doublet leptons and the charged right-handed leptons. In the present context, where the sterile neutrinos are relativistic, we can neglect the reactions that violate the generalised lepton number. This is because for a typical momentum mode the admixture of opposite chirality to a spinor of given helicity is of order M/T , such that the processes mediated by the Yukawa couplings Y that violate the generalised lepton number are suppressed by a relative factor of M 2 /T 2 compared those that conserve the generalised lepton number and that are accounted for in the present work. 16 We also note that there are contributions from the off-diagonal correlations in the sterile neutrinos that we attribute implicitly to the source term S a . In analogy with the terms proportional toΓ a in Eq. (B.52), we refer to the contribution involving q N ii in Eq. (C.1) as a backreaction term. The washout rate is complementary to the damping rates for sterile neutrinos, cf. Eqs. (20), and is given by The off-diagonal correlations of the sterile neutrinos give rise to the source for charge asymmetries in the doublet leptons, where in Eq. (C.1) we use the shorthand notation S a ≡ S aa . While in principle, there can also be off-diagonal correlations in the doublet charges, we set these to zero in the present context because at the temperatures we consider, processes mediated by the µ and τ Yukawa couplings erase these by the mechanism described in Ref. [100,129] corresponding to leptogenesis in the fully flavoured regime [130,131,132]. Because T ∼ |k| M ii for the typical momentum scale, we can focus on the limit of massless (ultrarelativistic) sterile neutrinos, where |k| ≈ sign(k 0 )k 0 . Further, it is useful to decompose sterile propagator in the relativistic regime in terms of helicity odd and even functions (B.43a) and (B.43b), iδS N (k) = 2πδ(k 2 ) −/ kδf even (k) + sign(k 0 )/ kγ 5 δf odd (k) . (C.6) Here we have used Eqs. (B.28) and (B.30a), where the off-shell correlators g 1h and g 2h are suppressed by a factor k 0 /M and where g 3h is related to g 0h . Additionally, the on-shell condition (B.34) has been used in the form of Substitution into the source term (C.4) yields This corresponds to the relativistic limit of the more general result derived in Ref. [115]. Provided we can neglect the term δf 0hij compared to the commutator term in Eq. (B.38) and we can also drop the backreaction term as well as the thermal masses, it follows and we recover the result from Ref. [64] for the source term. This gives rise to a first approximation for the asymmetry in the weak washout regime, which we improve upon in Section 3.
In terms of momentum averaged expressions, the source term can be written as such that in total, we obtain the differential equation for the evolution of the SM charges Eqs. (B.52) and (C.11) form a coupled system of differential equation for the active and sterile charges. In order to solve the whole system one can decompose Eq. (B.52) into even and odd parts and seek for numerical solutions. However, we can identify different parameter regions, such as the oscillatory and overdamped regime, where approximate analytic solutions can be found, as presented in Sections 3 and 4, respectively.

C.2 Spectator Processes
Standard Model processes redistributing charges during leptogenesis are called spectator effects and affect the final result for the baryon asymmetry [85,86]. In order to account for these it is useful to work with the asymmetries ∆ a = B/3 − L a defined in Eq. (12) which are conserved by all interactions other than those mediated by the Yukawa couplings Y between the active and sterile sectors. We then need to relate these asymmetries to the charge densities that appear on the right hand sides of the evolution equations (B.38). At temperatures below T 10 5 GeV, when the electron as the SM particle with the smallest Yukawa coupling finally reaches chemical equilibrium (see e.g. Ref. [87] for an overview of the equilibration rates of spectator processes), the SM Yukawa-mediated processes lead to the constraints Besides, weak and strong sphaleron processes force the relations where Q i denote left-handed quark doublets of flavour i, u i , d i are the corresponding right-handed electroweak singlets and the factor g s = 3 accounts for the three colour states. A common chemical potential for the weak doublets and colour triplets implies that the charge densities associated with the diagonal generators for weak and strong interactions vanish. Correspondingly a vanishing density of weak hypercharge leads to the condition g w Y φ q φ + a=e,µ,τ (g w g s Y Qa q Qa + g w Y a q a + g s Y ua q ua + g s Y da q da + Y ea q ea ) = 0 , (C.14) where we explicitly note the summation over the three active flavour indices. We can now can solve Eqs. (C.12, C.13, C.14) in order to obtain the desired relations between the charge densities of doublet leptons q 1,2,3 ≡ q e,µ,τ as well as of the Higgs bosons q φ and the asymmetries ∆ 1,2,3 ≡ ∆ e,µ,τ . These are conveniently expressed as q = A∆ and q φ = C∆. This way we obtain the matrices A and C given in equation (15)  and where q = (q 1 , q 2 , q 3 ) t as well as ∆ = (∆ 1 , ∆ 2 , ∆ 3 ) t are understood as column vectors in lepton flavour space. For completeness, we also define the column vector q N = (q N 1 , q N 2 , . . . , q Nns ) t for n s sterile neutrinos. Besides, in terms of ∆ we can express the baryon asymmetry as  [134,135]. In view of the sensitivity of the asymmetries from GeV-scale leptogenesis to spectator effects, it should be of interest to include this correction along with the time dependence of the rate of weak sphaleron transitions prior to their quench. Both corrections will lead to a temperature dependence in above conversion relations, a detailed study of which we leave to future work.

D.1 Time Scales in the Oscillatory Regime
For the validity of the approximations used to calculate the initial asymmetry in the oscillatory regime, the equilibration time given here by the inverse of the smallest eigenvalue of the decay matrix (B.41) needs to be much larger than the time by which the first oscillation is over. This oscillation time scale is determined by the difference of the squared masses In the coordinates we have chosen, Eq. (B.52) implies that the frequency of the oscillation ω vac induced by the vacuum term H vac N increases with z 2 , whereas the thermal contribution H th N results in a constant oscillation frequency ω th . For this reason the nonzero thermal oscillation may be of importance at early times when the vacuum oscillation has not started yet. However, one can show that ω vac is automatically larger than ω th at the time of the first oscillation z osc when imposing z eq z osc : with h th = 0.23. This implies that in the oscillatory regime the thermal effects may only have lead to a small fraction of a full flavour oscillation by the time when the first oscillation due to the vacuum masses already has been completed. Since the main part of the active charge is generated during the first oscillation, one can consider the contribution from the thermal masses as a perturbation. Furthermore, in the oscillatory scenario the effect of the temperature dependent Higgs vacuum expectation value can be neglected completely since it just becomes import at times z z osc . It is easy to show that the perturbative corrections to δn ij arising due to the presence of H th N vanish at order O(h th /γ av |Y * Y t |) as the leading order term of the out-of-equilibrium distribution is δn ij = −n eq δ ij and hence The first non-vanishing contribution from the thermal masses is of order O(h th /γ av |Y * Y t | 2 ), which can be neglected compared to the contributions coming from the vacuum masses δn ij of order O(|Y * Y t |), cf. Eqs. (31).

D.2 Momentum Dependence of the Source
In Section 3 we have calculated the active charge produced through the off-diagonal oscillations of the sterile neutrinos to order |Y * Y t | 2 with the simplification of fully momentum averaged expressions. We can go one step further and consider the momentum dependence of the vacuum term H vac N as in Eq. (B.40a) but still keep the replacement (B.48) in order to able to solve the remaining momentum integral analytically. For this reason we solve by analogy with Eq. (26) for the even and odd parts of the off-diagonal distributions δf ij whose solution to order |Y * Y t | can be obtained analogously while the integration over z remains unchanged, so that the active charge is given bỹ with a function that carries all momentum information Solving this integral exactly is beyond the scope of this paper since k ·Σ A N has a non-trivial momentum structure [88]. Nevertheless, we can use the momentum averaged replacement (B.48), which leaves us with a momentum integral that can easily be solved analytically: whereas we expect the error in Eq. (D.9) of to be of order one [88] and hence sufficient our purposes.

D.3 Sterile Charges in the Oscillatory Regime
In Section 3 we have pointed out that up to order |Y * Y t | 2 no sterile charge q N is generated by the off-diagonal oscillations. We will show in the following that this is true to all orders for n s = 2 sterile flavours, whereas this is not true for n s ≥ 3 since a non-vanishing contribution appears at O(|Y * Y t | 3 ). In order to do so, we introduce a function for i = j. Its derivative with respect to z reads (D.14) The deviations δn even ji and δn odd ji are determined by solving Eq. (26) for non-diagonal components (i = j). In case of n s = 2 flavours, one can express the anticommutators as since the diagonal entries of Y * Y t are purely real due to its Hermitian property. After some calculation we are left with which is true for all z. Additionally, this even results in a condition between n even and n odd : Whereas this holds for n s = 2 sterile flavours one can show that for n s ≥ 3, already at O(|Y * Y t | 3 ), there appears a non-vanishing contribution to F i . For that, we solve Eq. (26) for off-diagonal δn ij recursively to O(|Y * Y t | 2 ) by using solutions for δn at O(|Y * Y t |). This result can be used as an input for F i in Eq. (34), such that for n s = 3, we have: with | ijk | as the absolute value of the Levi-Civita-Symbol in order to account for (i = j, k = i, k = j). Thus, as a perturbative expansion in the Yukawa coupling Y , it is justified to assume zero initial sterile charge q N ii after the first oscillations in the oscillatory regime.