Leptogenesis in GeV scale seesaw models

We revisit the production of leptonic asymmetries in minimal extensions of the Standard Model that can explain neutrino masses, involving extra singlets with Majorana masses in the GeV scale. We study the quantum kinetic equations both analytically, via a perturbative expansion up to third order in the mixing angles, and numerically. The analytical solution allows us to identify the relevant CP invariants, and simplifies the exploration of the parameter space. We find that sizeable lepton asymmetries are compatible with non-degenerate neutrino masses and measurable active-sterile mixings.


Introduction
One of the interesting potential implications of (Majorana) neutrino masses is the generation of a matter-antimatter asymmetry in the Universe. It has been demonstrated that the generation of sizeable leptonic asymmetries, leptogenesis, is generic in extensions of the Standard Model that can account for neutrino masses [1]. In particular two new ingredients are essential for this mechanism to work: the existence of new weakly interacting particles that are not in thermal equilibrium sometime before the electroweak phase transition and the existence of new sources of CP violation.
Leptogenesis from the out-of-equilibrium decay of heavy Majorana fermions that appear in type I seesaw models [1] has been extensively studied (for a comprehensive review see e.g. [2]). The simplest version requires however relatively large Majorana masses > 10 8 GeV [3,4] (or > 10 6 if flavour effects are included [5]), which imply that this scenario would be very difficult to test experimentally. It is possible to have sizeable asymmetries for smaller masses if a large degeneracy exists, through resonant leptogenesis [6].
On the other hand, for Majorana masses in the GeV range, when the neutrino Yukawa couplings are small, another mechanism might be at work. In particular, the non-equilibrium condition takes place not in the decay, but in the production of the heavy sterile neutrinos. The small Yukawa couplings imply that some of the species might never reach thermal equilibrium and a lepton asymmetry can be generated at production and seed the baryon asymmetry in the Universe. This mechanism was first proposed by Akhmedov, Rubakov and Smirnov (ARS) in their pioneering work [7] and pursued, with important refinements in refs. [8,9]. For a recent review and further references see [10]. In most of these works, the case of just two extra sterile species is considered, which is also the limiting case of the so-called νMSM where there are three species, but one of them plays the role of warm dark matter (WDM) and is almost decoupled, having no impact in the generation of the lepton asymmetry. When the mechanism involves just two species, it has been found that the observed baryon asymmetry is only possible if the two states are highly degenerate in mass. This however was not the conclusion of the ARS paper.
The purpose of this paper is to explore systematically the parameter space in the case of three sterile species (which encompass the one with two neutrinos) as regards the baryon asymmetry, in particular we do not want to restrict the parameter space to have a WDM candidate. The model has many free parameters (only 5 out of the 18 parameters are fixed by the measured light neutrino masses and mixings) and the exploration of the full parameter space is challenging. Only with the help of approximate analytical solutions to the kinetic equations this task is feasible. The analytical solutions furthermore allows us to identify the relevant CP invariants and to reach regions of parameter space where the equations become stiff and very difficult to deal with numerically.
The paper is organised as follows. In section 2 we present the model, which is essentially a generic type I seesaw model, establish the notation and discuss on general grounds what are the CP reparametrization and flavour invariants we expect to find in computing any CP violating quantity such as any putative lepton asymmetry. In section 3 we present the kinetic equations that describe the production of sterile neutrinos and solve them analytically via a perturbative expansion in the mixing angles up to the third order. In section 4 we compare the analytical and numerical solutions for several choices of the parameters, and identify the region of parameter space where the analytical solution accurately describes the numerical one. In section 5 we use the analytical solutions and perform a Monte Carlo scan (using the software package MultiNest [11,12]) to find regions of parameter space that can reproduce the observed baryon asymmetry, and that are compatible with the measured neutrino masses and mixings. In section 6 we conclude.

Minimal Model of neutrino masses
We will concentrate on the arguably simplest model of neutrino masses that includes three right-handed singlets. The Lagrangian is given by: where Y is a 3 × 3 complex matrix and M a diagonal real matrix. The spectrum of this theory has six massive Majorana neutrinos, and the mixing is described in terms of six angles and six CP phases generically. One convenient parametrization for the problem at hand is in terms of the eigenvalues of the yukawa and majorana mass matrices together with two unitary matrices, V and W . In the basis where the Majorana mass is diagonal, M = Diag(M 1 , M 2 , M 3 ), the neutrino Yukawa matrix is given by: Without loss of generality, using rephasing invariance, we can reduce the unitary matrices to the form 1 : W = U (θ 12 , θ 13 , θ 23 , δ) † Diag(1, e iα 1 , e iα 2 ), V = Diag(1, e iφ 1 , e iφ 2 )U (θ 12 ,θ 13 ,θ 23 ,δ), Obviously not all the parameters are free, since this model must reproduce the light neutrino masses, which approximately implies the seesaw relation: where v = 246 GeV is the vev of the Higgs. On the other hand, the known neutrino masses and mixings do not give us enough information to determine the Majorana spectrum, not even the absolute scale. Very strong constraints can be derived from neutrino oscillation experiments for masses below the eV range [13,14,15,16]. Cosmology can exclude a huge window below 100 MeV [17,18,19,20,21,22], except maybe for one species that could be lighter provided the lightest active neutrino mass is below 3 × 10 −3 eV [20,21]. The GeV range is interesting because an alternative mechanism for lepton asymmetry generation could be at work [7,8,9]. Majorana neutrinos in this range are heavy enough to safely decay before Big Bang Nucleosynthesis, while they are light enough that they might have not reached thermal equilibrium by the time of the electroweak phase transition (EWPT), behaving as reservoirs of a putative lepton asymmetry. Our goal in this paper is to explore the full parameter space of this model allowed by neutrino masses, as regards leptogenesis. An essential condition will be that at least one 1 Although we use the same notation for the mixing angles and phases of W as those in the usual PMNS matrix, they should not be confused. 2 Note the unconventional ordering of the 2×2 rotation matrices in U .
of the sterile neutrinos does not reach thermal equilibrium before the EWPT. This can be ensured assuming a large hierarchy in the yukawas [7]: It is mandatory, however, to have an accurate analytical description, since the unconstrained parameter space is huge. We will solve the quantum kinetic equations in a perturbative expansion in the mixings in the next section. Since the lepton asymmetry is necessarily a CP-odd observable, on general grounds we can derive what are the expectations in terms of weak-basis CP invariants.

CP invariants
In [23], weak basis (WB) invariants sensitive to the CP violating phases which appear in leptogenesis, within the type I seesaw model, were derived. All of them should vanish if CP is conserved, and conversely the non-vanishing of any of these invariants signals CP violation. They must be invariant under the basis transformations: Defining h ≡ Y † Y , and H M ≡ M † M , a subset of the invariants can be written as: Since the I i are WB invariants, we can evaluate them in any basis. In the WB where the sterile neutrino mass matrix M is real and diagonal, one obtains: 10) j and, using the parametrization of eq. (2.1) It is explicit in the above expression that such unflavoured invariants depend only on the CP phases of the sterile neutrino sector, which are encoded in the unitary matrix W : one Dirac-type phase, δ and two Majorana-type phases α 1 , α 2 . Not surprisingly, these invariants are the relevant ones in unflavoured leptogenesis, i.e., in the conventional computation of the CP asymmetry generated by heavy Majorana neutrino decay neglecting flavour effects.
The combinations of W matrix elements which appear in Im(h 2 ij ) can be expressed in terms of the rephasing invariants defined in [24] as follows: (2.14) Notice that J W ≡ ±Im[W αi W * βi W * αj W βj ] is the Jarlskog invariant for the matrix W , while the quantities Im[(W αj W * αi ) 2 ] determine the Majorana phases, α 1,2 . When considering processes, such as heavy neutrino oscillations, where the Majorana nature does not play a role, only the Dirac phase δ will be relevant and therefore we expect to find just the Jarlskog invariant of the matrix W.
Since there are six independent CP-violating phases, it is possible to construct three more independent WB invariants, which would complete the description of CP violation in the leptonic sector. One simple choice are those invariants obtained from I i under the change of the matrix h byh ≡ Y † h Y , with h = λ λ † , being λ the charged lepton Yukawa couplings, i.e.,Ī and analogously forĪ 2 ,Ī 3 . The corresponding CP odd invariants are Im(h 2 ij ), which in the basis where also the charged lepton Yukawa matrix is real and diagonal can be written as: The lepton number (L) violating part of the flavoured CP asymmetries in leptogenesis depends on the above combinations [25]: wheref is an arbitrary function. Upon substitution of the neutrino Yukawa couplings as given in eq. (2.1) can be written as: These asymmetries contain the additional rephasing invariants of the form Im[W * βi V βα V * δα W δj ], which depend on the phases in the matrix V (δ, φ 1 , φ 2 ), showing that the flavoured CP asymmetries of leptogenesis are also sensitive to the CP phases in the V leptonic mixing matrix, besides those in W .
Alternatively, we choose to construct the WB invariants which will appear when the Majorana character of the sterile neutrinos is not relevant, i.e., L-conserving ones. These are given by: 20) where The L-conserving CP asymmetry in leptogenesis via heavy neutrino decay, as well as the CP asymmetries encountered in leptogenesis through sterile neutrino oscillations, are sensitive to the above combinations of Yukawa couplings [25]: where f is an arbitrary function, and can be written in terms of the rephasing invariants as: Notice that the crucial difference between the L-violating and the L-conserving CP asymmetries is that in L iα the combination of W matrix elements is such that all dependence on the Majorana phases α 1,2 disappears, as expected.
In the approximation of neglecting y 3 y 1 , y 2 , we obtain that so in principle we expect that the lepton asymmetry will depend on ten CP invariants, namely Im[W * 1i V 1α V * 2α W 2i ], with i = 1, 2, 3 and α = 1, 2, 3 and J W . However, they are not all independent. In ref. [24] it has been shown that in the minimal seesaw there are only six independent CP invariants that can be made out of the matrices V, W . Two of them correspond to the Majorana phases of W , α 1,2 , which as we have argued before will not contribute in the limit of small sterile neutrino Majorana masses that we are considering. Other two are the equivalent of the Jarlskog invariants for the matrices V, W and therefore determine the Dirac phases,δ, δ, respectively. The last two are of the form Im[W * 1i V 1α V * 2α W 2i ], for two reference values of i, α, that fix the additional phases φ 1,2 . Moreover, it can be shown that since we are neglecting the Yukawa coupling y 3 , the phase φ 2 of the matrix V does not appear in eq. (2.25), thus we are left with only three independent invariants.
The unitarity of the mixing matrices V, W implies that which allows to write the invariants Im[W * 1i V 1α V * 2α W 2i ] for α = 2 in terms of those with α = 1, 3, and the invariants for i = 2 in terms of the corresponding ones with i = 1, 3. By exploiting the identities (2.28) we can write for instance one of the invariants with β = 3 in terms of the invariant with α = 1 and the Jarlskog invariant for V , It is simpler, though, to write the results in terms of the following four invariants, even if only three are independent, expanded up to 3rd order in the small mixing angles θ ij ,θ ij : A generic expectation for the CP-asymmetry relevant for leptogenesis is Since the CP rephasing invariants are at least second order in the angles, we just need to take the diagonal elements in ∆ CP , to keep the result up to 3rd order. Then, in the limit y 3 = 0, we get: We will see in the next section that this is precisely the yukawa and mixing angle dependence we will find when solving the kinetic equations, which is a strong crosscheck of the result.

Sterile neutrino production
Our starting point is the Raffelt-Sigl formulation [26] of the kinetic equations that describe the production of sterile neutrinos in the early Universe. The density matrix is the expectation value of the one-particle number operator for momentum k: ρ N (k) for neutrinos, andρ N (k) for antineutrinos. We will assume that only sterile neutrinos and the lepton doublets are out of chemical equilibrium, but assume that all the particles are in kinetic equilibrium, using Maxwell-Boltzmann statistics: where ρ eq (k) ≡ e −k 0 /T , with k 0 = |k|, and µ a denotes the chemical potential normalised by the temperature. We will furthermore neglect spectator processes and the washout induced by the asymmetries in all the fields other than the sterile neutrinos and lepton doublets. We expect this approximation to give uncertainties of O(1) which for our purpose is good enough [27]. In [7], only the asymmetry in the sterile sector was considered, neglecting the feedback of the leptonic chemical potentials. In this case, the equations get the standard forṁ and the analogous forρ N with H → H * , where H is the Hamiltonian (we neglect matter potentials for the time being but we will include them later on) Γ is the rate of production/annihilation of sterile neutrinos in the plasma, which is diagonal in the basis that diagonalises the neutrino Yukawa's: where we assume y 3 = 0. In deriving eq.(3.2) it is assumed that the particles involved in the production/annihilation of the sterile neutrinos are in full equilibrium (all chemical potentials vanish), and that kinematical effects of neutrino masses are negligible. Note that only the matrix W appears in these equations and therefore any CP asymmetry generated can only be proportional to the invariant J W which depends at third order on the mixing angles of W .
In [8] it was correctly pointed out that the asymmetries in the sterile sector will be modified by the leptonic chemical potentials that will be generated as soon as sterile neutrinos start to be produced. Including the evolution of the leptonic chemical potentials has two important consequences: new sources of CP violation become relevant and washout effects are effective. Leptons are fastly interacting through electroweak interactions in the plasma and therefore it is a good approximation to assume they are in kinetic equilibrium.
An important question is what is the flavour structure of these chemical potentials. For T 10 9 GeV the Yukawa interactions of the tau and muon are very fast, which implies that µ will be diagonal in the basis that diagonalises the charged lepton Yukawa matrix, since no other interaction changing flavour is in equilibrium before the heavy neutrinos are produced. Note however that this is not the basis where the neutrino Yukawas are diagonal, the two are related by the mixing matrix V . As a result, when the evolution of the lepton chemical potentials is taken into account, the CP phases of the matrix V become relevant.
Adapting the derivation of [26] to this situation, we find that the evolution of the CP-even and CP-odd parts of the neutrino densities: ρ ± ≡ ρ N ±ρ N 2 and the lepton chemical potentials , µ α , to linear order in µ α , ρ − , satisfy in this case: , I α is the projector on flavour α and γ a,b N , γ a,b ν are the rates of production/annihilation of a sterile neutrino or a lepton doublet neglecting all masses, after factorizing the flavour structure in the Yukawas, .

(3.7)
A similar derivation can be found in [28] and we agree with their findings. These equations reduce to those in eq. (3.2) in the limit µ → 0 with Γ i = y 2 i (γ a N + γ b N ). Most previous studies have assumed that the rates are dominated by the top quark scatterings. In this case, the rates are given (in the Boltzman approximation) by the well-known result [29,30] The factor of 2 difference between the rates of the N and the ν is due to the fact that the lepton is a doublet and the sterile neutrino is a singlet. Note that there is a non-linear term of the form O(µρ + ), as first noted in [28]. More recently in [31], the equations have been written in terms of the µ B−Lα/3 chemical potentials, however not all the chemical potentials (e.g. higgs and top quark) have been included. A full treatment including all chemical potentials will be postponed for a future work, but we expect that including these spectator effects will change the results by factors of O(1).
In [30,32], it has been pointed out that the scattering processesLN ↔ W H get a strong enhancement from hard thermal loops and are actually the dominant scatterings. The results of [30,32] however do not include the chemical potentials of spectators, so it is not clear how to include them consistently in the above equations. We will neglect these effects in the following. Note however that the lepton flavour structure of these and of the top quark scatterings is the same.
It is easy to see also that total lepton number is conserved as it should: Two approximations are often used in solving these equations: 1) assume that the momentum dependence of ρ ± follows that of ρ eq , i.e. kinetic equilibrium for the sterile states, which implies r ± = ρ ± /ρ eq are constants and the integro-differential equations become just differential equations, 2) neglect the k 0 dependence of the rates by approximating The effect of these approximations has been studied numerically in [28] and the results do not differ too much. We will therefore adopt both approximations that simplify considerably the perturbative treatment.

Lepton asymmetries in the sterile sector
We are going to solve these equations perturbing in the mixing angles up to third order. We first consider the simpler case, neglecting leptonic chemical potentials and considering in turn the evolution in a static Universe and in the expanding case.

Static Universe
We start with eq. (3.2) and assume y 3 = 0. In this case, neither H nor Γ depend on time.
Defining ρ N ij /ρ eq ≡ a ij + ib ij and taking into account the hermiticity of ρ N we change the matrix equation into a vector equation: r ≡ (a 11 , a 22 , a 12 , b 12 , a 13 , b 13 , a 23 , b 23 , a 33 ). (3.11) At 0-th order the system of equations of eq. (3.2) can be rewritten aṡ 13) and the matrix A 0 is constant and has a block structure: The matrix can be easily diagonalised and exponentiated so the general solution to the equation is At the next order we have to keep O(θ ij ) in the Hamiltonian and translate the matrix form into the vector form: The equation for the first order correction to the density iṡ The solution at this order is therefore We can iterate this procedure to get the correction at order n: with solution We can define the evolution operator so that the solution can be written as As a first estimate of the leptonic asymmetry that can be generated, we are interested in ∆ρ 33 since this is the sector that will never reach equilibrium (in the absence of mixing) and therefore can act as reservoir of the leptonic asymmetry until the electroweak phase transition [7]. One can easily compute the solution of the eq. (3.21) up to order n = 3, which is the first order that gives a non-vanishing result, as expected from general considerations on CP invariants. The result at finite t is not particularly illuminating but the limit t → ∞ is rather simple: . A few comments are in order. We have not assumed any expansion in Γ i in this expression, only in the mixing angles. According to general theorems the equations should reach a stationary solution if all the eigenvalues of the matrix A 0 + A 1 + A 2 + ... are real and negative. However, because Γ 3 = 0, one of the eigenvalues of A 0 vanishes and it is lifted only at second order in perturbation theory, ∼ θ 2 i3 Γ i , therefore we expect the perturbative expansion should break down for t ∼ 1 , which is the time scale of equilibration of the third state. On the other hand, if θ is small, the perturbative solution should be accurate for times t ≥ Γ −1 1 (2) . Indeed this is precisely what we find comparing the perturbative and numerical solutions in figure 2.
The result is proportional to J W which is the only CP rephasing invariant that can appear in this case. The result vanishes if any two of the masses or the yukawa's are degenerate, since the CP phase would be unphysical in this case.

Expanding Universe
Let us turn now to the realistic case of an expanding Universe. As usual, we will consider the evolution as a function of the scale factor x ≡ a, in such a way that the Raffelt-Sigl equation becomes where H u (x) is the Hubble parameter, H u = 4π 3 g * (T )

45
T 2 M Planck , and y ≡ p T . Assuming for simplicity a radiation dominated Universe with constant number of degrees of freedom, during the sterile evolution time we can assume xT = constant that we can fix to be one. Therefore the scaling of the different terms is where M * P ≡ M Planck 45 4π 3 g * (T 0 ) and g * (T 0 ) is the number of relativistic degrees of freedom in the plasma during the sterile evolution.
Therefore the equation as function of x is: where we have defined The perturbative expansion works as in section 3.2.1, but now all the A n (x) are x-dependent: A n (x) with n ≥ 1 scale like the Hamiltonian, ie. x 2 , while A 0 (x) contains terms that scale with x 2 and others that do not depend on x. Fortunately, there is an important simplification in that A 0 (x) can be diagonalised by an x-independent matrix, therefore the path-ordered exponential can be easily evaluated. The result can be written in the same form of eq. (3.23), with the evolution operator given by At third order in the mixings, after algebraic simplifications and partial integrations, the result can be given in terms of integrals of the form where α i are combinations of ∆ ij and β i are combinations of γ i . Up to third order in the perturbative expansion only integrals with n ≤ 3 appear.
The integrals in the range [t 0 , t] can be approximated by the large t behaviour of the J 1n (α, t) functions, after resumming the Taylor series in β i . Further details are presented in appendix A. The finite t dependence of the asymmetry ∆ρ 33 is rather complicated, but the asymptotic value is non-zero and rather simple:

(3.35)
This can be simplified to where Comparing eq. (3.36) and eq. (3.24) we see that in the expanding case the asymmetry is cubic in γ i and not linear. Note that the dependence on the yukawa's is precisely that expected from a flavour invariant CP asymmetry. In fact this is effectively the situation in the expanding case, because the asymmetry is generated at times t γ −1 i and the dependence in the yukawa's in this regime is therefore perturbative. This is in contrast with the non-expanding case, where the asymmetry evolves all the way till t ∼ γ −1 i . To understand the reason behind this different behavior, it is useful to recall the definition of ∆ ij from eq. (3.28). Then, we see that ∆ ij x 3 1 implies ∆M 2 ij /(4T ) T 2 /M * P = H u (T ), therefore in this regime the sterile neutrino oscillations are much faster than the Hubble parameter and no asymmetry is produced anymore, since oscillations are averaged out. Thus in the expanding Universe the generation of the asymmetry occurs at x ∼ |∆ ij | −1/3 γ −1 i . Until now we have neglected the matter potentials, however given the suppression in three powers of γ of the leading result, there are corrections of same order coming from the potentials, and in fact they are numerically more important.
The equation including the potentials in the basis with diagonal neutrino Yukawas is: where The result for the asymmetry including the potentials is given by: The leading terms O(v 2 γ) at asymptotic times t γ −1 1,2 are: depends only on the ratios of mass differences and/or the ordering of the states. This result is parametrically the same as the result of [7] if we neglect the dependence of κ on the mass differences and has the dependence on the yukawas expected from eq. (2.34). Considering the naive seesaw scaling y 2 i ∼ 2 mν M i v 2 , for m ν ∼ 1 eV and assuming no big hierarchies or degeneracies, i.e.
The asymmetry is highly sensitive to the light neutrino mass. Note that we have pushed the value to the limit, a light neutrino mass in the less constrained 0.1 eV range would imply three orders of magnitude suppression. The asymmetry grows linearly with the mass of the heavy steriles. However, for masses larger than ∼ 10 -100 GeV lepton number violating transitions via the Majorana mass could washout further the asymmetry, an effect that requires a refinement of the formulation to be taken into account.

Lepton asymmetries in the active sector
The asymmetry generated ignoring the µ evolution depends only on the Dirac-type phase, δ, appearing in W as we have seen. However when the evolution of the leptonic chemical potentials is included, other phases contribute to the total lepton asymmetry. We will perform a perturbative expansion to third order in the mixings of both V and W matrices. The result at finite t θ 2 i3 (θ 2 i3 )γ −1 i can be written in the form: where all the four CP invariants appear, I CP = J W , I 1 , I 1 , I , given in eqs. (2.29). At finite t, the result for the functions A I CP is well approximated by and where we have definedγ i ≡ y 2 iγ N and ∆ v ≡ v 2 − v 1 , and the result for G 3 (t), G 41 (t), G 42 (t) are lengthier and reported in the appendix B. These results would get modified for γ i t 1 had we included the non-linear terms that modify the rate of thermalisation at large times. In these equations there is an implicit expansion up to third order in γ i (v i )/∆ 1/3 when ∆ 1/3 t 1, while the terms γ i (v i )t are resumed. In figure 3 we plot the functions A I (2) 1 (t) and A I (3) 1 (t), which depend only on one neutrino mass difference. We show two physical situations: one with very degenerate neutrinos and the other with no strong degeneracies.
These two invariants are the only ones relevant for the scenario that has been considered in most previous studies, where it has been assumed that only two sterile neutrinos have a role in generating the lepton asymmetry (see for instance [33] for a very recent analysis). This is the situation in the limit of complete decoupling of N 3 , ensured by the condition θ i3 = 0, implying that only the invariants I survive. In [8] an approximate analytical solution was obtained, expanding in the yukawa's, under the assumption that In this limit, the result of eqs. (3.46) and (3.47) can be simplified to we find: while for y 2 1 = y 2 2 /2 = 10 −14 (that would correspond to light neutrino masses in the eV range and heavy ones in the GeV range) we would have

(3.52)
Even if the CP invariants are of O(1), the asymmetry is too small unless there is a significant degeneracy between the two states [8]. It is important however to realise that the naive seesaw scaling is too naive and a full exploration of parameter space is necessary. (t) and A J W (t). They depend on the two neutrino mass differences, so we show three examples here: one in which there are no degeneracies, one where there are two almost degenerate states, and the case where the three states are almost degenerate. As in the previous case we see a large enhancement when only one of the mass differences is small and a further enhancement when the two are small compared to the absolute scale. In the case of A J W we find that there is a significant difference in the regime ∆ 1/3 ij t 1 if we plot A J W (t) truncated to the terms of O(y 6 i ). As we will see in the next section, the latter is much closer to the numerical result. The reason for this difference is that at small times, ∆ 1/3 t 1, only some terms of order O(y 8 i ) are kept in eqs. (3.46), while there is a strong cancellation if all had been included. Note however that this effect is only important at times where the asymmetry is suppressed and seems to affect only A J W .
It is interesting to note that even though the dependence on the yukawas of the functions A I CP (t) is different (fourth or sixth order), the maxima for all cases are roughly of the same order of magnitude. Note, however, that in the limit t γ −1 i , only the contribution of two invariants, J W and I where we kept only the leading terms O(y 4 ) proportional to I and we have used the result of eq. (3.50).
The first term in this expression corresponds to the expectation of [7], ie. the final asymmetry is proportional to that stored in the third sterile state, eq. (3.40), while the second term was missing in the simplified treatment of [7]. Note that they depend on different CP invariants.
1 , compared with the prediction, A I (2) 1 (t) (dashed red). Right: same for case 2 normalised to the invariant I The parameters are the same as in figure 3 for the degenerate case.

Numerical solution
In order to check the accuracy of the analytical solutions presented in the previous section, we have solved the differential equations numerically. As shown in [28], the momentum dependence does not change significantly the results so we will consider the averagemomentum approximation.
In figures 5-6 we compare the analytical and numerical solutions for the functions A I CP (t) in the highly degenerate case (the values of the mixing angles are of O(10 −2 )) . In order to isolate the appropriate invariant we make the following choices: • Case 1: θ i3 =θ i3 = 0 isolates I • Case 4:θ ij = 0 isolates I J W .
The numerical results normalised by the corresponding CP invariant are shown together with the predictions of the previous section. In the case of J W , we plotted the function A J W keeping only the terms of O(y 6 ) that is more accurate at small t and the full function at large t. The agreement in all cases is quite good. The differences observed at large t come from the non-linear terms in the equations. We also show the numerical results of the equations without them and find a very good agreement also at large t. Note that the approximation works well in the regime γt 1, that is in the strong washout regime of the fast modes.
Numerically it is very hard to go to regimes where the ratios γ/|∆| 1/3 become very small, since the system becomes stiff. On the other hand, there is no reason why the perturbative solution is not accurate in such regime. We will therefore assume this to be The parameters are the same as in figure 4 for the double degenerate case. the case in the following section and use the perturbative solution to perform a scan of parameter space.

Baryon asymmetry
The observed baryon asymmetry is usually quoted in terms of the abundance, which is the number-density asymmetry of baryons normalised by the entropy. After Planck this quantity is known to per cent precision [34]: The lepton asymmetries in the left-handed (LH) leptons generated in the production of the sterile neutrinos are efficiently transferred via sphaleron processes [35] to the baryons. The baryon asymmetry is given by Since we have neglected spectator processes in the transport equations, the B − L asymmetry is related to the chemical potentials computed in the previous sections by the relation where g * = 106.75 (which ignores the contribution to the entropy of the sterile states). Our estimate for the baryon asymmetry is therefore We have performed a first scan of the full parameter space of the model. Given the theoretical uncertainties mentioned in different sections of the paper, we have considered as interesting the points that can explain the baryon asymmetry within a factor of 5. For this we have used the analytical solutions, even though in some regions of parameter space they will not be precise, since they are based on a perturbative expansion on the mixing angles of the matrices V and W . We have considered however a few cases where the angles are not small and we find that the analytical solutions differ from the numerical ones only in some global numerical factor of a few , but the time dependence is very similar.
Even with an analytical expression the exploration of the large parameter space is a challenge. We have used the package Multinest [11,12] to perform a scan on the Casas-Ibarra parameters [36], where the Yukawa matrix is written as m light is a diagonal matrix of the light neutrino masses and R is a complex orthogonal matrix that depends on three complex angles z ij . We fix the light neutrino masses and mixings to the present best fit points in the global analysis of neutrino oscillation data of ref. [37] and leave as free parameters: three complex angles, the three phases of the PMNS matrix, the lightest neutrino mass as well as the heavy Majorana masses that are allowed to vary in the range M i ∈ [0.1, 100] GeV. In total thirteen free parameters. The scan searches for minima of the quantity | log 10 |Y B (t EW )/Y exp B || (in the range ≤ 1.5) and the MultiNest algorithm is optimised to sample properly when there are several maxima. For the determination of Y B we use the analytical results of the previous sections, for which the CP invariants are computed directly from the matrix elements of the V, W matrices that can be easily calculated by diagonalising the Yukawa mass matrix obtained in the Casas-Ibarra parametrization. Since the mechanism to work requires that at least one of the modes does not get to equilibrium before the electroweak phase transition we restrict the search to the range where one of the yukawa eigenvalues, y 3 , is much smaller than the others and the following conditions are satisfied Furthermore, since the kinetic equations neglect lepton number violating effects in the rates, we impose additionally the constraint We first consider a case where one of the sterile neutrinos is effectively decoupled from baryon number generation, that we can assume to be N 3 . This can be achieved with the choice of parameters: contribute. This case is the one that has been considered in most previous works on the subject [8,9,38,28,31], where the number of parameters is reduced to six: only one complex angle, two PMNS CP phases and two Majorana neutrino masses are relevant.
It is believed that a large degeneracy of the two sterile neutrinos is needed to obtain the correct baryon asymmetry. In figure 7  (blue,green,red). Successful leptogenesis is possible in a larger range of parameter space for IH than for NH. In the range shown our results agree reasonably well with those in ref. [39] for the IH, while the range for NH looks a bit smaller. We see that there are a significant number of points for which the degeneracy is mild for the IH. We have analysed more carefully some of these points by solving the full numerical equations. We find that even though these points correspond to cases where the angles in V, W are not small, the analytical and numerical solution agree very well and have the same t dependence as shown in figure 8. Note that the numerical solution is difficult at large times for non-degenerate solutions and the standard methods that we use fail. An optimised numerical method is needed to solve the stiffness problem and this will be studied elsewhere. It is very interesting to correlate the baryon asymmetry with observables that could be in principle measured such as the Dirac CP phase of the PMNS matrix, the amplitude of neutrinoless double beta decay or the active-sterile mixings that control the probability for the heavy sterile states to be observed in accelerators or in rare decays of heavy mesons. The effective mass entering the 0νββ decay is given by where M 0νββ are the Nuclear Matrix Elements (NMEs) defined in [40] 3 . The first term corresponds to the standard light neutrino contribution and the second is the contribution from the heavy states. U ei with i ≥ 4 is the active-sterile neutrino mixing.
In figure 9 we show the results for the active-sterile mixing as function of the sterile mass and compare them with present direct bounds and the prospects of SHiP [43] and LBNE near detector [44]. We show the result for M 1 but the one for M 2 is almost identical. We see that most of the parameter space for successfull baryogenesis is not excluded by present constraints and that the active-sterile mixings tend to be larger for the IH. A sizeable region in the range of the GeV could be explored in the future experiment SHiP in the case of the IH and by LBNE near detectors. It is interesting to note that the less degenerate solutions can not have very small active-sterile mixing, as shown in figure 10, where we plot the points on the plane deg ≡ |M 2 − M 1 |/(M 2 + M 1 ) versus the active-sterile mixing in the electron flavour. The degeneracy can be lifted to some extent at the expense of larger yukawa couplings which also imply larger mixings.
We have looked for direct correlations of the baryon asymmetry with the phases of the PMNS matrix. We have found that the distribution on the Dirac phase and the Majorana phase are flat. This is due to the fact that the complex angle can provide the necessary CP violation, even if the PMNS phases would vanish. The same is true for the effective mass of neutrinoless double beta decay, which depends on the Majorana phase. A dedicated scan is needed to quantify how the putative measurement of various observables could constrain the lepton asymmetry. This will be done elsewhere.
In the general case, N 3 is also relevant and the main difference with respect to the previous situation is that there is a significantly enlarged parameter space where degeneracy is not necessary. This was already found in refs. [46] for some points of parameter space. (green) for NH (up) and IH (down), with only two sterile neutrino species. The red bands are the present constraints [45], the solid black line shows the reach of the SHiP experiment [43] and the solid red line is the reach of LBNE near detector [44]. In figure 12 we show the points on the plane (∆M 12 , M 1 ) for the general case. The activesterile mixings are shown in figure 13. These mixings can be larger in this case, specially in the case of the NH. The SHiP prospects are therefore more promising in this context. As in the N = 2 case there is no direct connection between the asymmetry and the PMNS CP phases. On the other hand, the lightest neutrino mass is non-zero in this case, but the requirement that one yukawa needs to be significantly smaller than the others, eq. (5.6),  implies that the lightest neutrino mass must be small. In figure 11 we show the distribution of this quantity for those points that satisfy Y B ≥ Y exp B in the case of NH (the IH being very similar).

Conclusions
We have studied the mechanism of leptogenesis in a low-scale seesaw model that is arguably the simplest extension of the Standard Model that can account for neutrino masses. For Majorana neutrino masses in the GeV range, sizeable lepton asymmetries can be generated in the production of these states some of which never reach thermal equilibrium before the electroweak phase transition. Lepton asymmetries are efficiently transferred to baryons via sphaleron processes. This mechanism was proposed in [7,8] and studied in many (green) for NH (up) and IH (down), with three sterile species. The red bands are the present constraints, the solid black line shows the reach of the SHiP experiment [43] and the solid red line is the reach of LBNE near detector [44].
works, but a full exploration of parameter space in the general case of three neutrinos is lacking. To this aim we have developed an accurate analytical approximation to the quantum kinetic equations which works both in the weak and strong washout regimes of the fast modes (there is always a slow mode that does not reach thermal equilibrium before the EW phase transition). It relies on a perturbative expansion in the mixing angles of the two unitary matrices that diagonalise the Yukawa matrix. This analytical approximation allows us to identify the relevant CP invariants, and explore with confidence the regime of non-degenerate neutrino masses which is very challenging from the numerical point of view. We have used this analytical solution to scan the full parameter space using the MultiNest package to identify the regions where the baryon asymmetry is within an order of magnitude of the experimental value. We have performed first a scan in the simpler setting where one of the sterile neutrino decouples, which reduces the parameter space, and is the approximation that has been considered in most previous works on the subject, for example in the so-called νMSM. Although baryon asymmetries tend to be larger in the case of highly degenerate neutrinos, we find solutions with a very mild degeneracy that also correlate with a larger active-sterile mixing. These non-degenerate solutions appear for an inverted ordering of the light neutrinos. On the other hand we do not observe a direct correlation with other observables, such as the PMNS CP phases nor the neutrinoless double beta decay amplitude.
We have also performed a scan in the full parameter space, with the only requirement that one of the yukawa matrix eigenvalues is very small, and that one mode will not reach equilibrium before the electroweak transition, for the washout not to be complete. The main difference with the simpler case of two neutrinos is that the parameter space with successfull baryogenesis is significantly enlarged, in particular as regards non-degenerate spectra. Also the active-sterile mixings can reach larger values, particularly in the normal hierarchy case, improving the chances of future experiments such as SHiP or LBNE to find the GeV sterile neutrinos. There is much less difference in this case between normal and inverted neutrino orderings and also no direct correlation with the PMNS phases. On the other hand, the requirement of a small yukawa eigenvalue implies that the lightest neutrino mass cannot be large.
A number of refinements are needed to improve the precision of the determination of the baryon asymmetry. First a more precise determination of the scattering rates of the sterile neutrinos is required. Most previous studies, and this one, have included only topquark scatterings, but it has been pointed out recently that gauge scatterings are also very important. A correct treatment of these processes in the kinetic equations is necessary. Also the kinetic equations neglect effects of O((M i /T ) 2 ). Such effects are not so small for masses in the GeV near the electroweak phase transition and their effect should be quantified. Finally, spectator processes and the asymmetries of fields other than the sterile neutrinos and LH leptons have not been taken into account in the kinetic equations. A proper treatment could easily bring corrections of O(1). Finally, a more ambitious scan of parameter space should define more accurately the limits of eq. (5.6) for successfull baryogenesis. These effects will be studied in the future. We can factor out the α dependence and define: with

(B.4)
On the other hand, for the invariant I