Dark Matter in the minimal Inverse Seesaw mechanism

We consider the possibility of simultaneously addressing the dark matter problem and neutrino mass generation in the minimal inverse seesaw realisation. The Standard Model is extended by two right-handed neutrinos and three sterile fermionic states, leading to three light active neutrino mass eigenstates, two pairs of (heavy) pseudo-Dirac mass eigenstates and one (mostly) sterile state with mass around the keV, possibly providing a dark matter candidate, and accounting for the recently observed and still unidentified monochromatic 3.5 keV line in galaxy cluster spectra. The conventional production mechanism through oscillation from active neutrinos can account only for $\sim 43\%$ of the observed relic density. This can be slightly increased to $\sim 48\%$ when including effects of entropy injection from the decay of light (with mass below 20 GeV) pseudo-Dirac neutrinos. The correct relic density can be achieved through freeze-in from the decay of heavy (above the Higgs mass) pseudo-Dirac neutrinos. This production is only effective for a limited range of masses, such that the decay occurs not too far from the electroweak phase transition. We thus propose a simple extension of the inverse seesaw framework, with an extra scalar singlet coupling to both the Higgs and the sterile neutrinos, which allows to achieve the correct dark matter abundance in a broader region of the parameter space, in particular in the low mass region for the pseudo-Dirac neutrinos.


Introduction
The origin of neutrino masses and the nature of dark matter are two of the most pressing open questions of particle and astroparticle physics.
The basic structure of a 3-flavour leptonic mixing matrix, as well as the existence of two neutrino oscillation frequencies (∆m 2 ij ) -in turn suggesting that at least two neutrinos have non-vanishing masses -, are strongly supported by current neutrino data (for a review, see [1]). However, some of the current experiments (reactor [2], accelerator [3,4] and Gallium [5]) point to the existence of extra fermionic gauge singlets with masses in the eV range. This would imply that instead of the three-neutrino mixing scheme, one would have a 3 + 1-neutrino (or 3+more) mixing schemes (see, for instance, [6]).
In particular, sterile neutrinos with masses around the keV can be viable Warm Dark Matter (WDM) candidates. They can potentially solve, even if providing only a fraction of the total dark matter (DM) relic density some tensions with structure formation observations [11]. In addition, a sterile neutrino at this mass scale could in general decay into an ordinary neutrino and a photon which could be detected in cosmic rays. This last possibility has recently triggered a great interest in view of the indication, yet to be confirmed, of an unidentified photon line in galaxy cluster spectra at an energy ∼ 3.5 keV [12,13].
In order to account for massive neutrinos, one of the most economical possibilities is to embed a seesaw mechanism [14][15][16] into the Standard Model (SM). Low-scale seesaw mechanisms [17][18][19][20][21][22], in which sterile fermions with masses around the electroweak scale or even lower are added to the SM particle content, are very attractive scenarios. In particular, this is the case of the Inverse Seesaw mechanism (ISS) [17][18][19], which offers the possibility to accommodate the smallness of the active neutrino masses m ν for a comparatively low seesaw scale, but still with natural Yukawa couplings. The new (light) states can be produced in collider and/or low-energy experiments, and their contribution to physical processes can be sizeable in certain realisations. The ISS mass generation mechanism requires the simultaneous addition of both right-handed (RH) neutrinos and extra sterile fermions to the SM field content (#ν R = 0 and #s = 0) 1 . Being gauge singlets, and since there is no direct evidence for their existence, the number of additional fermionic singlets #ν R and #s is unknown (moreover the numbers #ν R and #s are not necessarily equal [23]). However, these scenarios are severely constrained: any (inverse) seesaw realisation must comply with a number of bounds. In addition to accommodating neutrino data (masses and mixings), they must fulfill unitarity bounds, laboratory bounds, electroweak (EW) precision tests, LHC constraints, bounds from rare decays, as well as cosmological constraints.
In [23], it was shown that it is possible to construct several minimal distinct ISS scenarios that can reproduce the correct neutrino mass spectrum while fulfilling all phenomenological constraints. Based on a perturbative approach, this study also showed that the mass spectrum of these minimal ISS realisations is characterised by either 2 or 3 different mass scales, corresponding to the one of the light active neutrinos m ν , that corresponding to the heavy states M R , and an intermediate scale ∼ µ only relevant when #s > #ν R . This allows to identify two truly minimal ISS realisations (at tree level): the first one, denoted "(2,2) ISS" model, corresponds to the SM extended by two RH neutrinos and two sterile states. It leads to a 3-flavour mixing scheme and requires only two scales (the light neutrino masses, m ν and the RH neutrino masses, M R ). Although considerably fine tuned, this ISS configuration still complies with all phenomenological constraints, and systematically leads to a normal hierarchy for the light neutrinos. The second, the "(2,3) ISS" realisation, corresponds to an extension of the SM by two RH neutrinos and three sterile states, and allows to accommodate both hierarchies for the light neutrino spectrum (with the inverse hierarchy only marginally allowed), in a 3+1-mixing scheme. The mass of the lightest sterile neutrino can vary over a large interval: depending on its regime, the "(2,3) ISS" realisation can offer an explanation for the short baseline (reactor/accelerator) anomaly [2][3][4][5] (for a mass of the lightest sterile state around the eV), or provide a WDM candidate (for a mass of the lightest sterile state in the keV range).
In this work, we investigate in detail this last possibility, conducting a thorough analysis of the relic abundance of the warm dark matter candidate, taking into account all available phenomenological, astrophysical and cosmological constraints. The conventional DM production mechanism, the so called Dodelson-Widrow mechanism [7], results in a tension with observational constraints from DM Indirect Detection (ID) and structure formation, since it can only account for at most ∼ 50% of the total DM abundance. A sizeable DM density can nonetheless be achieved when one considers the decay of the heavy pseudo-Dirac neutrinos. However this possibility is realised only in a restricted region of the parameter space. An extension of the model is thus needed in order to account for a viable DM in a broader portion of the parameter space.
The paper is organised as follows: in Section 2, after the description of the model -the "(2,3) ISS" realisation -, we address the prospects of the lightest sterile state as a viable DM candidate, which are stability, indirect detection and the dark matter generation mechanism. In Section 3, we consider all the relevant different astrophysical and cosmological constraints taking into account the effect of the heaviest sterile neutrinos (DM production from decays of heavy sterile states or possible entropy injection effects from a scenario with lighter sterile neutrinos) accounting as well for the indication of the monochromatic 3.5 keV observed line. Section 4 is devoted to an economical extension of the present model which succeeds in providing the remaining ∼ 50% of the dark matter relic abundance in a larger region of the parameter space. Our final remarks are given in Section 5. The numerical details regarding the production and evolution of the sterile neutrinos can be found in the Appendix.
2 Description of the model

The "(2,3) ISS" framework
The inverse seesaw mechanism can be embedded into the framework of the SM by introducing a mass term for neutrinos of the form: where C ≡ iγ 2 γ 0 is the charge conjugation matrix and n L ≡ ν L,α , ν c R,i , s j T . Here ν L,α , α = e, µ, τ , denotes the three SM left-handed neutrino states, while ν c R,i (i = 1, #ν R ) and s j (j = 1, #s) are additional right-handed neutrinos and fermionic sterile states, respectively. Since there is no direct evidence for their existence, and being gauge singlets, the number of the additional fermionic sterile states ν R,i and s j is unknown. In this work we will follow a bottom-up approach, focussing on the ISS realisations with the minimal content of extra right-handed neutrinos and sterile fermionic states providing a viable phenomenology. It is nonetheless worth mentioning that a realisation of the ISS with 3 RH neutrinos and 4 sterile states fulfilling all possible constrains has been recently found in the context of conformal EW symmetry breaking [24].
The neutrino mass matrix M has the form: where d, m, n, µ are complex matrices. The Dirac mass matrix d arises from the Yukawa couplings to the SM Higgs boson H = iσ 2 H, which gives after Electroweak Symmetry Breaking (EWSB): The matrices m and µ represent the Majorana mass terms for, respectively, right-handed and sterile fermions. Assigning a leptonic charge to both ν R and s with lepton number L = +1 [17][18][19], in order that the Dirac mass term −d * ν L ν R is lepton number conserving, the terms ν T R Cν R and s T Cs violate lepton number by two units. Within this definition, the entries of m and µ can be made small, as required for accommodating O(eV) masses for ordinary neutrinos through the ISS mechanism. This does not conflict with naturalness since the lepton number is restored in the limit m, µ → 0. In the following we will always work under the assumption that the magnitude of the physical parameters (entries of the above matrices) fulfill such a naturalness criterium |m|, |µ| |d|, |n| , which implies the condition #s ≥ #ν R for the mechanism to be viable [23]. Finally, the matrix n represents the lepton number conserving interactions between righthanded and new sterile fermions. Notice from Eq. (2) that we do not consider lepton flavour violating Yukawa terms.
The physical neutrino states ν I, I=1,...,3+#ν R +#s , are obtained upon diagonalization of the mass matrix M via the unitary leptonic mixing matrix U , and feature the following mass pattern [23]: • 3 light active states with masses of the form This set must contain at least three different masses, in agreement with the two oscillation mass frequencies (the solar and the atmospheric ones).
In order to be phenomenologically viable, the matrix M associated to a given ISS realisation must exhibit, upon diagonalization, three light, i.e.
O(eV), active eigenstates with mass differences in agreement with oscillation data and a mixing pattern compatible with the experimental determination of the PMNS matrix. Although the minimal ISS realisation satisfying these requirements is the "(2,2) ISS", corresponding to 2 right-handed and 2 additional sterile fermions, a detailed numerical study performed in [23] has shown that complying with all constraints requires an important fine-tuning.
On the contrary, a very good "fit" is provided by the "(2,3) ISS" realisation (2 right-handed and 3 additional sterile fermions). In this last case an additional intermediate state (with mass m 4 = m s = O(µ)) appears in the mass spectrum. Remarkably, and in order to comply with all constraints from neutrino oscillation and laboratory experiments, the coupling of this new state to the active neutrinos must be highly suppressed, thus leading to a dominantly sterile state, with a mass ranging from O(eV) to several tens of keV. As consequence of its very weak interactions, the lifetime of the lightest sterile neutrino largely exceeds the lifetime of the Universe and it can thus play a relevant rôle in cosmology.
In this work we will focus on the possibility that this sterile neutrino accounts, at least partially, for the Dark Matter component of the Universe, identifying the viable regions of the parameter space with respect to DM phenomenology of the "(2,3) ISS" model.
We nonetheless mention that, in the lightest mass region (eV scale), and although not a viable DM candidate, the lightest sterile neutrino could potentially accommodate a 3 + 1 mixing scheme explaining the (anti)-neutrino anomalies in short baseline, Gallium and reactor experiments [2][3][4][5].

Light sterile neutrino as Dark Matter
Before the analysis we will briefly summarize the main issues that should be addressed in order for the lightest sterile neutrino to be a viable dark matter candidate. Stability and Indirect Detection: The most basic requirement for a DM candidate is its stability (at least on cosmological scales). All the extra neutrinos of the ISS model have a non zero mixing with ordinary matter. As a consequence, the lightest one is not totally stable and can decay into an active neutrino and a photon γ. On the other hand, as already pointed out, its very small mixing makes the decay rate negligible with respect to cosmological scales. Nonetheless, a residual population of particles can decay at present times producing the characteristic signature of a monochromatic line in Xrays. This kind of signature is within reach of satellite detectors like CHANDRA and XMN which have put strong limits on the couplings between sterile and active neutrinos (due to the lack of detection of this kind of signal). Recently, the existence of an unidentified line in the combined spectrum of a large set of X-ray galactic clusters has been reported [12] and independently, in the combined observation of the Perseus Cluster and the M31 Galaxy [13]. These observations can be compatible with the decay of a sterile neutrino with a mass of approximately 7 keV. Confirmation of the latter result requires further observation, and most probably, higher resolution detectors like the forthcoming Astro-H. As we will show in the analysis, the "(2,3) ISS" model can account for this intriguing possibility; however, we will only impose that the sterile neutrino lifetime does not exceed current observational limits.

DM generation mechanism:
The second issue to address is to provide a DM generation mechanism accounting for the experimental value of its abundance. In the pioneering work by Dodelson-Widrow (DW) [7], it has been shown that the DM abundance can be achieved through active-sterile neutrino transitions 2 . This kind of production is always present provided that there is a non-vanishing mixing between active and sterile neutrinos; as a consequence, it is possible to constrain the latter as function of the neutrino mass by imposing that the DM relic abundance does not exceed the observed value. The "(2,3) ISS" framework allows for an additional production mechanism, consisting in the decay of the heavy pseudo-Dirac states. We will discuss this point at a subsequent stage.
Limits from structure formation: Sterile neutrinos in the mass range relevant for the "(2,3) ISS" model are typically classified as warm dark matter. This class of candidates is subject to strong constraints from structure formation, which typically translate into lower bounds on the DM mass. We notice however, that the warm nature of the DM is actually related to the production mechanism determining the DM distribution function. Sterile neutrinos -with masses at the keV scale -produced by the DW mechanism can be considered as WDM; this may not be the case for other production mechanisms.
In the next section we will investigate whether the "(2,3) ISS" can provide a viable DM candidate.
3 Dark matter production in the "(2,3) ISS" In this section we address the impact of the combination of three kinds of requirements on the DM properties on the "(2,3) ISS" parameter space. The results presented below rely on the following hypothesis: a standard cosmological history is assumed with the exception of possible effects induced by the decays of heavy neutrinos; only the interactions and particle content of the "(2,3) ISS" extension of the SM are assumed.
Regarding DM production we will not strictly impose that the relic abundance reproduces the observed relic abundance, Ω DM h 2 ≈ 0.12 [25], but rather determine the maximal allowed DM fraction f WDM within the framework of the "(2,3) ISS" parameter space.
The main production mechanism for DM is the DW, which is present as long as mixing with ordinary matter is switched on. In addition, the DM could also be produced by the decays of the pseudo-Dirac neutrinos. However, a sizeable contribution can only be obtained if at least one of the pseudo-Dirac states lies in the mass range 130 GeV -1 TeV. Moreover, the pseudo-Dirac states can also have an indirect impact on the DM phenomenology since, under suitable conditions, they can release entropy at their decay, diluting the DM produced by active-sterile oscillations, as well as relaxing the bounds from structure formation. As will be shown below, this effect is also restricted to a limited mass range for the pseudo-Dirac neutrinos.
For the sake of simplicity and clarity, we first discuss the case in which the heavy pseudo-Dirac neutrinos can be regarded as decoupled, and discuss at a second stage their impact on DM phenomenology.

Dark matter constraints without heavy neutrino decays
We proceed to present the constraints from dark matter on the "(2,3) ISS" model, always under the hypothesis that heavy neutrinos do not influence DM phenomenology.
Regarding the relic density, for masses of the lightest-sterile neutrino with mass m s > 0.1 keV, we use the results 3 of [26]: C α are active flavour-dependent coefficients 4 which can be numerically computed by solving suitable Boltzmann equations. In the case of a sterile neutrino with mass m s < 0.1 keV, we have instead used the simpler expression [9]: where sin 2 2θ = 4 α=e,µ,τ |U αs | 2 , with |U αs | being the active-sterile leptonic mixing matrix element. We have then computed the DM relic density using Eqs. (8,9) for a set of "(2,3) ISS" configurations satisfying data from neutrino oscillation experiment and laboratory constraints.
We have imposed f WDM = Ω DM /Ω Planck DM ≤ 1 thus obtaining constraints for m s and U αs . The configurations with DM relic density not exceeding the experimental determination have been confronted with the limits coming from structure formation. There are several strategies to determine the impact of WDM on structure formation, leading to different constraints; in fact most of these constraints assume that the total DM component is accounted by WDM produced through the DW mechanism. Notice that these constraints can be relaxed when this hypothesis does not hold and we will address this point in a forthcoming section.
In the following, and when possible, we will thus reformulate the bounds from structure formation in terms of the quantity f WDM which represents the amount of DM produced from active-sterile oscillation 5 .
The most solid bounds come from the analysis of the phase-space distribution of astrophysical objects. The WDM free-streaming scale is of the order of the typical size of galaxies; as a consequence, the formation of DM halos, as well as that of the associated galaxies is deeply influenced by the DM distribution function. According to this idea, it is possible to obtain robust limits on the DM mass by requiring that the maximum of the dark matter distribution function inferred by observation, the so called coarse grained phase space density, does not exceed the one of the fine-grained density, which is theoretically determined and dependent on the specific DM candidate. Using this method, an absolute lower bound on the DM mass of around 0.3 keV, dubbed Tremaine-Gunn (TG) bound [27] was obtained by comparing the DM distribution from the observation of Dwarf Spheroidal Galaxies (Dphs) with the fine-grained distribution of a Fermi-gas. A devoted study of sterile neutrinos produced by DW mechanism has been presented in [28], where a lower mass bound of the order of 2 keV was obtained. This limit can be evaded assuming that the WDM candidate is a subdominant component, while the DM halos are mostly determined by an unknown cold dark matter component. The reformulation of the limits in this kind of scenarios requires a dedicated study (an example can be found in [29]). In this work we conservatively rescale the results of [28] under the assumption that the observed phase-space density is simply multiplied by a factor f WDM . Moreover we have considered as viable the points of the "(2,3) ISS" model with m s < 2 keV, featuring a value f WDM 1%, which corresponds approximatively to the current experimental uncertainty in the determination of the DM relic density.
For masses above 2 keV another severe bound is obtained from the analysis of the Lyman-α forest data. From these it is possible to indirectly infer the spectrum of matter density fluctuations, which are in turn determined by the DM properties. The Lyman-α constraint is strongly model dependent and the bounds are related to the WDM production mechanism, and to which extent this mechanism contributes to the total DM abundance. In order to properly take into account the possibility of only a partial contribution of the sterile neutrinos to the total DM abundance, we have adopted the results presented in [30] where the Lyman-α data have been considered in the case in which sterile neutrinos WDM account for the total DM abundance, as well as in the case in which they contribute only to a fraction (the remaining contribution being originated by a cold DM component). More precisely, we have considered the most stringent 95% exclusion limit 6 , expressed in terms of (m s , f WDM ), and translated it into an exclusion limit on the parameters of our model 7 , namely the mass m s of the sterile neutrino and its effective mixing angle with active neutrinos θ s . We finally remark that WDM can be constrained also through other observations, as the number of observed satellites of the Milky way [32][33][34], giving a lower bound on the DM mass 5 The results presented are in fact approximative estimates. A proper formulation would require detailed numerical studies, beyond the scope of this work. 6 The limit considered actually relies on data sets which are not up-to-date. A more recent analysis [31] has put forward a stronger limit in the case of a pure WDM scenario, and thus the limits are underestimated. As it will be clear in the following, the final picture is not affected by this. 7 Notice that the Lyman-α method is reliable for DM masses above 5 keV. For lower values there are very strong uncertainties and it is not possible to obtain solid bounds. In [30] it is argued that the limit on fWDM should not significantly change at lower masses with respect to the one obtained for neutrinos of 5 keV mass. of approximately 8.8 keV. This last kind of limits however strictly relies on the assumption that the whole dark matter abundance is totally originated by a WDM candidate produced through the Dodelson-Widrow mechanism and cannot straightforwardly be reformulated in case of a deviation from this hypothesis; thus we have not been considered these limits in our study.
The inverse seesaw realisations passing the structure formation constraints have to be confronted to the limits from the X-ray searches, as reported in, for instance [33]. The corresponding constraints are again given in the plane (m s , θ s ), and can be schematically expressed by 8 : where [35]: with i running over the different active neutrino mass eigenstates (3 different final states in the decay are possible). In the above expression we have again accounted for the possibility that the sterile neutrino contributes only partially to the DM component by rescaling the limit with a factor f WDM 9 . The result of the combination of the three kinds of constraints applied in our analysis, namely dark matter relic density, structure formation and indirect detection, is reported in Fig. 1. As can be seen, the requirement of a correct DM relic density has a very strong impact, excluding a very large portion of the parameter space (grey region) at the highest values of the active-sterile mixing angles. Phase space density constraints rule-out most of the configurations with mass of the lightest sterile neutrino below ∼ 2 keV (blue region), a part a narrow strip (green region) corresponding to f WDM < 1%. In this last region, and although not ruled out, the "(2,3) ISS" model cannot solve the Dark Matter puzzle, at least in its minimal realisation. In the large mass region, namely above 2 keV, a further exclusion comes form Lyman-α and indirect detection bounds (respectively black and yellow region) reducing the allowed active-sterile mixing. A sizeable contribution to the DM relic density can be thus achieved in a small localised region (in red) of the parameter space, corresponding to masses of the lightest sterile neutrino in the range 2 -50 keV and for active-sterile mixing angles 10 −8 sin 2 2θ 10 −11 . We show in the right panel of Fig. 1 the maximal value of f WDM allowed by the cosmological constraints as function of the DM mass. As can be seen, the lightest sterile neutrino can only partially account for the DM component of the Universe with f WDM ∼ 0.43 in the most favourable case. The maximal allowed DM fraction increases for the lowest values of the mass until a maximum at around 7 keV, after which it displays a sharp decrease. This behaviour can be explained as follows: at lower masses, the Lyman-α bounds are the most effective and become weaker as the mass of the sterile neutrino increases, thus allowing The grey region corresponds to a DM relic density exceeding the cosmological value. The blue, black and yellow regions are also excluded by phase space distribution, Lyman-α and X-ray searches constraints, respectively. The green region corresponds to configurations not excluded by cosmology but in which the lightest sterile neutrino contributes with a negligible amount to the DM relic density. Finally, the red region corresponds to the "(2,3)ISS" configurations fulfilling all the cosmological constraints, and for which the contribution to the dark matter relic density from the light sterile neutrino is sizeable. On the right panel, maximal value of f WDM allowed by cosmological constraints as a function of the mass of the lightest sterile neutrino.
for larger f WDM . At the same time, the bounds from X-ray sources become stronger (since higher masses imply higher decay rates) thus reducing the allowed DM fraction as the mass increases.
Notice that the above analysis is valid within the assumption that the production of the lightest sterile neutrino occurs in the absence of a lepton asymmetry. Indeed, as firstly shown in [36], the production of sterile neutrinos can be resonantly enhanced (as opposed to the conventional DW production usually called non-resonant) in presence of a non-zero lepton asymmetry. In this case the correct dark matter abundance is achieved for much smaller active-sterile mixing angles, thus evading the limits from dark matter indirect detection; in addition, the resonant production alters the DM distribution function with respect to a non-resonant production, rendering it "colder" and thus compatible with Lyman-α constraints [37].
Interestingly a lepton asymmetry can be generated in frameworks featuring keV scale sterile neutrinos accompanied by heavier right-handed neutrinos. The entries of the active-sterile mixing matrix can in general be complex, and give rise to CP-violating phases; as a consequence, a lepton asymmetry can be generated by oscillation processes of the heavy neutrinos. In particular, it has been shown that a pair of quasi-degenerate right-handed neutrinos with masses of the order of a few GeV can generate a lepton asymmetry before the EW phase transition (which is converted to the current baryon asymmetry of the Universe) and then at much later times, the lepton asymmetry needed to provide the correct relic density for a keV scale sterile neutrino [38][39][40][41]. The "(2,3) ISS" model also features pairs of quasi-degenerate heavy neutrinos which can be of the correct order of mass. However, the lepton asymmetry needed to ensure the correct DM relic density, compatible with the bounds discussed, requires an extreme degeneracy in the heavy neutrino spectrum, of the order of the atmospheric mass differences. Such an extreme degeneracy is not achievable for the ISS model since the predicted degeneracy of the pair of heavy neutrinos is of O(µ), corresponding to around 1 keV for the cases under consideration 10 .

Impact of the heavy pseudo-Dirac states
The picture presented above can be altered in some regions of the parameter space due to the presence of the heavy neutrinos. Indeed, contrary to the DM candidate, they can exist in sizeable abundances in the Early Universe owing to their efficient Yukawa interactions, and influence the DM phenomenology through their decays. There are two possibilities. The first one is direct DM production from decays mediated by Yukawa couplings. The branching ratio of these processes is small when compared to that of other decay channels into SM states, since it is suppressed by the small active-sterile mixing angle, an efficient DM production can nevertheless be achieved through the so called freeze-in mechanisms if the pseudo-Dirac neutrinos are heavier than the Higgs boson. Significantly lighter pseudo-Dirac neutrinos, namely with masses below ∼ 20 GeV, can instead indirectly affect DM phenomenology. Indeed, they can be sufficiently long-lived such that they can dominate the energy density of the Universe, injecting entropy at the moment of their decay. We will discuss separately these two possibilities in the next subsections.

Effects of entropy injection
The conventional limits on sterile neutrino DM can be in principle evaded in presence of an entropy production following the decay of massive states dominating the energy density of the Universe [42]. A phase of entropy injection dilutes the abundance of the species already present in the thermal bath and, in particular, the one of DM if such an entropy injection occurs after its production. In addition, the DM momentum distribution gets redshifted -resembling a "colder" DM candidate -and suffering weaker limits from Lyman-α. This phase of entropy injection can be triggered in the "(2,3) ISS" model by the decay of the heavy pseudo-Dirac neutrinos if the following two conditions are realised: Firstly, at least some of the heavy sterile neutrinos should be sufficiently abundant to dominate the energy budget of the Universe. Secondly they must decay after the peak of dark matter production, but before the onset of Big Bang Nucleosynthesis (BBN). These two requirements will identify a limited region of the parameter space outside which the results of the previous subsection strictly apply.
All the massive eigenstates have Yukawa interactions with ordinary matter described by an effective coupling Y eff which is defined by: These interactions are mostly efficient at high temperature when scattering processes involving the Higgs boson and top quarks are energetically allowed; in addition, they maintain the pseudo-Dirac neutrinos in thermal equilibrium until temperatures of the order of ∼ 100 GeV, provided that Y 2 eff 10 −14 [43]. If this condition is satisfied, an equilibrium abundance of heavy pseudo-Dirac neutrinos existed at the early stages of the evolution of the Universe.
The Yukawa interactions become less efficient as the temperature decreases. At low temperature the transition processes from the light active neutrinos become important. For a given neutrino state, the rate of the transition processes reaches a maximum at around [44]: The transition rate of each neutrino at the temperature T max,I exceeds the Hubble expansion rate H if [44]: and thus, if this condition is satisfied, the corresponding pseudo-Dirac neutrinos are in thermal equilibrium in an interval of temperatures around T max,I . Notice that the picture depicted above assumes that the production of sterile neutrinos from oscillations of the active ones is energetically allowed; as a consequence it is valid only for neutrino masses lower than T max : As will be made clear in the following, neutrinos heavier than M I,max have excessively large decay rates to affect DM production and hence will not be relevant in the subsequent analysis. In Fig. 2, we present the typical behaviour of the pseudo-Dirac neutrinos in the regimes of high and low temperatures, which are dominated, respectively, by Yukawa interactions and active-sterile transitions. In the left panel of Fig. 2, we display the values of the mass and effective Yukawa couplings (Y eff ) of the lightest pseudo-Dirac state (the other heavy states exhibit an analogous behaviour), corresponding to a set of "(2,3) ISS" realisations compatible with laboratory tests of neutrino physics (red points). The green region translates the equilibrium condition for the Yukawa interactions. In the larger mass region, i.e. for masses significantly larger than 10 GeV, the value of the effective Yukawa coupling Y eff is always above the equilibrium limit and can even be of order one for higher values of the mass. In this region, the pseudo-Dirac neutrinos can have a WIMP-like behaviour and can be in thermal equilibrium until low temperatures. In the intermediate mass region, i.e. for masses between 1 and a few tens of GeV, equilibrium configurations are still present. However the values of Y eff are lower with respect to the previous case and the decoupling of the pseudo-Dirac neutrinos depends on the oscillation processes at low temperatures. Configurations for which Y eff is too small to ensure the existence of a thermal population of pseudo-Dirac neutrinos in the early Universe (they can be nonetheless created by oscillations at lower temperatures) are also present. This last kind of configurations are the only ones corresponding to masses below 0.1 GeV. We emphasise that the outcome discussed here is a direct consequence of the "(2,3) ISS" mechanism which allows to generate the viable active neutrino mass spectrum for pseudo-Dirac neutrinos with masses of the order of the EW scale, and for large values of their Yukawa couplings. For comparison, we display in the same plot the distribution of values of the effective Yukawa couplings of the WDM candidate as a function of the mass of next-to-lightest sterile state m 5 (blue points). As can be seen, the corresponding solutions are always far from thermal equilibrium due to the suppressed mixing U ν R ,4 . Figure 2: On the left panel: effective Yukawa couplings Y eff for the neutrino DM candidate (blue points) and of the lightest pseudo-Dirac particle (red points), as a function of the mass m 5 . The green region corresponds to values Y eff > √ 2 × 10 −7 , the limit above which the states are in thermal equilibrium. On the right panel: mixing of the electron neutrino with the lightest pseudo-Dirac state as a function of its mass. The yellow region corresponds to the kinematically forbidden values of the sterile mass, see Eq. (15). The red region denotes the solutions not in thermal equilibrium.
In the right panel of Fig. 2, we display the mixing (for small angles it is possible to approximate θ e5 U e5 ) of the lightest pseudo-Dirac state with the electron neutrino as a function of m 5 for the ISS realisations compatible with laboratory limits. The yellow region corresponds to the values of the sterile mass for which the DW production mechanism is kinematically forbidden (see Eq. (15)). The red region denotes the solutions which are not in thermal equilibrium.
Combining the results obtained from the two panels of Fig. 2, we can conclude that all the considered realisations in the relevant mass interval satisfy the equilibrium conditions. Consequently we can always assume the presence of an equilibrium population of the pseudo-Dirac states up to temperatures of the order of T max,I . We stress that T max,I is not the actual decoupling temperature that has been instead determined in for instance [38] and more recently in [41], and which turns out to be lower than T max ; however this affects only marginally our discussion.
As already pointed out, we will be interested in masses of the pseudo-Dirac neutrinos not exceeding 10 -20 GeV. For such a mass range we can safely assume that the neutrinos decouple when they are relativistic (see Eq. (13)) and that their decay occurs at a much later stage, when they become non-relativistic, as described in [44].
In this setup, the pseudo-Dirac states can dominate the energy budget of the Universe if their energy density, which is defined by exceeds the radiation energy density ρ r = π 2 30 g * (T )T 4 , where g * (T ) represents the number of relativistic degrees of freedom at the temperature T . Provided that the pseudo-Dirac neutrinos are sufficiently long-lived, this occurs at a temperature T given by: where we have taken g * (T D ) = 86.25 and m 5 is the mass of the lightest pseudo-Dirac neutrino.
In this scenario the decay of the pseudo-Dirac neutrinos is accompanied by a sizeable amount of entropy; the conventional radiation dominated era restarts at the reheating temperature T r,I [42], and the abundance of the species present in the primordial thermal bath is diluted by a factor S, which is defined as the ratio of the entropy densities of the primordial plasma at temperatures immediately below and above the reheating one. Notice that the above discussion corresponds to a simplified limit: in general the four pseudo-Dirac neutrinos have different masses and different lifetimes. In the "(2,3) ISS" model the pseudo-Dirac states appear as pairs with the mass splitting in each pair much smaller than the masses of the corresponding states. Identifying the mass scale of each pair as m and M , with m < M , we can write, to a good approximation 11 : where S m and S M are the dilution factors associated to the decays of the two pairs of pseudo-Dirac states, occurring at the two reheating temperatures T r,M and T r,m , given by: where Γ m,M is the decay rate of the heavy neutrinos. Notice that in the last term of each of the above equations, the effects of the first entropy dilution have been included in the abundance of the lightest pair of heavy neutrinos. The DM phenomenology is affected only when the pseudo-Dirac neutrinos dominate the Universe and decay after DM production. For keV scale DM, this translates into the requirement T r,m 150 MeV. On the other hand, a very late reheating phase would alter the population of thermal active neutrinos, leading to modifications of some quantities such as the primordial Helium abundance [46] and the effective number of neutrinos N eff , and producing effects in structure formation as well. By combining BBN and CMB data 12 it is possible to determine a solid bound T r,m > 4 MeV [48]. In addition, we have considered a (relaxed) limit of T r,m > 0.7 MeV by taking into account the possibility that this bound is evaded when the decaying state can produce ordinary neutrinos [49]. This choice is also motivated by the fact that after the decay of the neutrinos with mass M , the ratio ρ I /ρ r between the energy densities of the remaining neutrinos and of the radiation is of order 3 -5, implying that although subdominant, the radiation component is sizeable.
The requirement 4(0.7) MeV ≤ T r,m ≤ 150 MeV is satisfied only for a very restricted pseudo-Dirac neutrinos mass range. Indeed, sterile neutrinos can decay into SM particles through threebody processes mediated by the Higgs boson with a rate: which implies an excessively short lifetime for sterile neutrinos unless their masses are below (approximately) O(10 GeV), in such a way that the decays into third generation quarks are kinematically forbidden and Y eff can assume lower values. At these smaller masses, a sizeable contribution comes from Z mediated processes with a rate: where we have defined, for simplicity, effective mixing angles θ I , I = m, M between the pseudo-Dirac neutrinos and the active ones. We have reported in Fig. 3 the limit values of the lower reheating temperature T r,m as a function of the mass scale m and the effective mixing angle θ m , for three values of Y eff , namely 0.1, 10 −3 and 10 −6 . The regions above the red curves correspond to an excessively large reheating temperature which does not affect DM production. The light-grey (dark-grey) regions below the blue curves represent values of the reheating temperature in conflict with the conservative (relaxed) cosmological limit of 4 (0.7) MeV. For "natural", i.e. O(1), values of the effective coupling, the decay rate of the heavy neutrinos is dominated by the Higgs channel and tends to be too large except for a narrow strip at masses of 1 − 2 GeV. At lower values of Y eff the size of the region corresponding to the interval 4 (0.7) -150 MeV of reheating temperatures increases. The contours corresponding to the lower values of the reheating temperature are mildly affected by the values of Y eff since, at lower masses, the Higgs channel is suppressed by the Yukawa couplings of the first generation, in comparable amount with the Z channel. As can be seen from the left panel of Fig. 2, laboratory constraints favour values of Y eff < 10 −3 in the mass range 1 − 20 GeV thus favouring the possibility of an impact of the heavy neutrinos decays on the DM phenomenolog, in this region.
In Fig. 4 we have estimated the range of values of S in the allowed parameter space, see Eq. (18). In order to illustrate, we have chosen to vary m and θ m fixing M to be M = 2 × m and θ M = 10 −4 , which corresponds to a viable case for the "(2,3) ISS" mass spectrum; we have also fixed the effective Yukawa coupling as Y eff = 10 −6 in order to maximize the phenomenologically relevant region of the parameter space.
As can be seen from Fig. 4, we have only very moderate values of the entropy injection when the conservative lower limit of 4 MeV is imposed on the reheating temperature; values of S up to around 20 can be achieved once a weaker bound is considered.
Having determined the range of variation of the entropy dilution within the "(2,3) ISS" parameter space, we have reformulated the limits on the DM mass and mixing angle, as presented in the previous section, for the case where S > 1. The limits from DM relic density can be straightforwardly determined by simply rescaling it by a factor 1/S. The limits from Lyman-α are more difficult to address since this would require a different analysis (as for example [30,37]), which is computationally demanding and lies beyond the scope of this work. To a good approximation, one can assume a redshift factor of S 1/3 [44] for the DM momentum distribution and translate it into a modified limit for the DM mass given by m S s,Lyα ≥ m s,Lyα /S 1/3 , where m s,Lyα is the lower limit on the DM mass for a given value f WDM of the DM fraction, in the case where S = 1.
The X-ray limits remain unchanged with respect to the previous section since these rely on the DM lifetime. Nevertheless, the entropy injection leads to an indirect effect since a given pair (m s , θ) now corresponds in general to a lower relic density.
The outcome of our analysis is summarized in Fig. 5. In the left panel of Fig. 5, we display the cosmologically favoured region for several values of S ≤ 50 as compared to the case S = 1, represented by red points. Values of S larger than ∼ 20 are not within reach in the framework of the present model, but we have nonetheless extended our analysis up to these values in order to infer the maximal extension of the parameter space which could be achieved. The grey points are excluded by DM constraints unless its relic density is negligible (see previous subsection). As one can see, for S > 1 we have a larger range of allowed values for the active-sterile mixing angle; interestingly, the augmentation of S has a finite effect in enlarging the available parameter space. The grey (light grey) region corresponds to the lowest reheating temperature below 4 (0.7) MeV and is thus in tension with the cosmological bounds. In the region above the red curve the entropy injection takes place before the DM production.
On the right panel of Fig. 5, we display the maximal value of f WDM for several values of S (we also display for comparison the case corresponding to S = 1). As one can see, there is only a marginal increase, namely from 0.43 to 0.48, of the maximal allowed value of f WDM . On the other hand, the maximal DM fraction is achieved for smaller values of the allowed DM mass, namely ∼ 2 keV, as opposed to values around ∼ 7 keV in the case where S = 1. We finally notice that the maximal DM fraction is not an increasing function of S but on the contrary, a maximum achieved at S = 20 is followed by a sharp decrease. The reason of such a behaviour is mostly due to the X-ray exclusion. Indeed, as already pointed out, any fixed value f WDM imposes a condition on (m s , sin 2 2θ) which is not sensitive to the mechanism accounting for the DM generation (more specifically, the value of S in our case). Since the dark matter generation mechanism also depends on (m s , sin 2 2θ), the interplay with the X-ray exclusion, as well as the effect of entropy injection, favours larger mixing angles (thus maximizing the production of dark matter) and lower values of the mass (which in turn minimize the DM decay rate). Our analysis shows that f WDM = 1 could be achieved for S 10. This is not sufficient to relax the Lyman-α bound down to the value m s = 2 keV because of the scaling of the latter limit as S 1/3 .
We finally point out that for very high values of S, sizeable values of f WDM could be achieved for very large mixing angles, already excluded by indirect dark matter detection. This is at the origin of the saturation of the cosmologically favoured region observed in the left panel of Fig. 5.

Dark Matter Production from heavy neutrino decays and the 3.5 keV line
As already mentioned, the pseudo-Dirac neutrinos can produce dark matter through their decays. These processes are mediated by Yukawa interactions and the decay rate is proportional to Y eff sin θ, and thus suppressed with respect to the decay channels into only SM particles by the active-sterile mixing angle. A sizeable DM production can be nonetheless achieved through the so called freezein mechanism [50]. It consists in the production of the DM while the heavy neutrinos are still in thermal equilibrium and, to be effective, requires that the rate of decay into DM is very suppressed, such that it results lower than the Hubble expansion rate. In our setup, this condition can be expressed as: Y eff sin θ < 10 −7 .
The dark matter relic density depends on the decay rate of the pseudo-Dirac neutrinos into DM as follows: where the sum runs over the pseudo-Dirac states and g I represents the number of internal degrees of freedom of each state. For pseudo-Dirac neutrinos lighter than the Higgs boson, DM production occurs through three-body processes whose rate is too suppressed to generate a sizeable amount of DM. On the other hand, the above analytical expression is not strictly applicable for heavier pseudo-Dirac neutrinos since the mixing angle θ depends on the vacuum expectation value (vev) of the Higgs boson and is thus zero above the EW phase transition temperature. To a good approximation, the correct DM relic density is determined by multiplying Eq. (22) by ε 2 (m I ), where the function ε(m I ) is given by: with f (x I ) describing the evolution of the Higgs vev v(T ) with the temperature and which can be in turn approximated, according the results presented in [51], by: where T EW ≈ 140 GeV is the temperature associated to the EW phase transition. As shown in Fig. 6, the function ε 2 (m I ) sharply decreases with the mass of the pseudo-Dirac neutrino since most of the FIMP (Feebly Interacting Massive Particle) production occurs around the mass of the decaying particle. As a consequence, we can have sizeable production of DM only for masses of the decaying particles not too much above the scale of the electroweak phase transition while DM production is negligible for masses of the pseudo-Dirac neutrinos above the TeV scale. Using the expression of the rate associated to the process N I → h + DM: the DM relic density is given by: It is then clear that the correct DM relic density can be achieved with a suitable choice of the parameters. It is worth noticing that this production mechanism is complementary to the DW one, which is always active provided that there is a nonzero active-sterile mixing. We have reported in Fig. 7 the (observed) value Ω DM h 2 = 0.12 of the DM abundance, assuming for simplicity the same mass m 5 and effective Yukawa couplings Y eff for the 4 heavy pseudo-Dirac states, for different values of the DM mass and considering the maximal value of sin θ allowed by cosmological constraints -including thus the corresponding contribution from DW production mechanism-. The displayed red points correspond to configurations of the "(2,3) ISS" model in agreement with all laboratory constraints. Those configurations corresponding to pseudo-Dirac states far from thermal equilibrium, and thus not accounting for a freeze-in production mechanism, are delimited by a blue region. The shape of the lines can be understood as follows: for pseudo-Dirac masses comparable with the Higgs one, the kinematical suppression in Eq. (26) is significant, requiring sizeable Yukawas; for m I 200 GeV the dependence on m I is weaker, and the curve reaches a plateau, while for m I 400 GeV the suppression due to the function ε (m I ) becomes significant requiring larger Yukawas, eventually violating the freeze-in condition Y eff sin θ < 10 −7 for m I 1.0 TeV. Figure 7: Viable configurations (continuous lines) for the heavy pseudo-Dirac masses m 5 and the corresponding effective Yukawa couplings Y eff accounting for the observed dark matter abundance of light sterile neutrinos via the freeze-in production mechanism, with masses and mixings for the sterile neutrinos compatible with cosmological bounds. The red points denote different realisations of the "(2,3) ISS" model. In the blue region the production is not effective since the pseudo-Dirac states are out of thermal equilibrium. The lines end when the condition Y eff sin θ < 10 −7 is violated. The yellow line accounts for the still unidentified monochromatic 3.5 keV line in galaxy cluster spectra [12,13].
The requirement of light sub-eV active neutrino masses together with µ ≈ keV and M R ≈ v, implies values for the Yukawa couplings in the appropriate range to accounting for the observed DM abudance (m ν ≈ µY 2 v 2 /M 2 R , see Eq. (7)). We emphasize here that this is not the case for a type-I seesaw realisation since in this case the relation m ν ≈ Y 2 v 2 /M R < 1 eV implies Y 10 −6 if M R ≈ v, and thus the contribution from the freeze-in process is not sufficient to account for the total DM abundance.
Among the lines displayed in Fig. 7, we have highlighted in yellow the one corresponding to the following DM mass and mixing angle, which can account for the monochromatic 3.5 keV line observed in the combined spectrum of several astrophysical objects [12,13]. The results presented in Fig. 7 do not take into account the possible constraints from structure formation. As will be made clear in the next section, the limits discussed above should be sensitively relaxed since the DM produced through the freeze-in mechanism has a "colder" distribution with respect to the DW mechanism. A reformulation of the corresponding limits is beyond the scope of this work, especially in the case in which the DM production receives sizeable contributions both from DW and decay of the pseudo-Dirac neutrinos. We argue nonetheless that the parameters accounting for the keV line can be compatible with bounds from structure formation since for this choice (of parameters), the DM abundance is entirely determined by the decay of the heavy neutrinos (the DW contribution for that value of the mixing angle is less than 4%) and the corresponding distribution function is "colder" with respect to the one of a resonantly produced DM, which results compatible with the observational limits [52].
To summarize the results obtained and discussed in this section, one can state that in the absence of effects from the heavy pseudo-Dirac neutrinos, the "(2,3) ISS" model can, in the most favourable case, account for to approximatively the ∼ 43% of the total DM density for a mass of approximatively 7 keV. This percentage slightly increases up to 48%, for a DM mass of around 2 keV, once accounting for an entropy dilution factor of 5 -20 which can be possible for masses of the pseudo-Dirac neutrinos of 3 -10 GeV. The total DM component can be accounted for only in the region m h < m I < 1 TeV, when the DM can be produced through the freeze-in mechanism, although the compatibility with structure formation should be still addressed. In order to also reproduce the correct relic density for masses of the sterile neutrinos below the Higgs boson mass, it is necessary to extend the particle field content of the model; for this purpose, we will propose in the following section a minimal extension of the "(2,3) ISS" model.

Dark Matter Production in minimal extension of the "(2,3) ISS" model
In order to achieve the correct dark matter relic density in the pseudo-Dirac states low mass regime, we consider a minimal extension of the "(2,3) ISS" model. This consists in the introduction of a scalar field Σ, singlet under the SM gauge group, interacting only with the sterile fermionic states and the Higgs boson. There are of course several other possibilities, see for instance [53]. In this minimal extension, the part of the Lagrangian where the new singlet scalar field is involved reads: We consider that the field Σ has a non-vanishing vev Σ that would be at the origin of the Maiorana mass coupling µ which can thus be expressed as: For simplicity we will limit the scalar potential to the following terms (see e.g. [54] for a more general discussion): Following a pure phenomenological approach, we will consider values of the portal coupling λ HΣ from order of 10 −2 , corresponding to limits from effects on the Higgs width [55], down to very low values, i.e. O 10 −8 or 10 −9 (see for instance [56] and references therein for some examples of theoretically motivated models with extremely suppressed λ HΣ ). We will assume for simplicity that the scalar singlet field is heavier that the Higgs boson, m Σ > 200 GeV, and assume m Σ ≤ Σ in order to avoid non perturbative values of λ HΣ .
The DM density is generated by the decay of Σ and is thus tied to the abundance of the latter, which in turn depends on the efficiency of the process ΣΣ ↔ hh triggered by the portal like coupling λ HΣ . A proper description of the DM density requires the resolution of a system of coupled Boltzmann equations for the DM number density, as well as for the abundance of the Σ field and possibly for the heavy pseudo-Dirac neutrinos, which also interact with Σ -also including effects of entropy release. Further details of this computation can be found in the appendix.
In the following we will present analytical expressions which describe, to a good approximation, the DM production mechanism. For simplicity we will assume that the pseudo-Dirac neutrinos are in thermal equilibrium (the case of non-equilibrium configurations substantially coincides with the studies already presented in [54,56]) and with lifetimes such that the effects of entropy injection are not relevant. As will be made clear, pseudo-Dirac neutrinos have a non trivial impact on DM production. We will thus for definiteness discuss two specific mass regimes, namely the case in which all the pseudo-Dirac neutrinos are lighter than Σ and the case in which they have instead similar or greater masses.
At high enough values of λ HΣ , the pair annihilation processes ΣΣ ↔ hh maintain the field Σ into thermal equilibrium 13 . Indeed, by comparing the 2 → 2 rate, associated to the thermally averaged cross-section σv ∼ 10 −2 × λ 2 HΣ m 2 Σ , with the Hubble expansion rate, the field Σ can be considered to be in thermal equilibrium in the Early Universe for λ HΣ ≥ λ HΣ where: On the contrary, its decay rate into DM is always suppressed compared to the Hubble rate due to the low value of the couplings h αα (see Eq. (29)). The DM can thus be produced through the freeze-in mechanism from the decays of Σ and its corresponding abundance can be expressed as: where: is an effective coupling taking into account all the decays Σ → N I + DM, I = 4, 8 which are kinematically open 14 . The contribution to the relic DM density reads: On general grounds, out-of-equilibrium -i.e. after chemical decoupling -decays of Σ may also contribute to DM production and the corresponding contribution to the DM density can be schematically expressed as: In the above equation, T f.o. is the standard freeze-out temperature of Σ and b I represents the number of DM particles produced per decay for a given decay channel. The branching ratio of the decay of Σ into DM is given by: where sin α ∝ λ HΣ represents the mixing between Σ and the Higgs boson 15 . Due to the very low couplings h αα , the branching ratio of the decay of Σ into DM is very suppressed with respect to the branching ratio of the decay into two fermions induced by the mixing with the Higgs boson, even for low values of the mixing itself. Furthermore, the total lifetime of the scalar field is comparable with the freeze-out timescale. Consequently, the out-of-equilibrium production is sizeable for λ HΣ ∼ λ HΣ when the scalar field features an early decoupling, i.e.
Finally, the DM relic density can be estimated as: In the case in which λ HΣ λ HΣ (see Eq. (31)), the field Σ is too feebly interacting to be in thermal equilibrium in the Early Universe. It can be nonetheless produced in sizeable quantities by freeze-in and then decay through out-of-equilibrium processes [56]. The field Σ is produced by the 2 → 2 processes mediated by the portal interactions as well as by the 2 → 1 processes, N I N I → Σ, I = 5, 8 if the pseudo-Dirac neutrinos are lighter than Σ. The abundance of Σ can be expressed as: with The two terms inside the parenthesis refer to the contributions from the 2 → 1 (where the sum over I, J runs over the pseudo-Dirac neutrinos in thermal equilibrium) and the 2 → 2 processes, respectively. The DM abundance is given, analogously to Eq. (35), by B Y SFI Σ . In the case where the pseudo-Dirac neutrinos are heavier than Σ, the DM generation process in the regime λ HΣ ≥ λ HΣ proceeds along the same lines as described before. There is however the additional contribution from the decays N I → h DM , given by Eq. (26) as well as a further freeze-in contribution from the decays N I → ΣDM given by: In the regime where λ HΣ ≤ λ HΣ , the 2 → 1 production channel for the field Σ is replaced by the production from the decays of the pseudo-Dirac neutrinos through the processes, if kinematically open, N I → N J Σ, I = 5, 8, J = 4, I − 1. In this scenario the abundance of Σ reads: The validity of the assumptions leading to the analytical approximations above has been confirmed by (and complemented) by numerically solving the Boltzmann equations for the system Σ−DM abundances. We display in Fig. 8 the evolution of the abundances of Σ and of the DM, for a set of values of the the coupling λ HΣ and for fixed values of m s , m Σ and Σ to 5 keV, 500 GeV and 1 TeV, respectively. On the left panel the masses of the pseudo-Dirac pairs have been fixed to, respectively, 10 and 20 GeV, while on the right panel the chosen values are 500 GeV and 1 TeV. For this last case we have fixed the coupling Y eff of the pseudo-Dirac neutrinos with the Higgs boson and the DM-active neutrino mixing angle θ to, respectively, 0.01 and 10 −6 in such a way that the freeze-in production from the decays N I → h + DM gives a subdominant contribution, not exceeding 30%.
At higher values of λ HΣ , the abundance of Σ traces its equilibrium value and the DM production occurs prevalently through the freeze-in mechanism. At lower values of λ HΣ the out-of-equilibrium production becomes important thus increasing the total DM relic density.
For λ HΣ < λ HΣ , the abundance of Σ does not follow the equilibrium value but is an increasing function of time (as consequence of the freeze-in production from the 2 → 2 scatterings as well as from the 2 → 1 processes, for pseudo-Dirac neutrinos lighter than Σ, or from decays of the pseudo-Dirac neutrinos themselves in the opposite case) until its decay, which occurs at later timescales compared to the case of high values of λ HΣ . We report in Fig. 9 the contours of the cosmological value of the DM relic density in the (m s , m Σ ) plane for several values of the coupling λ HΣ and for the masses of the pseudo-Dirac neutrinos considered in Fig. 8. As already pointed out, for λ HΣ = 10 −3 the DM relic density is determined by the freeze-in mechanism and thus increases with the DM mass while decreasing with respect to m Σ . For λ HΣ = 10 −6 the out-of-equilibrium production is instead the dominant contribution implying that Ω s h 2 ∝ m s m Σ . For λ HΣ = 10 −7 and λ HΣ = 10 −9 , in the case of light pseudo-Dirac neutrinos, the relic density is again proportional to the ratio m s /m Σ , as expected from Eq. (38). In the case of heavy pseudo-Dirac neutrinos, the regime λ HΣ λ HΣ is substantially dominated by the freeze-in production of Σ from the decays of the pseudo-Dirac neutrinos and its subsequent out-of-equilibrium decays. The dependence on m Σ shown in Fig. 9 is due to the kinematical factor in Eq. (40). We notice that the correct DM relic density, for the chosen set of parameters, is achieved for DM masses between 1 -15 keV. These results can be straightforwardly generalized in the case of entropy production from the pseudo-Dirac neutrinos. Indeed, as can be seen from Fig. 8, the DM production typically occurs at earlier stages compared to the ones at which sizeable entropy production is expected (see previous section). As a consequence the instantaneuous reheating approximation can be considered as valid and we can just rescale the DM relic density by a factor S. In this case, the correct DM relic density is achieved for higher values of the DM masses.
We emphasise, as already done in the previous section, that a complementary contribution to the DM relic density from DW production mechanism is in general present. The DM production related to the decays of Σ allows to achieve the correct relic density without conflicting with the Xrays limits since it does not rely on the mixing with the active neutrinos; this is not the case for the bounds from structure formation. However, applying the limits on DM from structure formation is a very difficult task in our scenario since different DM production channels coexist, originating different dark matter distribution functions. A proper treatment would require to reformulate the bounds case by case by running suitable simulations, which lies beyond the scope of the present work. We will nonetheless provide an approximate insight of how the latter bounds are altered, with respect to the conventional DW production mechanism, by taking some representative examples.
In the following discussion, we consider the case in which the DM is produced by the decays of the field Σ, either through freeze-in or through out-of-equilibrium decays. An approximate reformulation of the limits from structure formation can be obtained by comparing the average momentum of DM at the keV scale with the one corresponding to DW production and by rescaling the lower limit on the DM mass with the shift between these two quantities. The DM distribution function in the various cases of production from decay has been determined in e.g. [54] and [57]). The dark matter produced through freeze-in is typically generated at temperatures of the order of the mass of Σ. Its average momentum depends only on the temperature and can be simply expressed, at temperatures of the order of the keV, as [58]: sensitively lower than the corresponding result (of ∼ 2.83) in the case of DW production. A similar result holds as well in the case of DM produced through freeze-in from the decays of the pseudo-Dirac neutrinos 16 . In the case in which the DM is prevalently produced out-of-equilibrium, the timescale of production varies with the lifetime of Σ and the distribution function tends to be warmer as the latter increases. We report in Fig. 10 the lower limit 17 from Lyman-α on the DM mass, obtained by applying our approximate rescaling to the limit presented in [31], for some scenarios of DM production mechanism, namely freeze-in and out-of-equilibrium production for different decay rates of Σ Figure 10: Lower limit on the DM mass from Lyman-α as a function of the entropy dilution factor S. The limit refers to the cases of dominant freeze-in production from the decay of the scalar field Σ and dominant production from out-of-equilibrium decays for the parameter Λ = 0.01, 0.1, 1 (defined in Eq. (42)). We also report for comparison the corresponding limit in the case of dominant DW production mechanism.
parametrised through the dimensionless quantity: where, As can be noticed, in the most favourable cases, namely freeze-in or out-of-equilibrium production with Λ ≥ 1 (corresponding to m Σ 600 GeV), the limit from Lyman-α is relaxed to approximately 5 keV, in absence of entropy injection, and to further low values for the case S > 1. We remark again that these results assume that all the DM is produced by the decay of Σ (into heavy neutrinos). Interestingly, in the case in which a sizeable contribution from DW production mechanism is also allowed, the DM distribution would feature a Warm and a "colder" component and thus the "(2,3) ISS" model could potentially realise a mixed Cold + Warm DM scenario. This would constitute an intriguing solution to some tensions with observation from structure formation (see e.g. [60] and references therein). This possibility should be thouroughly investigated by means of numerical simulations since the analytical estimates presented above are not valid for multicomponent distributions. This will be thus left for a dedicated study.

Conclusion
In this study we have considered the possibility of simultaneously addressing the dark matter problem and the neutrino mass generation mechanism. We have focused on the truly minimal inverse seesaw realisation -the "(2,3) ISS" model-fulfilling all phenomenological and cosmological requirements and which provides a Warm Dark Matter candidate (for a mass of the lightest sterile state around the keV).
We have conducted a comprehensive analysis taking into account the several possibilities of neutrino mass spectra. In most of the parameter space the DM can be produced only through active-sterile transitions according to the DW production mechanism, accounting, in the most favourable case, for at most ∼ 43% of the relic DM abundance, without conflict with observational constraints. This situation can be improved for two specific choices of the spectrum of the heavy pseudo-Dirac neutrinos. Firstly, one can consider the case of moderately light, i.e. ∼ 1 − 10 GeV, pseudo-Dirac neutrinos. These states can dominate the energy density of the Universe and produce entropy at the moment of their decay, altering the impact of DM on structure formation. However the constraints from dark matter indirect detection are still too severe and the allowed DM fraction is increased only up to ∼ 50%. The second possibility relies upon relatively heavy, ∼ 130 GeV − 1 TeV, pseudo-Dirac pairs, which can produce the correct amount of DM through their decays. In this kind of setup it is also possible for the "(2,3) ISS" to account for the recently reported 3.5 keV line in galaxy cluster spectra. In the final part of this work, we have proposed a minimal extension of the "(2,3) ISS" with the addition of a scalar singlet (at the origin of the lepton number violating masses of the sterile fields) which allows to achieve the correct DM relic density for generic values of the masses of pseudo-Dirac neutrinos. The latter can still participate, at various levels, to the DM production.
a system of two coupled Boltzmann equations whose general form is: +DW.
The first equation traces the time evolution of the field Σ. The first row on the right-hand side represent the decay of Σ, respectively into at least one DM particle and only into pseudo-Dirac neutrinos, if kinematically allowed. Since the latters are assumed in thermal equilibrium this second term is balanced by a term accounting for inverse decays and thus vanishes if Σ is in thermal equilibrium. On the contrary the first term is not balanced by an inverse decay term since the DM has too weak interactions to be in thermal equilibrium and then can be assumed to have a negligible abundace at early stages; this originates the freeze-in production channel. The second row represents the decays, if kinematical allowed, of the pseudo-Dirac neutrinos into Σ and another neutrino. The factor (n I − n I,eq ) assumes that Σ is in thermal equilibrium and disappears if the pseudo-Dirac neutrinos are as well in thermal equilibrium. Here we have again distinguished the decay term into DM, which is non balanced by the inverse process, and the decay term into final thermal states (this distinction holds only if Σ is in thermal equibrium. In the regime λ HΣ λ HΣ the second row should be replaced by the term I Γ N I n I ). The last term finally represents the annihilation processes of Σ. B and B I represent the effective branching fractions of decay of, respectively, Σ and the pseudo-Dirac neutrinos. Γ and σv represent the conventional definitions of the thermal averages [61]: where the function F (x) is determined by numerically solving the integral above. The second equation traces the DM number density. The first two rows represent the DM production from, respectively, Σ and the pseudo-Dirac neutrinos. The term labelled DW represents instead the contribution associated to production from oscillation processes. In the parameter space of interest the two production processes, decay and oscillations, occur at well separated time scales; as a consequence we can drop the DW term from the equations and possibly add its contribution to the final relic density.
In order to account possible effects of entropy injection from the decays of the pseudo-Dirac neutrinos the system above should be completed with a third equation accounting for the non conservation of the entropy (see e.g. [45]). On the other hand it has been shown that the pseudo-Dirac neutrinos can dominate the energy budget of the Universe and inject sizeable amount of entropy only at very late times, compared to the DM production from decay which occurs at temperature close to the mass scale of Σ (a possible exception is the case λ HΣ λ HΣ ). To a good approximation we can thus stick on a system of the form (43) and apply a posteriori possible entropy effects.
For simplicity we will describe two specific examples, namely all the pseudo-Dirac neutrinos lighter or heavier than Σ. In the first case all the source terms associated to the decays of the pseudo-Dirac neutrinos can be dropped. Moving to the quantities Y Σ,DM = n Σ,DM /s and x = m Σ /T as, respectively, dependent and independent variables, the system reduces to: where Y Σ,eq = 45 4g * π 4 x 2 K 2 (x) and H is the Hubble expansion rate H = 4π 3 g * . This last expression assumes that during the phase of DM generation the Universe is radiation dominated, this is reasonable since we have shown in the main text that the number density of the heavy neutrinos tends to dominate at rather low temperatures.
The numerical solution of this system has been presented in the left panel of Fig. 8 for some sample values of the relevant parameters.
The analytic expressions provided in the text correspond instead to suitable limits in which this set of equations can be solved analytically. In the regime λ HΣ ≥ λ HΣ the right-hand side of the equation for the DM is dominated, at early times, by the annihilation term and we have simply Y Σ = Y Σ,eq . In this regime we have only to solve the equation of the DM substituting Y Σ,eq on the right-hand side. The equation can be straightforwardly integrated for: K 2 (x) Y I , Y I = Y I,eq and we can again fix Y Σ = 0 on the right-hand side. Analytical solutions are reliable if the timescales of production and decay of Σ are well separated, otherwise one should refer to the numerical treatment.