Sterile neutrino dark matter via GeV-scale leptogenesis?

It has been proposed that in a part of the parameter space of the Standard Model completed by three generations of keV...GeV right-handed neutrinos, neutrino masses, dark matter, and baryon asymmetry can be accounted for simultaneously. Here we numerically solve the evolution equations describing the cosmology of this scenario in a 1+2 flavour situation at temperatures $T \le 5$ GeV, taking as initial conditions maximal lepton asymmetries produced dynamically at higher temperatures, and accounting for late entropy and lepton asymmetry production as the heavy flavours fall out of equilibrium and decay. For 7 keV dark matter mass and other parameters tuned favourably, $\sim 10\%$ of the observed abundance can be generated. Possibilities for increasing the abundance are enumerated.


Introduction
The idea of accounting for dark matter through keV-scale sterile neutrinos [1,2] is strongly constrained by now (for a review see, e.g., ref. [3]). The non-observation of γ-rays from putative sterile neutrino decays restricts their Yukawa couplings to be very small, |h Ia | < 10 −12 . With such small couplings a sufficient number of sterile neutrinos can be produced in the Early Universe only if the production is enhanced through a resonant mechanism [2], requiring the presence of large lepton asymmetries. Some time ago, it was pointed out [4] that this scenario could be embedded in a framework in which two generations of GeV-scale righthanded neutrinos first generate a baryon asymmetry [5,6], and then continue to generate lepton asymmetries, which subsequently boost dark matter production [7][8][9].
In the most detailed dark matter computation carried out so far [9], it was assumed that lepton asymmetries are produced first, at T > ∼ 5 GeV, whereas dark matter production is only active at T < ∼ 5 GeV. However, if the mass scale of the heavier sterile neutrinos is M > ∼ 2 GeV, they decay at T ≪ M/π < ∼ GeV, and these decays may produce further lepton asymmetries [4,10]. Then lepton asymmetry generation and dark matter production may proceed simultaneously, and need to be accounted for within a unified framework.
The purpose of the present paper is to assume that the initial lepton asymmetries have been dynamically produced by two generations of GeV-scale right-handed neutrinos. In a recent work [11], we showed that in this case lepton asymmetries > ∼ 10 3 times larger than the baryon asymmetry can arise. Furthermore, the lepton asymmetries have an intriguing structure, being evenly distributed amongst all flavours and settling into a stationary state (see also ref. [12]). We now follow that state down to lower temperatures, at which the GeV-scale right-handed neutrinos freeze out and decay. This non-equilibrium dynamics modifies the expansion rate of the universe and may also produce new lepton asymmetries. The question is whether this could help to boost the asymmetries that, according to ref. [11], were too small to have a substantial effect in the dark matter context. For the dark matter sector itself we fix the mass to the prototypical 7 keV scenario, with the corresponding Yukawa couplings pushed to the maximal allowed range as suggested by supposed observations [13,14]. 1 The presentation is organized as follows. The rate equations applying to the 1+2 sterile neutrino system are summarized in sec. 2. In sec. 3 we explain how the falling out of equilibrium of the "heavy" GeV-scale flavours modifies the expansion of the universe, and transcribe the rate equations to this situation. Subsequently, the heavy part of the rate equations can be simplified as explained in sec. 4, whereas the "light" keV-scale part may experience resonant enhancement, cf. sec. 5, which prohibits any substantial simplification. Parameter choices are justified in sec. 6, and numerical results are presented in sec. 7. A brief summary and outlook conclude this investigation in sec. 8.

Review of rate equations for the 1+2 flavour situation
The theory we work with is described by the Lagrangian where M I ≥ 0 are Majorana masses;φ = iσ 2 φ * is a Higgs doublet; a L , a R are chiral projectors; ℓ a = (ν e) T a is a left-handed lepton doublet of generation a; h Ia are the components of the neutrino Yukawa matrix; and summations over indices are left implicit.
We consider the situation M 1 ∼ keV ≪ M 2,3 ∼ GeV, and assume the 2nd and 3rd generations to be almost degenerate in mass. The average mass is denoted by M H ≡ (M 2 + M 3 )/2. The heavy flavours I = 2, 3 and the associated Yukawa couplings |h Ia | < ∼ 10 −7 are chosen to reproduce the active neutrino mass differences and mixing angles, whereas the first generation has much smaller Yukawa couplings, |h 1a | < 10 −12 , as is suitable for playing a role in dark matter physics.
The density matrix of the hierarchical 1+2 flavour system is expressed as and similarly for other objects. Here denotes a symmetrization/antisymmetrization with respect to helicity (±), and off-diagonal heavy-light components of ρ ± have averaged out up to effects suppressed by 1/M H . The lepton asymmetry in flavour a is denoted by n a , whereas n B is the baryon asymmetry. The evolution equation for lepton asymmetries can be split into the contributions of the light and heavy flavours, 2 (2.4) where n F denotes the Fermi distribution, k ≡ d 3 k (2π) 3 , and ω I ≡ k 2 + M 2 I . The light and heavy components of the density matrix evolve aṡ The coefficients associated with the light flavours read where (...) L indicates the use of a "light" mass M 1 , whereas the heavy coefficients read where (...) H stands for a "heavy" mass M H . We have denoted leptonic chemical potentials byμ a ≡ µ a /T , and expressed the dependence on neutrino Yukawa couplings through 14) whereas Q ± (a) = [Q (a+) ± Q (a−) ]/2 denote symmetrization and antisymmetrization with respect to helicity. The coefficients Q andQ parametrize the C-even and C-odd parts, respectively, of "absorptive" reactions (i.e. real processes), whereas U andŪ parametrize "dispersive" corrections. Specifically, where Π R a is a retarded correlator associated with the operator j a =φ † a L ℓ a to which the sterile neutrinos couple; τ = ± denotes helicity; I = L, H refers to the flavour; and Q andQ can be extracted by symmetrizing and antisymmetrizing in chemical potentials, respectively.
For a practical determination of Q,Q, U,Ū , we have generalized the computations of refs. [11,18,19] to arbitrary kinematics (i.e. not only the ultrarelativistic regime πT ≫ M but also πT ∼ M or πT ≪ M ), restricting however still to the approximation M ≪ m W in the treatment of 2 ↔ 2 scatterings below the electroweak crossover.
In order to close the system, the chemical potentials appearing in eqs. (2.7)-(2.15) need to be re-expressed in terms of the number densities appearing on the left-hand side of eq. (2.4). This requires the determination of "susceptibilities". We follow the approach in appendix A of ref. [11], simplifying the formulae by restricting to T 2 ≪ v 2 , where v ∼ 246 GeV, but adding charged lepton and light quark masses according to ref. [9]. Hadronic contributions are smoothly switched off at low T by a replacement N c → N c,eff , as proposed in ref. [18].

Non-equilibrium expansion
The GeV-scale flavours that are responsible for leptogenesis at T ∼ 130 GeV, freeze out and subsequently decay when πT ≪ M H . These non-equilibrium decays release entropy [20], an effect which has been argued to be substantial for M H ≃ 1...10 GeV [21], and which therefore needs to be included in dark matter and baryogenesis computations. (When πT ≫ M H , the GeV-scale flavours already have a small effect on the energy and entropy densities, however this is on the percent level and thus insignificant on our resolution.) Denoting by a(t) the cosmological scale factor and by m Pl = 1.22091×10 19 GeV the Planck mass, and assuming a flat universe, Friedmann equations can be expressed aṡ where e is the energy density, p is the pressure, and H is the Hubble rate. We write where e T and p T are the Standard Model energy density and pressure at a temperature T , whereas e H and p H represent the contribution of the heavy right-handed neutrinos. If we denote by a co-moving momentum mode and by kt ≡ 3 the corresponding phase space integral, the energy density and pressure carried by the heavy flavours can be expressed as We now insert eq. (3.3) into eq. (3.2), and move the thermal terms to the left-hand side. Making use of de T = T ds T and e T + p T = T s T , where s T is the Standard Model entropy density, we find In order to proceed, it is helpful to express the phase-space integrals in eq. (3.5) in terms of the time-independent variable k (cf. eq. (3.4)), because then ρ + II appears in a form for which a time-evolution equation is available. For a combination appearing in eq. (3.6) this implies A derivative with respect to t now operates on two terms, ρ + II as well as the last piece, Once p H from eq. (3.5) is inserted, the contribution from eq. (3.8) cancels against the contribution from p H in eq. (3.6). In total, then, At this point we make use of the equation of motion of ρ + II . It follows from eq. (2.6) that, to a good approximation, Inserting eq. (3.10) into eq. (3.9) and going subsequently back to co-moving momenta as integration variables, we get We see that entropy is generated only if ρ + II falls out of equilibrium. Let us simplify the setup by making use of so-called momentum averaging. Even though not associated with any formally small expansion parameter, this turns out to represent a reasonable approximation in many cases [11,22,23]. We integrate eq. (3.10) over k and then change variables into k t , which leads to Now introduce the ansatz where Y + II is a yield parameter, and denote . (3.14) Then eqs. (3.11) and (3.12) become As a final step, the evolutions of Y + II and s T a 3 can be decoupled from each other, by inserting eq. (3.15) into (3.16). Moreover, introducing where c 2 s is the speed of sound squared. From eqs. (3.15) and (3.18) the Jacobian J can be solved for, For M H = 0.2 GeV the process takes place at temperatures lower than those shown. Right: the (inverse) Jacobian from eq. (3.19), normalized to the value within the Standard Model [24], with . Two competing effects in J lead to a non-monotonous behaviour: an increase of H, and a decrease due to the non-equilibrium term placed just after 3H in eq. (3.19).
Then the basic equations become Numerical solutions for s T a 3 , obtained from eqs. (3.20) and (3.21), are shown in fig. 1(left). In fig. 1(right) we show the corresponding J from eq. (3.19), normalized to its Standard Model value. According to fig. 1(left), entropy release substantially reduces any yields that were generated at Even if only the evolution of Y + II is directly coupled with the evolutions of s T a 3 and J (cf. eqs. (3.20), (3.21)), the results in fig. 1 also influence other evolution equations. The redshift factor from eq. (3.4) can be expressed as As we have replaced t as an integration variable through x = ln(T max /T ) with the help of J , the co-moving momentum will from now on be denoted by k T . Defining Y ≡ n/s T , evolution equations for particle densities and phase space distributions from sec. 2 are transcribed aṡ

Simplified treatment of heavy flavours
It was indicated in the previous section that for the heavy flavours it is advantageous to resort to momentum averaging. Moreover, it is convenient to go over into an interaction picture. In order to simplify the notation of eqs. (3.23) and (3.24), let us denote Then the role of the free Hamiltonian is played by diag ω 2 , ω 3 1 − H + H 1 , where the averaging ... 1 is defined in eq. (3.14). From here we can subtract the trace part without loss of generality. The remaining upper diagonal appearing in eq. (2.6) is defined as 2) becomes large at low temperatures (recall that ω I are normalized by J , which decreases like the Hubble rate, as ∝ T 2 ). In this situation the fast oscillations between the heavy sterile neutrinos, induced by H fast , can be "integrated out". Working to leading order in 1/ H fast as described in ref. [23], we find that in this regime Here we have complemented the momentum average in eq. (3.14) through The helicity-symmetric diagonal components of the density matrix evolve according to eq. (3.20), whereas the other components obey

Resonant contribution in light flavour
The question arises whether momentum averaging could also be adopted for f ± . This is, however, hindered by the possible appearance of a "resonance" in the coefficients Q ± (a)L ,Q ± (a)L , which parametrize the evolution of f ± through eqs. (2.7)-(2.9). The resonance originates through the helicity-conserving indirect contribution, which for M 1 ≪ k T has the form . (5.1) The function b has a C-even and C-odd part; the latter, which is proportional to chemical potentials, is denoted by b| µ ≡ c. At low temperatures the C-even part is to a good approximation proportional to ω 1 [25]. Therefore we may write b =b ω 1 + c .

(5.2)
The functionb is positive, whereas c is odd in the interchange µ i ↔ −µ i . Therefore, after extracting the C-even Q (a−) and the C-oddQ (a−) from eq. (5.1) by symmetrizing and antisymmetrizing in chemical potentials, respectively, both contain one appearance of For small ω 1 Γ u this can be approximated as This is qualitatively different from non-resonant contributions, which are proportional to Γ u . We observe from eqs.

Parameter values and initial conditions
We start by considering the benchmark point from ref. [11], tuned to produce the observed baryon asymmetry as well as maximally large low-temperature lepton asymmetries, within a specific slice of the parameter space. The most important parameters are M H ≈ 0.7732 GeV, ∆M = 10 −11 GeV, Im z = −0.15, where the last one refers to the Casas-Ibarra parameter fixing the absolute value of the neutrino Yukawas [26]. The small | Im z| implies |h Ia | ≃ 2 × 10 −8 . The corresponding active-sterile mixings, |h Ia |v/( √ 2M I ) ≃ 4.5 × 10 −6 , are tiny and thus challenging to constrain in (future) experiments.
The initial conditions for the evolution are set at a temperature T ≈ 5 GeV, where the system is to a good approximation in a stationary state [11]. Taking also into account that rate coefficients are dominated by helicity-conserving contributions at low temperatures, eqs. (4.3)-(4.6) imply that whereμ ave ≡ 1 3 aμa and X − eq ≡ 1 2 . To be optimistic, we multiply lepton asymmetries obtained in ref. [11] by a factor two, leading to the initial condition Y a − Y B /3 ≃ −6.2 × 10 −7 for all a, which fixes the chemical potential appearing in eq. (6.1) asμ ave ≃ −6.6 × 10 −5 . The initial baryon asymmetry is set at the observed value Y B = 0.87 × 10 −10 ; in view of entropy dilution, it should be taken to be somewhat larger at the beginning, but this has little effect on our considerations here, and is also easy to achieve in practice [11].
In addition to the benchmark point , we have carried out further scans like in ref. [11], and again multiplied the corresponding lepton asymmetries by a factor two. This leads to two further parameter points which serve to illustrate the dependence on M H . Initial lepton asymmetries can be kept large by decreasing M H , but there is not much room here, given that according to refs. [27,28]

Numerical solution
Important ingredients characterizing the solution of the rate equations are equilibration rates, which determine how efficiently different components of the density matrix approach their would-be equilibrium values. As an example, consider the dimensionless combination appearing in eq. (4.5) but for simplicity normalized to the thermal Hubble rate rather than J , The result is shown in fig. 2(left). We observe that the system can follow equilibrium (i.e. thatΓ H > ∼ 1) when T > ∼ 2 GeV, but at T < 2 GeV there is a period when this should not happen. At very low temperatures, rates are dominated by vacuum decays, and the system again approaches equilibrium. For the light flavour, we show an equilibration rate from eqs. (2.5), (2.9), evaluated at a fixed comoving momentum, in fig. 2(right) (Γ L ≡ 2 a φ + (a)11 Q + (a)L ). Given that in the full rangeΓ L ≡ Γ L /(3c 2 s H T ) ≪ 10 −3 , the light flavour never comes near thermal equilibrium. For benchmark (i.e. M H ≈ 0.8 GeV), the solutions of the rate equations for lepton asymmetries and the density matrix of the heavy sector are shown in fig. 3. At first the density matrix follows the equilibrium form, but at T ≪ 1 GeV the equilibrium form starts  Middle: helicity-symmetric components of the density matrix, compared with the equilibrium value Y + eq . Right: the fraction of dark matter that Y + 11 ≡ k T f + /s T accounts for, cf. eq. (7.2). The decrease at low T is due to entropy dilution.  to decrease as mass effects become important. The actual solution cannot immediately follow this change, given that the equilibration rate has become small.
We note from fig. 3 that even if the density matrix deviates from equilibrium at low temperatures, there is no substantial re-generation of lepton asymmetries taking place in this regime. The reason is that the rate coefficients are so small that the source terms, cf. the last lines of eq. (4.3), remain inefficient.
Making us of Ω dm h 2 ≈ 0.12 [29] and ρ cr /[h 2 s(T 0 )] = 3.65 eV [30], where s(T 0 ) is the current entropy density, the fraction of dark matter carried by the lightest right-handed neutrinos Due to a smaller mass, the right-handed neutrinos do not decay efficiently, which leads to the problem that they may carry too much energy density at late times (cf., e.g., refs. [27,28]). Dark matter abundance is in the same ballpark as in fig. 3, however entropy dilution has not started yet.
can be expressed as We observe from fig. 3(right) that intermittently about 8.5% of the total dark matter abundance could be accounted for, before entropy dilution kicks in at late times. The yields of the helicity asymmetries are illustrated in fig. 4. Helicity asymmetries remain modest (Y − II ≪ Y + II , I ∈ {1, 2, 3}), which is a manifestation of the fact that thermal production dominates over resonant production (one helicity state is produced from neutrinos, the other from antineutrinos). The dark matter phase space spectra, which are strongly tilted towards the IR compared with kinetically equilibrated fermions, are shown in fig. 4(right).

Summary and outlook
The purpose of this paper has been to update our previous sterile neutrino dark matter analysis [9] by fixing the initial lepton asymmetries to maximal values that can be produced by the dynamics of GeV-scale right-handed neutrinos [11]. The parameters of the latter are constrained to be responsible for generating the active neutrino masses and mixing angles [26]. We permit for the generation of further lepton asymmetries in the low-temperature decays In this case the initial lepton asymmetries obtainedà la ref. [11] are small, but novel asymmetries are generated while Y + 22 , Y + 33 are out of equilibrium (the suppression by H fast in eq. (4.3) is originally moderate in this case, ∼ 1/500). However there is not much effect on the dark matter abundance. . The ratio Ω 1 /Ω dm at T = 1 MeV as a function of the initial lepton asymmetry. Some further entropy dilution is expected at T < 1 MeV, particularly for M H = 0.2 GeV, but the cosmological background is simultaneously becoming more complicated, as active neutrinos decouple and big bang nucleosynthesis starts. In any case, obtaining Ω 1 ≈ Ω dm would require lepton asymmetries about two orders of magnitude larger than those found in ref. [11]. of the right-handed neutrinos, 4 by including both the light and heavy sterile flavours in the set of rate equations, and track the modification of the universe expansion caused by the energy density carried and entropy released by the heavy flavours. In addition we resolve both helicity states of the sterile neutrinos; this is important for heavy flavours given that initial lepton asymmetries are correlated with helicity asymmetries [11,12], and for the light flavour given that resonant production [2] affects one helicity state only. Even though we do find rich dynamics in the heavy sector (cf. figs. 3, 5, 6), the dark matter abundance does not vary greatly between the cases, reaching typically less than 10% of the observed value. The reason for this behaviour can be understood as follows. The dependence of Y + 11 on lepton asymmetries must be quadratic at small Y a , given that energy density is a C-even quantity. A strongly growing dependence only sets in at |Y a | ≫ 10 −6 , cf. fig. 7. This is associated with the dominance of resonant production, which in the language of eq. (5.4) requires c 2 > 2bM 2 1 . Given that the asymmetries obtained in ref. [11] are below this level, dark matter production takes place predominantly through normal thermal processes. Therefore, our dark matter results are rather insensitive to the heavy sector, apart from its influence through the expansion of the universe, as depicted in fig. 1. In order to account for 100% of dark matter, initial lepton asymmetries should be a factor ∼ 10 2 larger than those found in ref. [11], i.e. of the order |Y a | ∼ (2...10) × 10 −5 , depending on entropy dilution.
Even if we have failed to account for all of dark matter through the dynamics of keV...GeV scale sterile neutrinos, the results could change in the future, for several reasons (in line with the basic minimalistic premise of our study, we keep the Lagrangian of eq. (2.1) intact for this discussion, without any extra non-SM fields, and assume a standard cosmology): (i) As the dark matter production is largely thermal rather than resonant, it is proportional to the rate coefficient Γ u , which contains large hadronic uncertainties [31]. It would be interesting to estimate or constrain Γ u through lattice simulations.
(ii) We have chosen the Yukawa couplings of the heavy flavours to be as small as possible, in order to diminish lepton number washout and therefore to have maximal initial asymmetries [11]. However, as the initial asymmetries have little influence in any case, the Yukawas could be made larger, without spoiling baryogenesis (cf., e.g., refs. [11,32,33]). Then the heavy flavours would stay closer to equilibrium and decay faster, producing less entropy. Even though we do not expect substantial variations of the dark matter abundance from here, a comprehensive study of the heavy flavour Yukawas would be welcome. This should also include the search for potential "atypical" CP phases where late-time lepton asymmetries might be anomalously large.
(iii) We have restricted ourselves to the SHiP window, M H < M B ≃ 5 GeV [34], but nature may have chosen otherwise. Increasing M H in the analysis of ref. [11], we find that initial lepton asymmetries would be smaller then. However, as anticipated in refs. [4,10], novel lepton asymmetries are produced later on (cf. fig. 6). Therefore it seems promising to explore what happens with larger values of M H . Then, however, 2 ↔ 2 scatterings entering the rate coefficients need to be addressed without resorting to the approximation M H ≪ m W , which poses a significant technical challenge. The initial temperature should be chosen in the regime T ≫ M H , i.e. larger than here.
(v) Last but not least, the observational status of the dark matter sterile neutrino remains unclear. Here we have relied on the indications in refs. [13,14], however other parameter values could be studied within our framework, and might change the conclusions.
To summarize, we have established a framework which permits to study sterile neutrino dark matter production in a non-degenerate 1+2 flavour situation, with ongoing lepton number violation and with the heavy flavours falling out of equilibrium and gradually decaying. As a proof of concept, we have excluded several SHiP-like benchmarks as an explanation for all of dark matter. Broader parameter scans may help to bridge the gap.