Fermion number violating effects in low scale leptogenesis

The existence of baryon asymmetry and dark matter in the Universe may be related to CP-violating reactions of three heavy neutral leptons (HNLs) with masses well below the Fermi scale. The dynamical description of the lepton asymmetry generation, which is the key ingredient of baryogenesis and of dark matter production, is quite complicated due to the presence of many different relaxation time scales and the necessity to include quantum-mechanical coherent effects in HNL oscillations. We derive kinetic equations accounting for fermion number violating effects missed so far and identify one of the domains of HNL masses that can potentially lead to large lepton asymmetry generation boosting the sterile neutrino dark matter production.


Introduction
Though the canonical Standard Model (SM) has been completed by the discovery of the Higgs boson and may be a valid effective quantum field theory all the way up to the Planck scale (for recent discussions see [1][2][3]) it is inconsistent with a number of observations. They include the non-zero neutrino masses, the presence of Dark Matter (DM) in the Universe, and its baryon asymmetry (BAU). Perhaps, the most minimal way to address all these problems on the same footing is to extend the SM by three right-handed neutrinos with masses below the Fermi scale [4,5]. These new fermions N I , I=1,2,3 (following the Particle Data Group [6] we will call them Heavy Neutral Leptons or HNLs for short) are singlets with respect to the SM gauge group and thus are allowed to have Majorana neutrino masses. The lightest of these particles, N 1 , may play a role of Dark Matter [7,8]. Two others (N 2 and N 3 ), if (almost) degenerate, can produce the baryon asymmetry of the Universe [9], [5] and explain non-zero neutrino masses and mixings at the same time. This model was dubbed the νMSM for "Neutrino Minimal Standard Model" [4]. For a number of computations of baryon asymmetry in this model see [10][11][12][13][14][15][16][17][18][19][20].
The most conservative scenario of the Universe evolution, which does not require any new physics beyond the νMSM, proceeds as follows. First, the Universe is inflated by the SM Higgs field [21] and heated up due to Higgs field oscillations to temperatures T ∼ 10 14 GeV [22][23][24]. The Higgs inflation prepares the initial conditions for the Hot Big Bang [22] at T ∼ 10 14 GeV: baryon and lepton numbers of the Universe are equal to zero, and the number densities of HNLs at this time are zero as well. The Email addresses: Shintaro.Eijima@epfl.ch (Shintaro Eijima), Mikhail.Shaposhnikov@epfl.ch (Mikhail Shaposhnikov) particles N 2 and N 3 enter into thermal equilibrium below the sphaleron freeze-out temperature T sph ≃ 130 GeV and produce baryon asymmetry of the Universe in a set of processes which include their coherent oscillations, transfer of lepton number from HNLs to active leptons and back [9], [5], and rapid anomalous sphaleron transitions [25]. The lighter HNL -N 1 -DM sterile neutrino never equilibrates and is mainly produced at temperatures T DM ∼ 100 − 300 MeV by transitions from the ordinary neutrinos to N 1 [7,8,[26][27][28][29][30]. The combination of X-ray and Lyman-α bounds on the DM sterile neutrino excludes the "nonresonant" Dodelson-Widrow mechanism [7] for their production, which operates in the cosmic plasma with small lepton asymmetries. In other words, to get enough DM particles N 1 , the processes involving N 2,3 should produce [11] sufficiently large lepton asymmetry ∆L/L > 2 × 10 −3 which must be present at temperatures T DM . This is needed to boost the production of N 1 due to the resonant mechanism proposed by Shi and Fuller in [8] and developed in a rigorous way in [27][28][29][30]. The production of this large lepton asymmetry must take place below the sphaleron temperature T sph , otherwise the baryon asymmetry will be too large [11].
The estimates of the equilibration rates of N 2,3 in [11,15,16] and in more recent works [31,32] based on careful thermal field theory computations showed that for all parameter choices consistent with observed pattern of neutrino masses and oscillations the HNLs N 2,3 enter in thermal equilibrium at some temperature T in exceeding tens of GeV and go out of thermal equilibrium at temperatures T out < T in which can be as small as 1 GeV. This has led to the conclusion that the equilibrium period between T in and T out erases all the lepton asymmetry which could have been generated at freeze-in temperature T in , requiring that the large lepton asymmetry needed for effective dark matter production must be created at T < T out .
The analysis made in [16] demonstrated that a large lepton asymmetry can indeed be generated in the scattering processes involving N 2,3 at the freeze-out temperature T out and below it in out-of-equilibrium decays of N 2,3 . This asymmetry does not exceed ∆L/L ≃ 3 × 10 −2 , leading to the conclusion that the mass of the DM sterile neutrino must lie in the interval from 1 to 50 keV, to be consistent with the Lyman-α and phase density constraints coming from observations of dwarf galaxies [33,34]. As for N 2,3 , their physical masses should be between 1.5 GeV and ≃ 80 GeV (the W -boson mass) and be extremely degenerate, ∆M phys /M < 10 −15 [11,15,16,35,36] 1 . The latter condition comes from the requirement that the period of N 2 ↔ N 3 oscillations should be comparable with the age of the Universe at the time of lepton asymmetry production, to insure the resonance [11,35]. The minimal scenario that has been proven to work, albeit under the requirement of a strong fine-tuning (of the order of 10 −4 [16]) between two different contributions to the physical mass difference: one coming from the Yukawa couplings and the Higgs condensate, and another from Majorana masses of N 2,3 .
The aim of the present paper is to show that the part of the lepton asymmetry generated at T in can in fact survive until the temperatures of sterile neutrino DM production ∼ 100 MeV, in-spite of the fact that HNLs are well in thermal equilibrium between T in and T out . Qualitatively, this comes about because of the following reasons. In the symmetric phase of the electroweak theory the transfer of asymmetry from active to sterile sector and back occurs mainly via the processes with fermion number conservation (we attribute positive fermion number to left-handed neutrinos and to right-handed HNLs) with the rate Γ + . The rate of fermion number non-conserving processes Γ − is suppressed by a kinematic factor (M/k) 2 , where M is the HNL mass, and k ∼ 3T is the typical momentum of fermions in the plasma.
On the contrary, in the Higgs phase, at temperatures of the order of tens GeV, the dominant reaction is induced by the mixing term between ν's and N 's and has a rate Γ − exceeding that of the Universe expansion at T out < T < T in . It proceeds with fermion number non-conservation: left handed neutrinos go into left-handed anti-HNLs and vice-versa.
In a large portion of the νMSM parameters the reactions with fermion number conservations are faster and give the main contribution to baryogenesis at the sphaleron freeze-out temperature T ≃ 130 GeV. However, in a specific domain of the NHL masses and couplings, the rate of these processes never exceeds the Hubble rate. Thus, the asymmetry in this almost conserved number is protected from dilution, in-spite of the fact that HNLs are equilibrated due to the processes with fermion number vi-1 These restrictions are not required if the DM sterile neutrino is produced by some new interactions not contained in the renormalizable νMSM Lagrangian [37][38][39]. olation 2 . Moreover, at T ≃ T in the rate Γ + can be close to the Hubble rate, meaning that large asymmetry in this number can be generated. To understand whether it is indeed produced would require numerical solution of our integro-differential kinetic equations for many parameters of the νMSM, which is not attempted here.
The paper is organised as follows. In Section 2 we will derive kinetic equations accounting for helicity structure of HNL interactions. In Section 3 we analyse the different rates and identify the range of νMSM parameters which may potentially lead to large lepton asymmetries surviving until small temperatures where the production of DM sterile neutrino takes place. In Section 4 we summarise our results.

Kinetic equations and helicity
In this section we derive the kinetic equations describing the evolution of HNLs density matrix and lepton densities. To elucidate their structure, we will consider first the temperatures well below the electroweak scale, deeply in the Higgs phase. We also neglect for the time being the subtleties related to electric neutrality of the plasma [41][42][43][44] and consider the system with HNLs and active neutrinos only, this will be corrected towards the end of this Section. We take a pair of almost degenerate HNLs N 2 and N 3 , the generalisation to the case of arbitrary number of N is straightforward. It is convenient to unify N 2 and N 3 in one Dirac spinor, as has been done in [11], and consider the νMSM Lagrangian where L SM is the SM part, Ψ = N 2 + N c 3 is the HNL field in the pseudo-Dirac basis, M = (M 3 + M 2 )/2 and ∆M = (M 3 − M 2 )/2 are the common mass and Majorana mass difference of HNLs, respectively, h αI is a matrix of Yukawa coupling constants and Φ is the temperature dependent Higgs vacuum expectation value 3 , which is 174.1 GeV at zero temperature. The HNL field Ψ is given in terms of creation and annihilation operators by plus, particles minus, anti-particles where σ = ± describes the HNL helicity. The attribution of fermion numbers to a σ , b σ leading to fermion number conservation in the limit M → 0, ∆M → 0 is shown in Table 1.
We will consider L int as a perturbation and work in the second order in Yukawa couplings and first order in ∆M assuming for power counting that M ∆M ∼ h 2 Φ 2 , where h is a typical value of the Yukawa couplings. The quadratic mixing term L int leads to communication between sterile sector of HNLs and the rest of the SM, ensuring creation and destruction of HNLs, their coherent oscillations, lepton number non-conservation, and transfer of asymmetries from active flavours to sterile and back.
To construct kinetic equations, we work in the Heisenberg picture of quantum mechanics and use the ideas of [45] in what follows 4 . The derivation is performed by using four creation operators a † σ (k), b † σ (k) for HNLs and two (for each generation α) neutrino and antineutrino operators, . Let ρ be a (time-independent) density matrix of the complete system. The HNL abundances, including coherent quantum-mechanical correlations between N 2 and N 3 , are given by the averages The neutrino number densities are a † να (k)a ν β (k) and b † να (k)b ν β (k). Symbolically, all these number-density operators will be denoted by Q 0 i . With the use of our generic notation Q 0 i , the abundances are given by q 0 The time derivative of q 0 i is readily found: where H is the total Hamiltonian of the system under consideration. Since we neglected the existence of charged leptons for the moment, the commutator of Q 0 i with the SM Hamiltonian is zero, the only non-trivial contribution comes from the interaction Hamiltonian H int , associated with L int defined in (1). It is easy to see that the commutators [H int , Q 0 i ] are the quadratic polynomials with respect to creation and annihilation operators, containing all possible terms, to list just a few: . These operators are multiplied by the first power of Yukawa couplings or by ∆M . Denote by Q 1 i these binomials, and by q 1 i = Tr[Q 1 i ρ]. Now we continue the chain, and write an equation for every q 1 i which appeared at the first step: and then repeat the procedure again and again. In this way we get an infinite chain of kinetic equations, which includes the averages of higher and higher polynomials in creation and annihilation operators. To truncate the system, we proceed as follows. The total Hamiltonian of the system can be written as where H 2 is the quadratic part including HNLs and active neutrinos, and H SM int describes the SM interactions. As an example, let us take an operator . Its commutator with the total Hamiltonian contains 3 terms. The first one is where E N (k) = √ k 2 + M 2 and ǫ ν (k) = k are the energies of HNL and active neutrino respectively (the small neutrino mass can be safely neglected here). The sec- and is of the order of Fermi constant G F and of 4th order in creation and annihilation operators. It accounts for neutrino interaction in the medium. The first two terms contain the operators that have already showed up in the first kinetic equation (3), while the third one contains new operators. To find their time evolution would require next steps in the iterative procedure.
To deal with the third term, we note that the neutrino interactions in the plasma can be accounted for by modification of neutrino energy ǫ ν (k), replacing it with the temperature dependent dispersion relation E ν (k) (related to the real part of neutrino propagator Σ) and by attributing to it an imaginary part γ ν (k) > 0 (associated with absorptive part of Σ). These considerations suggest the following modification of the commutation relations: [a † να (k), where the signs in front of γ ν (k) are chosen in such a way that they correspond to damping rather than an instability. The similar rules apply to antineutrinos. These substitutions effectively account for the third term which can now be removed. The system of kinetic equations for q 0 i and q 1 i is now complete. It can be simplified even further as the active neutrinos are well in thermal equilibrium at all temperatures we are interested it, γ ν /H ≫ 1, where H is the Hubble rate. Again, we take Q 1 1 as an example. The equation for it now reads and has an approximate slow varying solution All q 1 i can be found in this way and inserted into equation (3) for q 0 i . As a result, we get the kinetic description in terms of q 0 i only. The realisation of this program requires a straightforward but tedious computation, which we have done with the use of DiracQ Mathematica package [46]. By introducing the notations where ρ eq ν and ρ eq N are equilibrium distribution functions of neutrinos and HNLs, and "1" is the unity matrix, we arrived to the following result 5 The effective Hamiltonian is where σ 1 is the Pauli matrix. The production rates for HNLs are 5 We do not account here the expansion of the universe, but it can be easily accommodated.  and those for active neutrinos are The coefficients are given by and the matrices of Yukawa coupling constants are where E ν = k − b L and function b L is often called neutrino potential in the medium. It has been computed in a number of papers in different limits [47,48]. The neutrino damping rate as well as b L can be taken from a recent work [32]. The parts of the Hamiltonian and production rates with subscript "+" and "−" are associated with the fermion number conserving and violating operators, respectively. In Figure 1 they are expressed diagrammatically, where the vertexes (denoted by the cross) in (a) come from the structures in H int containing the product of two creation (or annihilation) operators, as for example in (36) and those in (b) coming from a product of creation and annihilation operators as, for instance: The structure of the kinetic equations is exactly the same as it was first found in [5] and elucidated in [11] 6 , but now with all the terms expressed explicitly through parameters of the theory. The same set of equations (15)(16)(17)(18) was used in [10,[12][13][14][15][16][17][18][19][20], for analysis of baryon asymmetry generation in the νMSM, but with different choices of kinetic coefficients. The novel results of our work are the formulas for the Hamiltonian and different rates (28)(29)(30)(31)(32)(33) neatly separating the effects of the processes with and without fermion number non-conservation, the goal which was attempted already in [11] but not achieved at that time.
Now, we make these equations more realistic, accounting for the presence of charged fermions in the plasma, equilibrium character of electroweak reactions, and eventually sphaleron transitions. For this end we introduce leptonic numbers ∆L α of every generation, being a sum of asymmetries in neutrinos and charged leptons, integrated over momentum, and consider ∆ α = ∆L α −1/3∆B, where ∆B is the baryon asymmetry. Due to weak interactions the asymmetries in neutrinos are rapidly transferred to charge leptons and distributed among different momenta, whereas when sphalerons are operating there is a rapid transfer of lepton number to baryons. However, the rate of ∆ α change is proportional to HNL Yukawa couplings and corresponds to a slow process. A well defined procedure accounting for equilibrium character of weak reactions and electric neutrality of the plasma [41][42][43][44] allows 6 The kinetic equations derived in the seminal work [9] unfortunately are not correct as they do not contain the transfer (last) terms in Eqs. (17,18) and the equations for lepton asymmetries (15,16).
to write kinetic equations for ∆ α instead of ρ να and ρν α : where µ α are the chemical potentials for ∆ α and f f = 1/(e k/T +1) is the distribution function for massless fermion.
In the Higgs phase for m b T m W , where m b and m W are the masses of b-quark and W boson, the chemical potentials are given by where ∆B 0 is the freeze-out baryon number and ω αβ is the susceptibility matrix. Note that ∆ α are the momentumindependent quantities (neutrinos and charged leptons of a given generation are well in thermal equilibrium and thus their number densities are described well by the Fermi distribution with the chemical potentials for lepton numbers), but the density matrix for HNLs does depend on k. Yet another effect coming from charged leptons is the flavour dependence of the neutrino dispersion relations and neutrino damping rates. It can be neglected for temperatures exceeding the τ -lepton mass.
These completes the discussion of kinetic equations for HNLs and active flavours deeply in the Higgs phase, where the main contribution to active-sterile transition comes from the mixing terms contained in L int .
At the temperatures in the region of the electroweak crossover the structure of evolution equations remains the same, but a number of kinetic coefficients has to be modified. In particular, one has to add the "direct processes" (in terminology of Ref. [32]) involving HNLs and active neutrinos. These processes occur with fermion number conservation and their contribution does not vanish when Φ → 0. In fact, they dominate in baryogenesis around the sphaleron freeze-out in a part of the νMSM parameter space. The account for them results in the following modifications: , (42) where γ direct and Im Π R is the rate of the direct production of HNLs in the symmetric phase, which comes mainly from 2 ↔ 2 interactions. This contribution to γ + dies out in the Higgs phase; this is accounted for by a function K(m h ) [32]. The contribution to γ + from the Higgs decay to N ν is given by γ ph [32], f b (ǫ) = 1/(e ǫ/T − 1) is the bosonic distribution function.
In the symmetric phase these expressions simplify a lot: where h + is nothing but the Weldon high-temperature correction [49]. The last point is the modification of the susceptibility matrix in the symmetric phase, where sphalerons are in thermal equilibrium 7 : This formula is extracted from [41]. For determination of precise temperature dependence of ω αβ over the electroweak crossover see Ref. [32]. The relation between Ξ defined in Ref. [32] and our ω αβ reads ω αβ = 1 6T 2 Ξ −1 . The set of equations derived in this section allows to follow the system from very high to sufficiently small temperatures T ≃ 1 GeV, and address both baryon asymmetry generation around 130 GeV and late time lepton asymmetry production. At even smaller temperature one should take into account the flavour dependence of neutrino degrees of freedom in the medium, and include the decays and inverse decays of HNL which were omitted in our equations.

Thermal equilibrium and approximately conserved numbers
The kinetic equations (15)(16)(17)(18) allow to address the question of existence of approximately conserved quantum numbers. Let us consider two combinations of the HNL and active flavour asymmetries, where Making the corresponding linear combinations of Eqs. (15)(16)(17)(18), we find: The remarkable property of these relations is that the rates of the L ± change is proportional to the corresponding γ ± functions, which have very different behaviours as a function of temperature. In particular, the rate associated with γ + may never come into thermal equilibrium. We plot the ratio of the different HNL production rates to the Hubble rate H in Figures 2 and 3 for the normal hierarchy (NH) case and in Figures 4 and 5 for the inverted hierarchy (IH). The Γ ± rates defined in Eqs. (23,24) are 2 × 2 matrices, so we diagonalise each of them at any given temperature and denote by Γ ±,i with i = 1, 2 marking the corresponding eigen-values. In addition, we average the rates over momentum with the use of the Fermi distribution function, f F , where n F is the fermion number density. The rates depend on quite a number of parameters, the most important being the HNL mass M and the imaginary part of a complex mixing angle ω appearing in the Casas-Ibarra parametrisation of the HNL-neutrino mixing matrix [51]. The quantity X ω ≡ exp(Im ω) shown in the figures is a free parameter which does not change the active neutrino masses (we fix them with the available neutrino data). The amplitude of Yukawa couplings scales as F αI ∝ X ω for large Imω > 0 or F αI ∝ X −1 ω for Imω < 0. The parameter X ω is related to the introduced previously ǫ in [10,11] as ǫ = 1/X 2 ω . More specifically, the relation between Yukawa couplings and neutrino parameters, when ∆M is negligibly small, is given by where m 3 and m 2 are the heaviest and second-heaviest active neutrino masses and v = 174.1 GeV is the Higgs expectation value at zero temperature. It is the largest eigenvalue of Γ ± which determines the approach to thermal equilibrium for quantum numbers L ± . In general, the equilibration rates are larger for larger HNL masses, and smaller if X ω is close to one. In Figures 2 and 4 we show the dependence of the rates on temperature for different HNL masses.
In Figures 3 and 5 the behaviour of the rates as a function of temperature for different values of X ω is shown. The large (or small) X ω increases the magnitude of Yukawa couplings, leading to faster equilibration.
Let us comment on the splitting between the eigenvalues of the rates. From Eqs. (57-60) for large (or small) X ω the hierarchy between the diagonal components of the rates (23,24,32,33) gets large, which is reflected in the splitting of eigen-values. Note that the off-diagonal components of them are independent on X ω . The splitting is clearly visible in the NH case even if X ω ≃ 1 due to the presence of non-diagonal terms (59,60) which are of the order of atmospheric neutrino mass scale m atm ≃ 5 × 10 −2 eV. For the IH case the eigen-values of the rates at X ω ≃ 1 are almost degenerate due to the smallness of the off-diagonal components suppressed by (m sol /m atm ) 2 where the solar neutrino mass scale is m sol ≃ 9 × 10 −3 eV.
In the symmetric phase at T > 160 GeV, Γ − = 0 and the lepton number L − is conserved. The generated asymmetry in the active neutrino sector is the same as that in the HNL sector with an opposite sign. In the temperature range relevant for baryogenesis, i.e. above the sphaleron freeze-out T ≃ 130 GeV but below T = 160 GeV both Γ − and Γ + rates are present. However, the inspection of the figures shows that the rate Γ + dominates over Γ − , meaning that the discussion above is approximately valid, with small corrections from "-" contributions. As a result the previous computations of the baryon asymmetry performed in [5,[10][11][12][13][14][15][16][17][18][19][20] neglecting Γ − are legitimate.
When the temperature goes down below T ≃ 130 GeV the generation of the baryon asymmetry stops but of the lepton asymmetry continues [11]. Eventually, the rate Γ − starts to dominate over Γ + . The HNLs enter in thermal equilibrium at temperature T in corresponding to the intersection of the largest rate with the horizontal line Γ/H = 1, and go out of thermal equilibrium at T = T out < T in found in a similar way. Typically, for the set of masses considered, the highest rate is achieved at T ≃ 10 − 20 GeV. The maximum of the rate Γ − always exceeds the rate of the Universe expansion, as has been already demonstrated in [11,16,31,32].
The most interesting range of parameters that can potentially lead to generation of large leptonic asymmetries at T ≃ T in corresponds to relatively small HNL masses and X ω ≃ 1. Indeed, the Figures 2-5 show that the maximum of the ratio Γ +,2 /H does not exceed O(1) for M 2 GeV and X ω ≃ 1. At the same time, this ratio is close to 1 for these parameters, meaning that the large asymmetry in L + can potentially be generated but protected from washout until small temperature. The character of equilibrium which is led by the "-" reactions at smaller temperatures but has an (effective) L + conservation will ensure the relation between the asymmetries in active flavours and HNLs given by 8 Note that only the asymmetries in active flavours contribute to the resonant DM production. We leave the analysis of produced lepton asymmetries for a future publication.

Conclusions
In this paper we have derived kinetic equations for low scale leptogenesis, accounting for helicity effects in HNL interactions with active flavours. They are valid both for high temperatures around the sphaleron freeze-out (to describe the baryogenesis), and for smaller temperatures to account for lepton number generation. We showed that the total lepton number can be considered as approximately conserved in the certain domain of νMSM parameters, and thus can potentially be generated at T ≃ T in and survive until small temperatures making the resonant production of DM possible. It is still remains to be seen whether sufficiently large lepton asymmetry is indeed generated, as these depends on CP-violating effects. The work in this direction is in progress.
Quite remarkably, the requirement of existence of approximately conserved leptonic number makes the region of small HNL masses preferable, allowing to search for HNLs at SHiP or SHiP -like experiments. At the same time, it lies close to the lower bound coming from the re-