The new νMSM (ννMSM): radiative neutrino masses, keV-scale dark matter and viable leptogenesis with sub-TeV new physics

We study the scenario in which the Standard model is augmented by three generations of right-handed neutrinos and a scalar doublet. The newly introduced fields share an odd charge under a ℤ2 parity symmetry. This model, commonly known as “Scotogenic”, was designed to provide a mechanism for active neutrino mass generation as well as a viable dark matter candidate. In this paper we consider a scenario in which the dark matter particle is at the keV-scale. Such particle is free from X-ray limits due to the unbroken parity symmetry that forbids the mixing between active and right-handed neutrinos. The active neutrino masses are radiatively generated from the new scalars and the two heavier right-handed states with ∼ O\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{O} $$\end{document}(100) GeV masses. These heavy fermions can produce the observed baryon asymmetry of the Universe through the combination of Akhmedov-Rubakov-Smirnov mechanism and recently proposed scalar decays. To the best of our knowledge, this is the first time that these two mechanisms are shown to be successful in any radiative model. We identify the parameter space where the successful leptogenesis is compatible with the observed abundance of dark matter as well as the measurements from the neutrino oscillation experiments. Interestingly, combining dark matter production and successful leptogenesis gives rise to strict limits from big bang nucleosynthesis which do not allow the mass of dark matter to lie above ∼ 10 keV, providing a phenomenological hint for considered low-scale dark matter. By featuring the keV-scale dark matter free from stringent X-ray limits, successful baryon asymmetry generation and non-zero active neutrino masses, the model is a direct analogue to the νMSM model proposed by Asaka, Blanchet and Shaposhnikov. Therefore we dub the presented framework as “The new νMSM” abbreviated as ννMSM.


Introduction
The Standard model (SM) is a remarkably accurate theory. Its particle content and interactions are almost flawlessly mapping Nature's choice. However, there are still several open questions and in particular those that stand out are: • What is the origin of neutrino mass?
• Does dark matter (DM) interact with SM particles?
• Through which mechanism was the observed asymmetry between matter and antimatter generated?
The indisputable answer to any of these questions would indicate a tremendous progress for physics, in particular for the community striving to discover "new physics", as it is by now clear that the solution does not lie within the SM.
We attempt to address all of the above questions within the model proposed by Ma [1], dubbed "Scotogenic", in which the SM field content is supplemented by three right-handed neutrinos and a scalar doublet, all odd under a postulated Z 2 parity symmetry. The model was originally envisioned to account for small active neutrino masses and the electroweak scale DM produced via standard thermal freeze-out. Instead, we put in focus keV-scale for JHEP08(2018)067 the mass of the fermionic DM particle. We consider two viable regimes for DM production and ensure that the studied parameter space is fully consistent with leptonic mixing parameters and the observed neutrino mass squared differences.
In addition, we scrutinize the generation of BAU (baryon asymmetry of the Universe) from the produced lepton number asymmetry (leptogenesis). Such asymmetry transfer is viable due to the existence of non-perturbative sphaleron processes at high temperatures. We consider two complementary mechanisms, namely BAU production from right-handed neutrino oscillations introduced by Akhmedov, Rubakov and Smirnov (ARS) [2] as well as recently proposed asymmetry generation from scalar decays [3].
The motivation for this work stems from the νMSM model proposed by Asaka, Blanchet and Shaposhnikov [4][5][6]. In the νMSM, the SM particle content is extended with only three right-handed neutrinos, where the lightest one is a keV-scale DM produced via neutrino oscillations [7,8] due to the mixing between active and right-handed neutrinos. The heavier two GeV-scale right-handed states generate active neutrino masses in the seesaw type-I model [9][10][11][12] as well as produce BAU through the ARS mechanism. When confronted with the current experimental data, the νMSM is seriously challenged. In particular, the vast portion of the viable parameter space for DM is excluded by the combination of structure formation and X-ray limits [13][14][15].
A keV-scale DM candidate and ARS mechanism for baryogenesis are also prominent characteristics in our scenario. Furthermore, the relevant fermionic Yukawa and mass terms in these two models differ only in the employed scalar doublet -the SM Higgs doublet is considered in νMSM whereas our model hinges on the existence of a second Z 2 -odd scalar doublet. In contrast with νMSM, in the Scotogenic setup keV-scale DM does not mix with the active sector due to the imposed Z 2 symmetry. Hence, the DM parameter space opens up due to the absence of astrophysical X-ray limits [16][17][18][19]. Let us note that, in such a framework, the controversial 3.5 keV line [20][21][22][23][24] does not have a DM origin and the atomic physics explanation is favored [19,[25][26][27].
As in the original νMSM, we successfully identify the parameter space corresponding to non-zero neutrino masses, correct DM relic abundance and the observed amount of baryon asymmetry in the Universe. Hence, we simultaneously address all three of the aforementioned questions. The strongest constraint on our model, arising from the measurements of primordial abundances of light nuclei, forbids the DM candidate to be heavier than ∼ 10 keV.
The structure of the paper is as follows. In section 2 we introduce the model and discuss relevant theoretical and experimental limits. In particular, we show how the neutrino masses are generated and discuss conditions which ensure the correct low-energy neutrino phenomenology (mixing angles and mass squared differences) in our numerical scans. In sections 3 and 4 we identify the viable parameter space for DM and BAU in the model, respectively. Results from these two sections are combined in section 5 with the purpose of finding the regions where DM and BAU can be addressed simultaneously. In section 6 we discuss the implications from structure formation and feasibility of probing this model at various experimental facilities. Finally, in section 7 we conclude.

JHEP08(2018)067 2 The model
We supplement the Standard Model particle content with an additional scalar doublet Σ = (σ + , σ 0 ) T , and three right-handed neutrinos N i . We study the spectrum in which all new scalars are heavier than the right-handed neutrinos, of which the lightest one, N 1 , is a keV-scale DM candidate. We assume that all these newly introduced degrees of freedom have an odd (−) charge under a Z 2 parity symmetry, whereas the SM particles have an opposite, even (+) charge. The scalar sector is therefore equivalent to the one in the inert doublet model [28] with the potential that reads [1,29] where Φ = (φ + , φ 0 ) T is the SM Higgs doublet containing the Higgs boson with mass equal to 2λ 1 v 2 , where v = 246/ √ 2 GeV denotes the vacuum expectation value of φ 0 . By introducing an additional doublet, the scalar sector contains four additional degrees of freedom with masses Here, m ± , m S and m A are the masses of the charged, CP-even and CP-odd scalar, respectively. In the remainder of the paper we will denote charged scalars with σ ± and neutral CP-even (CP-odd) scalar with S (A).
In the fermion sector, the presence of a Z 2 symmetry forbids the "traditional" lepton portal y φNΦ † L + h.c., (2.3) where L is the SM lepton doublet and y φ is a corresponding Yukawa coupling. However, the Yukawa interaction between N i , leptons and Σ field is allowed. After adding a Majorana mass term for right-handed neutrinos, the relevant lepton sector Lagrangian reads where α denotes the SM lepton generations and m N i is the mass of i-th right-handed neutrino. Without loss of generality we take the right-handed neutrino mass matrix in the diagonal form. Note that the replacement Σ → Φ transforms the Yukawa term in eq. (2.4) into the one given in eq. (2.3).

Relevant constraints
Before discussing a realization of nonzero neutrino masses we present the most relevant theoretical and experimental constraints for the considered model.
The most relevant limit arises from the measurement of primordial abundances of light nuclei [30][31][32], known as the big bang nucleosynthesis (BBN). These precisely measured values are sensitive to energy injection into the plasma, for example from the decays of long-lived particles during the BBN epoch. In our model, late decays of N 2 to keV-scale DM (N 1 ) could spoil the BBN predictions. The decay rate for such process is [33] Γ( where l α are the charged SM leptons of generation α. The BBN limit is relevant chiefly due to the tiny y 1α couplings, necessary to keep DM out of the thermal equilibrium with SM bath. The findings from the recent analysis presented in ref. [32] allow us to infer the masses and abundances of N 2 consistent with respect to BBN constraint. A quantitative discussion of the impact of these constraints on our model parameter space is presented in section 5.
Due to the absence of mixing between active and right-handed neutrinos, X-ray limits [14] on keV-scale DM are non-existent within the presented model which opens up a viable parameter space for light right-handed neutrinos. However, one also needs to take into consideration the parameter space excluded by structure formation bounds (Lyman-α forests) and Milky Way satellite counts [34]. Throughout this paper, we will show results consistent with the limits arising from the most conservative scenario where the keV-scale DM is assumed to inherit a thermal distribution function with p /T ≈ 3.1. In section 6 we discuss how this spectrum can be made colder. From refs. [15,[34][35][36][37] we infer that in order to be in accord with the Milky Way satellite count, m N 1 6 keV is viable. The existing limits from Lyman-α forests are stronger but rather controversial due to effects stemming from inter-galactic medium [37,38] and therefore we do not adopt them.
• Lepton Flavor Violation and Scalar Potential.
The scotogenic model predicts lepton flavor violation processes [39] of the type l α → l β γ and l α → 3l β . In table 1 we summarize the bounds given in ref. [40].
The quartic couplings in the scalar potential receive constraints from the requirement of the vacuum stability [29,41] where λ 1 > 0 and λ 2 > 0.
We find the other constraints [42], arising for instance from the compatibility with the electroweak precision data (Peskin-Takeuchi parameters) [43], to be much weaker with respect to the aforementioned ones.

JHEP08(2018)067
LFV process BR upper bound Table 1. Constraints on LFV processes, taken from [40]. The left column indicates rare LFV processes and the right one shows the corresponding upper limits on the branching ratios (BR).

Active neutrino masses
The Scotogenic model is known as one of the most minimal realizations [44][45][46][47][48][49][50][51][52] of nonvanishing neutrino masses established at oscillation experiments [53,54]. In order to obtain the general formula for the neutrino mass matrix in the flavor basis, it is required to calculate the self-energy contributions to the active neutrino propagator, arising from the exchange of both S and A. The Majorana neutrino mass matrix reads [1,39] ( 7) where the summation index i denotes the right-handed neutrino generations. For later convenience, we introduced Λ which abbreviates all the mass matrix components apart from the Yukawa couplings. As we will demonstrate in section 3, in order not to overclose the Universe, the Yukawa couplings of N 1 (y 1α ) must be suppressed with respect to the second and third generation ones. Therefore, N 1 does not effectively yield a contribution to eq. (2.7) which consequently sets the lightest active neutrino to be massless.
For obtaining the Yukawa couplings that are in accord with the active neutrino mass squared differences and mixing parameters, we employ the Casas-Ibarra parametrization [55]. Phenomenologically, following the above discussion on the relevance of only N 2 and N 3 states for active neutrino mass generation, the corresponding 2 × 3 Yukawa submatrix reads Here, Λ 23 = diag(Λ 2 , Λ 3 ) and R is an orthogonal matrix parametrized with one complex angle where ω and ξ are real parameters. The neutrino mass matrix in the mass basis is either for the inverted one. Here, m 2 sol and m 2 atm are solar and atmospheric mass squared differences [56], respectively. Throughout this work we assume a mass spectrum with normal ordering.
The leptonic mixing matrix U PMNS is parametrized with mixing angles θ ij and phases δ i as [56]  Having discussed the construction of the phenomenologically viable Yukawa matrices, we can now estimate the constraints on the parameter space from LFV processes given in table 1. The most stringent bound is coming from µ + → e + γ process with the corresponding one-loop contribution approximately given by (2.14)

Dark matter production
The main goal of this section is to describe the two mechanisms for generating the abundance of keV-scale DM. First, in section 3.1, we discuss DM production via its feeble interactions with thermalized neutral (S, A) and charged scalars (σ ± ). More precisely, the DM abundance gradually increases from the decays of heavy Z 2 -odd scalars, a mechanism known as DM "freeze-in" [57] (see [58,59] for most recent studies). We calculate the magnitude of the coupling required to obtain the observed DM energy density in the Universe. Furthermore, DM is produced via decays of heavier right-handed neutrinos which "freeze-out" [60] from the thermal bath. The viability of such an option is explored in section 3.2. These two production mechanisms do not exclude one another and, in fact, both can significantly contribute in general scenarios.

Production via scalar decays
It is well-known that the standard thermal "freeze-out" of keV-scale DM is not feasible within the considered model. 1 Therefore, we are led toward the "freeze-in" production through which the abundance of non-thermalized N 1 is built up from the decays or annihilations of particles that are in the thermal bath. One may easily infer [33] that, in this model, the dominant production arises from the decays of neutral and charged scalars (A, S, σ ± ) which are in the thermal bath due to the strong gauge interactions. The relevant processes are where ν α correspond to SM neutrinos of generation α. The general form of the Boltzmann equation for the DM production via "freeze-in" reads [57] dn Here, n N 1 is the DM number density, H is the Hubble parameter, Π k abbreviates the phase space factor d 3 p k 2E k (2π) 3 for particle k, |M| 2 Σ→N 1 L is the squared matrix element for the decay of a scalar into N 1 and SM lepton, and f Σ is the distribution function of scalar particles. For brevity, we jointly denote the summation over all scalar decay channels (see eq. (3.1)) with Σ → N 1 L.
The formula given in eq. (3.2) can be simplified to the form [57] dn N 1 dt where m Σ is the mass of a decaying scalar, and K 1 is the modified Bessel function of the second kind. By relating the number density of N 1 with the corresponding yield Y = n N 1 /s, where s is the entropy density, as well as using dT = −H T dt we obtain the following differential equation After changing the integration variable from T to x = m ± /T , as well as expanding eq. (3.4) over all Σ components, we reach the expression where we explicitly inserted formulae for H, s and the partial widths of the scalars [33]. For simplicity, we assumed that y 1α coupling is identical for each flavor α, therefore denoted y 1 . In eq. (3.5), g * is the number of the degrees of freedom in the thermal bath at the time of DM production, M Pl is non-reduced Planck mass, and r A (r S ) denotes the scalar mass ratio m A /m ± (m S /m ± ).

JHEP08(2018)067
10 -2 10 -1 10 0 10 1 10 2 10 -10  After an explicit integration where we assume the reheating temperature to be much higher than any other mass scale in our model, we obtain the following expression for the present DM yield By taking into account all SM and new degrees of freedom, we have g * ≈ 114.25 at O(10 2 − 10 3 ) GeV temperature where the freeze-in occurs.
For the purpose of estimating the required order of magnitude for y 1 , let us further simplify eq. (3.6) by assuming r A = r S = 1 which corresponds to the λ 4,5 → 0 limit. We arrive at By using the well-known relation between DM yield and the relic abundance we finally obtain Ωh 2 FI ≈ 0.12 In figure 1, we present the numerical solutions of eq. (3.5) for different choices of y 1 and fixed m N 1 . We observe that the dominant production occurs at O(1) x. We point out that there is no strong dependence of Y FI on the values of quartic couplings. Hence, any physical choice of these parameters (obeying the constraints from eq. (2.6)) yields a very similar result to the one presented in figure 1. We show the DM freeze-in for m N 1 = 6 keV, JHEP08(2018)067 the value in accord with the structure formation limits given in [15]. For larger values of DM mass, coupling y 1 needs to be smaller according to eq. (3.9).
Let us finally remark that N 2 and N 3 could in principle also be produced via freeze-in from scalar decays. However, if that was the case, from eq. (3.9) we infer that the corresponding Yukawa couplings (y 2α , y 3α ) would then have to be 10 −8 . Such small couplings would not be sufficiently large to generate eV-scale neutrino masses (see eq. (2.7)). The required strength of (y 2α , y 3α ) puts N 2 and N 3 in thermal equilibrium with SM particles.

Production via N 2 decays
As already discussed in section 2, the Z 2 -odd sector consists in total of three right-handed neutrinos, neutral (A and S) and charged scalars σ ± . The mass spectrum is assumed to be . This spectrum has several appealing features. Since scalars are heavier than right-handed neutrinos, potential breaking of parity symmetry induced by renormalization group evolution is evaded [63]. Additionally, models with a low scale DM, in comparison with conventional O(10 2 ) GeV DM, manifest improvement in the predictions for the small scale structure [64].
Due to small Yukawa couplings of N 1 (y 1α y 2α , y 3α ) (see section 3.1), the decay rate for the process N 2 → N 1 L αLβ is much smaller in comparison to the processes involving heavier Z 2 -odd particles. We have numerically checked that such decays occur only after N 2 undergoes a successful thermal freeze-out. Therefore, the amount of DM produced from N 2 decays is determined by the freeze-out abundance of N 2 , which we calculate in the following.
The Boltzmann equation for each of the considered Z 2 -odd particles ( N 2 , N 3 , σ ± , A, S, commonly denoted as χ i , i = 1 . . . N ) is given as [65] where the terms in the first, second and third row on the right hand side represent (co)annihilations, scatterings and decays, respectively. Here, is the cross section for scattering with SM bath (denoted with X and Y ), and Γ ij = X Γ(χ i → χ j X) are decay rates.

JHEP08(2018)067
Since all heavier Z 2 particles eventually decay into N 2 , we may define the total N 2 number density as 2 n = i n i . (3.13) Considering the size of the coupling constants in our setup, all the χ i are in the thermal equilibrium at O(1) TeV temperatures with the number density equal to where g i denotes the number of internal degrees of freedom of species i. By employing eqs. (3.12) and (3.13) and as well as the relation [65] n i n n eq i n eq , (3. 15) we arrive at the standard Boltzmann equation for the evolution of N 2 number density [66] dn dt = −3Hn − σ eff v (n 2 − n 2 eq ), (3.16) where σ eff v is given by [65] Note that in eq. (3.16) there are no scattering terms as well as decays. After employing the summation given in eq. (3.13) these contributions get removed as they do not change the total number of Z 2 -odd particles. Changing variables to dimensionless quantities Y = n/s and y = m N 2 /T yields We have evaluated eq. (3.18) with micrOMEGAs [67] and compared with the output of MadDM [68]. The Yukawa couplings of N 2 and N 3 are generated as described in section 2.2. Therefore, in the code, only the points in the parameter space that are fully consistent with the neutrino mixing angles and mass squared differences are evaluated.
In figure 2, we show the time evolution of Y for several different couplings and righthanded neutrino masses. As can be seen, the yield can reach high values for large λ 5 . This is because, in such case, the Yukawa couplings y 2α are small (see eq. (2.7)) and N 2 does not remain in thermal equilibrium for long (freeze-out takes place at higher temperatures with respect to the usual x ≈ 20).
It is important to point out the coannihilation processes which are dominant if the scalars and N 2 have similar masses. For m ± > m N 2 0.8 m ± , Y is generally lowered by several orders of magnitude with respect to the general case with hierarchical fermion and scalar masses. This can also be seen in figure 2 where the green (purple) curve indicates the case with m N 2 = 0.85 m ± (m N 2 = 0.95 m ± ). We will particularly emphasize the importance JHEP08(2018)067  Figure 2. Evolution of the total N 2 yield for different values of λ 5 in a scenario with vastly different scalar and N 2 masses (red and blue curve) as well as in the case when the masses are similar (green and purple dashed curve). In the latter case, due to larger interaction strength, N 2 remains in thermal equilibrium longer, yielding orders of magnitude smaller N 2 abundance with respect to the former case.
of coannihilations when exploring viable joint parameter space for DM and leptogenesis (see section 5). The contribution to the DM relic abundance from N 2 decays is where the relation between Ωh 2 N 2 and Y matches the one in eq. (3.8). Taking into account the complementary contribution from scalar decays given in section 3.1, the requirement for having the amount of DM in accordance with the measurements [69] is As already discussed, we construct the Yukawa couplings of N 2 and N 3 in agreement with the results from neutrino oscillation experiments. Such a choice already unambiguously fixes Ωh 2 N 2 →N 1 contribution, meaning that if it is already overshooting the observed value of 0.12, the corresponding parameter choice is excluded. Otherwise, the Yukawa couplings of N 1 could be accommodated in such a way to satisfy eq. (3.20). Our findings in section 3.1 indicate y 1α 10 −8 .
In figure 3, results for Ωh 2 N 2 →N 1 are shown as a function of m N 2 . For m N 2 500 GeV one may infer that in the case of larger λ 5 couplings, Ωh 2 N 2 →N 1 poses the dominant contribution to DM relic density since, due to the weakness of N 2 N 2 → SM SM channels, N 2 undergoes freeze-out very early. On the other hand, for higher m N 2 , additional annihilation channels involving new scalars become dominant, increasing the overall cross section, and correspondingly reducing Ωh 2 N 2 →N 1 . In contrast, for very tiny λ 5 (pink curve), N 2 annihilates efficiently even in the absence of coannihilation processes, and therefore the Log

Leptogenesis
In order to generate the observed BAU we study the option where initially an asymmetry in the lepton sector is produced [70]. The production of lepton asymmetry is dubbed leptogenesis. In the simplest realizations, the models featuring leptogenesis incorporate the seesaw type-I mechanism and rely on CP violating decays of hierarchical heavy right-handed neutrinos. The produced lepton asymmetry can be partially converted to the baryon asymmetry due to the existence of non-perturbative processes in thermal equilibrium [71]. These processes, called sphalerons, violate B+L, but conserve B−L numbers and thus allow for an asymmetry conversion between the two sectors. When the temperature drops below T s = 131.7 GeV [72], the sphalerons decouple from the thermal bath and the asymmetry conversion ceases.
The drawback of the above mentioned realization is that in order to generate the observed BAU the masses of the right-handed neutrinos have to be at least of O(10 8 ) GeV [73,74] which is not experimentally reachable. 3 It is, however, possible to lower the needed mass scale below TeV. In this paper we present such an option in the JHEP08(2018)067 frame of the Scotogenic model. We rely on two competing mechanisms to create a sufficient baryon asymmetry.
Firstly, we incorporate oscillations among the right-handed neutrinos (ARS) [2,[80][81][82][83] which transfer the asymmetry to the SM leptons via lepton portal (see eq. (2.4)). Secondly, we also take into account decays of the new Σ doublet which at finite temperatures serves as an additional source of CP violation. This process is suppressed by a mass insertion factor (m N 2 /T ) 2 and was neglected in the past literature. However, recently it has been shown that scalar decays are important [3] and can even dominate ARS in some regions of the parameter space [84]. Hence, we take both ARS and scalar decays into account.
It is worth noting that several attempts to implement leptogenesis in the Scotogenic model had already been made (see e.g. [85,86]). In refs. [87,88], the authors were able to achieve a low scale leptogenesis without imposing mass degeneracies between right-handed neutrinos.
For studying leptogenesis, we apply the procedure presented for seesaw type-I mechanism in ref. [84] to the Scotogenic model. In other words, we replace SM Higgs doublet with Σ in the lepton portal term. Only the two right-handed neutrinos, namely N 2 and N 3 participate in the lepton asymmetry production. In the density matrix formalism, the production of a lepton asymmetry via both ARS mechanism and scalar decays is automatically included. Adopting the notation from [84], we work with yields ρ associated with the number density matrices of right-handed neutrinos Here, n is a 2 × 2 matrix with diagonal elements associated to the densities and the offdiagonal entries parametrizing the mixing between N 2 and N 3 fields. The states with the opposite helicity are indicated with the overline. In order to calculate the lepton asymmetry δY l , for each flavor l, we solve the following set of coupled differential equations

JHEP08(2018)067
where the parameter z equals z = T s /T and T s ∼ 131.7 GeV is the temperature at which sphalerons decouple. Integration is performed between z start = T sph /T = 10 −3 and z stop = 1 with the initial conditions ρ αβ (z start ) = 0,ρ αβ (z start ) = 0, δY l (z start ) = 0, where indices α and β run from 1 to 2 due to the relevance of only two right-handed states. In order to obtain the final baryon asymmetry δY B we evaluate − 2 3 l=e,µ,τ δY l [89]. The most general form for the diagonal matrix elements of E N is where g = 1 for right-handed neutrinos. In the relativistic case, the expression in eq. (4.5) reduces to [84] where we introduced the level of degeneracy between the masses of N 2 and N 3 defined via relation m N 3 = m N 2 (1 + δ M ). Note that, as we will show later explicitly, δ M 10 −8 for a successful leptogenesis [84].
In the numerical evaluation, we have employed the general expressions for all the terms involving right-handed neutrinos given in eqs. (4.2), (4.3) and (4.4), i.e. we did not take the relativistic approximation. This is because our findings indicate that successful leptogenesis occurs when m N 2,3 T s , for which the relativistic approximation does not hold. In contrast, the equilibrium number density for SM leptons is evaluated in the relativistic approximation since the heaviest SM lepton is two orders of magnitude lighter with respect to T s . In eqs. (4.2), (4.3) and (4.4) we denote the reaction densities for ARS and Σ decays with γ LC and γ LV , respectively, where LC and LV are abbreviations for the lepton number conserving and violating processes. Washout terms are labeled with extra subscripts -W C or W Q where the former (latter) indicates lepton asymmetry loss due to classical (quantum) effects. The derivation of the reaction densities and the washout terms is outlined in appendix A. Here, we would like to emphasize the scalar mass treatment in such calculation. First of all, in evaluating the integrals over the phase space, for simplicity and the possibility of analytical evaluation we take m ± = m S = m A . We have nevertheless checked numerically that introducing the splitting between these masses, arising from λ 4,5 = 0 (see eq. (2.2)), yields an insignificant change of reaction densities and washout terms which does not influence the final value of baryon asymmetry.

JHEP08(2018)067
In addition to the bare masses for new scalars, which are O(10 2 −10 3 ) GeV, the thermal corrections are important. We therefore solve eqs. (4.2), (4.3) and (4.4) in two different temperature regimes. For the low temperatures, we neglect thermal corrections and use while for high temperatures we adopted the expression for thermal corrections, calculated in refs. [90] and [91] where g and g are SU(2) L and U(1) Y coupling constants, respectively. In total, we are using the following prescription for the scalar mass otherwise.
(4.11) The boundary temperature separating the two regimes is determined by the condition that the thermal mass is equal to 3 times the bare mass. We checked numerically that varying the boundary temperature by O(1) numbers does not affect the final results. Below this boundary we drop thermal corrections to the leptons, considering them to be effectively massless at low temperatures. Let us note that in both regimes, the phase space for Σ decays is kinematically open, i.e. scalar masses are always larger than the total mass of decay products (right-handed neutrino and SM lepton). We would also like to stress that we have checked that in the relativistic limit our results match those presented in ref. [84].
The temperature dependence of the γ terms (reaction densities and washout terms) is shown in figure 4 where we compare the relativistic regime (neglecting all bare masses and m N , solid lines), with the general aforementioned approach (taking the scalar mass as given in eq. (4.11) and accounting for the effects from non-zero right-handed neutrino masses, dashed lines). We note the suppression of γ LV terms at low temperatures. This is because, in the absence of thermal effects, it is much less probable to have an on-shell mediator particle [3]. The γ LC factors feature an opposite effect -get larger at z ∼ 0.1 where the phase space for scattering processes between L and Σ increases. This is because m ± is fixed to the bare mass and the lepton mass is gradually decreasing due to the dropping temperature. We also observe that all γ terms in both panels become Boltzmann suppressed at z ∼ 1.
Radiative generation of neutrino masses introduces a suppression factor of λ 2 5 /16π 2 which lifts up the overall strength of (y 2α , y 3α ) in comparison to the corresponding values in seesaw type-I model. Numerically, these couplings can not be made smaller than O(10 −6 ) in our model. Such interaction strength may pose a problem for generating observed BAU due to strong washout of the generated lepton asymmetry at late times. Let us note that, in principle, leptogenesis can be successfully achieved with such large couplings [84] when ξ (see eq. (2.9)) is large. In addition, as pointed out in ref. [83], the large imaginary JHEP08(2018)067 entries (induced by ξ) may lead to certain cancellations between different terms in the Boltzmann equations, preventing the generated lepton asymmetry from strong washout effects. However, this reasoning can not be applied to our model because pushing ξ to large values drastically increases the strength of washout effects. In fact, even the asymmetry produced at very late times will be washed out efficiently before sphaleron decoupling.
In order to reduce the Yukawa couplings as much as possible we choose the quartic couplings as given in table 2, consistent with the limits from the stability of the Higgs potential (see section 2.1). With couplings chosen in such a way we find that successful leptogenesis may only be generated with m N 2 O(10 2 ) GeV. For the right-handed neutrino masses smaller than T s , the washout effects at z 0.1 remove the vast majority of the generated asymmetry. In contrast, if m N 2 T s , the right-handed neutrino abundance becomes Boltzmann suppressed at z 10 −1 which strongly suppresses washout integrals given in appendix A. While at such late times the asymmetry production is suppressed in a similar way, we find that the strong production at z 1 suffices for generating the observed δY B = 0.86 × 10 −10 [69] baryon asymmetry. Throughout the paper, without the loss of generality, we use only positive values of ξ. This parameter enters exponentially (exp ±ξ ) in the Yukawa matrix and thus its negative values would not change qualitatively the overall picture for BAU generation.  Figure 5. Evolution of ρ 11 (z), ρ 12 (z) and δY B (z) for two different values of δ M . We have fixed ω = π/4, δ = −π/2, ξ = 2 and α 1 = α 2 = 0. The asymmetry for δ M = 10 −8 is produced earlier but also more strongly washed out with respect to δ M = 10 −10 case. We found no significant dependence of the generated asymmetry on the values of Dirac (δ) and Majorana (α 1 and α 2 ) CP phases.
In figure 5 we show the deviation from the equilibrium value of N 2 number density (ρ 11 (z)/ρ N eq − 1) and the mixing between N 2 and N 3 (ρ 12 (z)). The evolution of the baryon asymmetry δY B is also presented for m N 2 = 200 GeV, µ 2 = 800 GeV and with the following values of angles and phases ω = π/4, δ = −π/2, α 1 = α 2 = 0. The kink at z ≈ 0.05 is due to the change of regimes for the thermal masses (see eq. (4.11)). The narrow feature at slightly smaller temperatures indicates the sign change of δY B . A significant deviation of ρ 11 and ρ 22 from ρ N eq as well as avoiding very tiny off-diagonal matrix elements (ρ 12 and ρ 21 ) are crucial for the successful leptogenesis. This is only ensured for m N 2 T s , as discussed above. We would like to stress that there is no significant dependence of the generated lepton asymmetry on α 1 and α 2 . In addition, note that δ M has to be tiny for generating the asymmetry, implying strong level of degeneracy between heavier right-handed neutrino masses. From the shape of ρ 12 (z) and δY B (z) curves for δ M = 10 −8 (left panel) and δ M = 10 −10 (right panel) we infer that in the former case, the washout effects are effective throughout a longer time period. This can be understood from the temperature scale z osc at which the oscillation among right-handed neutrinos is most effective. It is given by [2] z osc = T s 2 45/(4π 3 (4.12) For smaller δ M , the asymmetry is produced later which leads to a weaker washout and consequently larger δY B . In figures 6a and 6b we show the generated baryon asymmetry in ξ − M N 2 plane for two fixed values of λ 5 . From this figure one can also qualitatively infer the dependence of the generated baryon asymmetry on the Yukawa interaction strength since different values of λ 5 imply different y 2,3 α . As expected, there is only a mild dependence on m N 2 which does not alter the Yukawa strength dramatically. From figure 6c one can deduce that there is a saturation effect for δY B 4 × 10 −9 , where the dependence of the asymmetry on λ 5 flattens. Furthermore, one should note that in case of large λ 5 , δY B is generally very weakly dependent on this quartic coupling. Finally, figure 6d summarizes the interplay between λ 5 JHEP08(2018)067 (a) δYB in (ξ, mN 2 ) plane for λ5 = 2.
(c) δYB in (log λ5, mN 2 ) plane, ξ = 1.  and ξ. Intuitively, it is clear that the largest asymmetry production happens for low ξ and high λ 5 values, whereas going to the opposite regime favors washout, hence leading to a reduced final asymmetry. Interestingly, ξ has a stronger impact than the quartic coupling, despite the fact that in the considered parameter range both quantities vary the size of the Yukawa coupling by almost two orders of magnitude. The reason for this behavior is that while λ 5 just sets the overall strength of the Yukawa couplings, ξ enters in the Yukawa matrix in a specific pattern, leading to rather non trivial effects.
From figure 7 we observe that when considering only ARS, |δY B | decreases roughly by one order of magnitude with respect to the values from figure 6a. What can also be inferred is the stronger dependence on m N 2 in absence of scalar decays. We find that it is very important to account for the asymmetry production via new scalar decays in addition to ARS. Actually, turning either of the two production mechanisms off yields a loss of a substantial amount of the generated baryon asymmetry, as evident from figure 8. Obviously, there is a strong interplay between both regimes, because their individual contributions are rather small and do not trivially add up to yield the final result. Each of the regimes individually leads to a similar final asymmetry contribution. Initially, the asymmetry generation is governed exclusively by the ARS mechanism and Σ decays start to compete at z ≈ 0.005. Interestingly, the intermediate asymmetry produced by Σ decays is one order of magnitude larger compared to the one in ARS, but tends to get washed out more strongly. At z ≈ 0.05, δY B changes its sign and later fades away due to washout effects. At even smaller temperatures all regimes feature similar asymmetry evolution (governed by washout effects) until sphaleron processes eventually decouple. Again, including both regimes does not translate into summing up their individual contributions. The asymmetry generation in the presence of both mechanisms shows a qualitatively different behavior, featuring a peak at larger temperatures. Non-linear combination of ARS and scalar decays explicitly demonstrates the absence of the weak-washout regime in which the contributions would add trivially.

Identifying the viable joint parameter space for DM and leptogenesis
In sections 3 and 4 we separately scrutinized DM and BAU production. Therefore, the question arises whether it is possible to find regions in the parameter space for which the observed values of DM relic abundance and |δY B | are simultaneously reached. Generally, these two mechanisms have conflicting requirements on the strength of the Yukawa coupling. In shown together with the general case when the both production mechanisms are included (yellow solid). We take ξ = 3 and m N2 = 400 GeV. It is evident that the combination of both production mechanisms suffices to account for BAU, in contrast with the individual contributions.
order not to overproduce DM from decays of N 2 (section 3.2), sufficiently large Yukawa interactions are required. On the other hand, leptogenesis relies on weak interactions as otherwise the washout effects would easily destroy any generated lepton asymmetry.
A key method for curing this contrast is to impose coannihilations between the righthanded neutrinos and the scalars (especially the lightest scalar σ ± ) by choosing them to have similar masses. The coannihilation case was already discussed in section 3 (see in particular figure 2). Such a regime opens up new scalar annihilation channels which do not rely on the strength of the second and third generation Yukawa couplings. Therefore, a huge suppression of the relic density can be achieved for Yukawa couplings set low enough to generate a significant lepton asymmetry.
We now revisit BBN constraints on the N 2 decays introduced in section 2.1. By requiring N 2 to decay before t BBN ∼ 1 sec we obtain |y 2α | 2 6.3 · 10 −7 m ± 1 TeV After a more careful analysis of the processes influencing the primordial abundances of light elements this limit gets significantly relaxed. Following ref. [32], we infer that for Y m N 2 O(10 −9 ), where Y is the DM yield (see eq. (3.18)), the BBN limit given in eq. (5.1) is relaxed by roughly three orders of magnitude. This is because, in our model, N 2 has only leptonic decays and the strongest effect from charged leptons on the primordial abundances of the light nuclei is coming from the photodissociation process which is most effective at approximately 1000 seconds after the Big Bang.
We calculated the relic density Ωh 2 N 2 for different choices of µ 2 and also evaluated δY B in a range of m N 2 and ξ, taking into account m N 2 < m ± . Then, we determined which region in the parameter space is ruled out due to low y 2α . We conservatively adopted the JHEP08(2018)067 . Allowed region in the (ξ, m N2 ) parameter space for fixed µ 2 = 600 GeV (left) and µ 2 = 800 GeV (right). The BBN exclusion limits are shown in red, while the blue shaded region does not produce a sufficiently large baryon asymmetry. We observe that there is a substantial region consistent with BBN limits in which the correct amounts of DM and baryon asymmetry can be obtained.
In figure 9 we present the results for two different choices of µ 2 . In both panels, the blue region is indicating excluded parameter space due to the insufficient amount of generated asymmetry and the red region is excluded by BBN.
There are several things we would like to point out. First, the final baryon asymmetry only mildly depends on the involved particle masses and in general larger mass scales will only lead to a slight decrease of δY B as can be seen from the slope of the BAU line in figure 9. In contrast, choosing a higher mass scale weakens the BBN bounds, e.g. for a DM mass of 6 keV we need Ωh 2 N 2 8.5 for µ 2 = 600 GeV, while for µ 2 = 800 GeV the limit is relaxed to Ωh 2 N 2 88.7. These values also suggest that the N 2 decay contribution (Ωh 2 N 2 →N 1 ) to the DM relic abundance is negligible and thus y 1 is fixed by the requirement that DM is completely produced via N 1 freeze-in. Larger scalar masses lead to stronger Yukawa interactions in order to maintain the correct relic density (see eq. (3.9)).
In figure 9, the blue region indicates the region excluded for δ M ≤ 10 −11 . A smaller level of degeneracy between right-handed neutrinos, such as δ M = 10 −10 , would correspond to the shift of this exclusion to the left (toward BBN bound). For δ M ≥ 10 −10 and m N 1 ∼ 6 keV we find no parameter space free from either BBN or BAU exclusions.
We have also observed the dependence of BBN limits on the maximum allowed DM mass. For µ 2 in the range [300, 1000] GeV, we found a maximal value of the DM mass of 9.4 keV, which is consistent with structure formation bounds. Generally, choosing higher degeneracies will open up the available parameter space, thus allowing larger DM masses. However, even in such cases we estimated the maximal allowed DM mass to be at most O(10) keV.

JHEP08(2018)067
In summary, we identified the parameter space in which the produced DM abundance and BAU are in accord with the observed values. We have seen that the most stringent constraints arise from BBN considerations, which can be relaxed by employing coannihilation processes between N 2,3 and Σ which effectively put an upper bound on the allowed mass for the DM. We found that the DM production in our model is mainly driven by the freeze-in.
In this section we have shown that the three biggest challenges of particle physics, pointed out in section 1, can be successfully and simultaneously solved within the considered model.

Detection prospects
In section 2.1 we discussed the limits from structure formation on keV-scale DM. Here, we wish to point out that these limits are relaxed in this model which opens up some parameter space for such DM candidate. Let us illustrate how DM can be made colder, i.e. less constrained from structure formation considerations. We discuss the most relevant case where the observable DM relic density is dominantly set through the freeze-in from the decays of heavy scalar particles (see section 3.1), as the strong production from N 2 decays is in tension with BBN limits. The production of DM occurs at the mass scale of decaying charged and neutral scalars. Due to the production at such high temperatures, cooling of DM particles is efficient. Namely, the effective temperature of the DM sector, when compared to photons, is reduced by the amount of entropy dilution factor which is ≈ 2.9 [24,92].
Having the absence of X-ray signal, we reach the conclusion that the testability of our model is currently limited only to the searches at the LHC as well as the facilities probing lepton flavor violation processes. The limits coming from the latter are discussed in section 2.1 and consistently taken into account throughout the paper. In this section, therefore, we mainly comment on the LHC prospects, which were studied for this model in [93] where the dominant production of DM is assumed to be via freeze-in through scalar decays. Among others, the authors are considering the regime where scalar particles are heavier than all three generations of right-handed neutrinos, which matches our setup.
The answer to the question which LHC search has the strongest sensitivity depends on the mass difference between m N 2 and m N 3 . In section 4, we showed that baryogenesis can be achieved via generation of a lepton asymmetry, where the crucial ingredient was the approximate degeneracy between N 2 and N 3 states. In that case, the LHC searches are limited to σ ± → N 2,3 l ± α , i.e. channels including charged leptons in the final state. Charged scalars σ ± are produced in pairs either via gluon fusion or Drell-Yan processes [93], and therefore the expected signature consists of two prompt charged leptons and missing energy due to elusive right-handed neutrinos which are decaying to N 1 only after leaving the detector.
Let us note that with an assumption of a complementary baryogenesis mechanism that goes beyond our framework (for instance decays of very heavy singlets [94] or electroweak baryogenesis [95]) for which tiny δ M is not required, the discovery potential at LHC increases. Namely, if we assume a mass gap between N 2 and N 3 to be 10 GeV there is another viable search in this model -displaced vertices [96]. Even though decays of heavier right-handed neutrinos into DM (N 1 ) do not happen within the detector, the decay JHEP08(2018)067 N 3 → N 2 l ± α l ∓ β can be rapid enough. Then, the displaced leptons from the initial scalar decay and the subsequent N 3 decay may be observed. Limits for both leptons+ / E T and displaced vertices are presented in [93]. Let us emphasize that displaced vertex search can constrain y 3α and y 2α couplings up to the level of 10 −4 , which is two orders of magnitude stronger with respect to the limits from LFV assuming no strong hierarchy between the entries of the Yukawa matrix. The existence of the upper limit is a consequence of the requirement for the resolution of displaced vertices. It is worth mentioning that the limits from LHC searches attenuate very quickly above EW scale. For more details we refer the reader to ref. [93].
Let us also remark that the presence of extra charged scalars implies the radiative contribution to the decays of SM Higgs particle into a pair of photons or a photon and a Z boson [97,98]. Specifically, for the diphoton channel, in order to get constructive (destructive) interference between SM and new physics contribution, one requires negative (positive) Higgs portal couplings. Significant deviation from measured values [40] are achieved for rather large values of such couplings [99,100] which are avoided in our analysis.
We finalize this section with a comment on the detection prospects for N 2,3 discovery at future lepton colliders as well as SHIP [101], that were proposed recently [80,102,103]. At lepton colliders, right-handed neutrinos are assumed to be generated in the process e + e − → N 2,3 ν α [104], for which the mixing between active and right-handed states is required. The mixing is necessary also for SHIP where right-handed neutrinos would get produced in the decays of heavy hadrons. Hence, due to the absence of the mixing induced by an exact Z 2 parity symmetry, right-handed neutrinos in this model are not testable at these facilities.

Summary and conclusions
In this work we studied the generation of neutrino masses, dark matter and baryon asymmetry of the Universe within the so called Scotogenic model, where three right-handed neutrinos and an additional scalar doublet, all odd under a Z 2 parity symmetry, are added to the SM. Active neutrino masses are obtained radiatively via loops involving new particles. We considered a mass spectrum where all scalar masses are at or below the TeV-scale. Furthermore, we invoked a hierarchy in the right-handed neutrino mass spectrum, choosing m N 1 ≈ O(1) keV and m N 2,3 ≈ O(100) GeV, where the lightest state, N 1 , is a keV-scale DM particle.
For DM production, we looked at two complementary contributions: first, we examined freeze-in of N 1 from the decays of new scalars. Second, we also took three-body decays of "frozen-out" N 2 into account. We were able to derive the correct DM density in a wide parameter region.
As the baryon asymmetry is concerned, we studied leptogenesis from the combination of right-handed neutrino oscillations and scalar decays. We showed that it is possible to derive a significant asymmetry in case of highly degenerated right-handed neutrino masses. We have established that it is crucial to set the mass of heavier two states to O(10 2 ) GeV in order to prevent late time washout.
Finally, in finding the joint parameter space for DM and BAU, we have shown that the BBN bound plays a major role and rules out a large portion of the parameter space, effectively forbidding DM mass to exceed ∼ 10 keV. Nevertheless, by imposing coannihilla-JHEP08(2018)067 tions between fermions and scalars, we were successful in identifying the regions where both DM and BAU are produced in right amounts. Hence, in the considered model, we have simultaneously solved the three greatest mysteries in particle physics.

Acknowledgments
We thank Joachim Kopp for the collaboration in the initial stage of this project. We are greatly indebted to Daniele Teresi for cross checking our leptogenesis implementation in the initial phase of this work as well as helping us understand several key leptogenesis features, in particular the production from scalar decays. VB would also like to thank Mikhail Shaposhnikov and Kai Schmitz for very insightful and useful discussions on the presented model. Work in Mainz is supported by the Cluster of Excellence Precision Physics, Fundamental Interactions and Structure of Matter (PRISMA -EXC 1098).

A Phase space integration and reaction densities
Here we present the results for the reaction densities and washout terms, obtained by performing the phase space integration. In the following, we make the abbreviation E i ≡ E i (p i ) for SM lepton (E L ), right-handed neutrino (E N ) and scalar (E Σ ) energies. The right-handed neutrino and SM lepton momenta are denoted with k and p, respectively. We also use the approximation m N 2 m N 3 ≡ m N . The phase space integral is evaluated as where θ 12 is the angle between the outgoing lepton and right-handed neutrino and f (cos θ 12 ) ≡ p 2 + k 2 + 2pk cos θ 12 + m 2 ± − k 2 + m 2 N − p 2 + M 2 L . (A.2)

JHEP08(2018)067
After a further evaluation we reach the expression with All reaction densities and washout terms include After inserting appropriate distribution functions we finally obtain JHEP08(2018)067 Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.