FIMP realization of the scotogenic model

The scotogenic model is one of the simplest scenarios for physics beyond the Standard Model that can account for neutrino masses and dark matter at the TeV scale. It contains another scalar doublet and three additional singlet fermions (N_i), all odd under a Z_2 symmetry. In this paper, we examine the possibility that the dark matter candidate, N_1, does not reach thermal equilibrium in the early Universe so that it behaves as a Feebly Interacting Massive Particle (FIMP). In that case, it is found that the freeze-in production of dark matter is entirely dominated by the decays of the odd scalars. We compute the resulting dark matter abundance and study its dependence with the parameters of the model. The freeze-in mechanism is shown to be able to account for the observed relic density over a wide range of dark matter masses, from the keV to the TeV scale. In addition to freeze-in, the N_1 relic density receives a further contribution from the late decay of the next-to-lightest odd particle, which we also analyze. Finally, we consider the possibility that the dark matter particle is a WIMP but receives an extra contribution to its relic density from the decay of the FIMP (N_1). In this case, important signals at direct and indirect detection experiments are generally expected.


Introduction
The identification of the dark matter particle stands as one of the most pressing problems in fundamental physics today. Its solution requires physics beyond the Standard Model but it is not yet known what this new physics is. Most of the models studied in the literature assume that the dark matter consists of Weakly Interacting Massive Particles (WIMPs) and that its relic density is the result of a freeze-out process. One advantage of this scenario is that it naturally yields a relic density of the same order as the observed dark matter density -the so-called WIMP-miracle. In addition, WIMP models generally give rise to signals, in direct and indirect detection experiments as well as at colliders such as the LHC, that are within the reach of current experiments. Up to now, however, such signals have not been found and strong bounds on many of these models have been derived. If this situation persists for the next few years, the WIMP paradigm would likely have to be abandoned [1] and dark matter would have to be explained in some other way.
A simple and appealing alternative to the WIMP framework is provided by FIMP (Feebly Interacting Massive Particle) dark matter [2]. Its basic idea is that, in contrast to WIMPs, the dark matter interacts so weakly that it does not reach thermal equilibrium in the early Universe. Thus, its relic density is not the result of a freeze-out. Instead, the dark matter particles are slowly produced via decays or scatterings of the particles in the thermal plasma -a process dubbed freeze-in-but they are never abundant enough for their annihilations to be relevant. Consequently, the dark matter abundance steadily increases as the Universe cools down until the so-called freeze-in temperature is reached, and it remains constant afterward. Due to its feeble interactions, FIMPs do not give rise to observable signals neither in direct nor in indirect dark matter detection experiments. These experiments, therefore, provide an unambiguous way of testing or falsifying this scenario: if a signal were detected, one could immediately conclude that dark matter does not consist of FIMPs. Currently, FIMPs provide a viable and attractive framework to account for the dark matter.
Several explicit realizations of the FIMP framework have already been investigated [3,4,5,6,7,8,9,10,11,12,13]. FIMPs are necessarily singlets under the Standard Model (SM) gauge group so the two simplest extensions of the SM that incorporate a FIMP include a new singlet scalar [3,8,9] or a new singlet fermion [13], both of which give rise to an interesting phenomenology. In this paper, we will study a richer realization of the FIMP framework based on the scotogenic model (also known as the radiative seesaw model) [14]. This model is one of the simplest scenarios for physics beyond the SM that can simultaneously account for neutrino masses and dark matter at the TeV scale. It contains another scalar doublet and three additional singlet fermions (N i ), all odd under a Z 2 symmetry. Even though the phenomenology of this model has been extensively studied in a number of previous works -see e.g. [15,16,17,18,19,20,21,22,23]-, none of them considered the possibility of FIMP dark matter. The basic idea is that the couplings of one of the singlet fermions, N 1 , are so small that it does not reach thermal equilibrium in the early Universe and is instead produced via freeze-in. We show that dark matter production is dominated by the decays of the odd scalars and study the dependence of the resulting abundance with the different parameters of the model. In particular, the viable parameter space for FIMP dark matter is precisely determined and it is shown to span a wide range of dark matter masses, from the keV to the TeV scale. Besides freeze-in, the dark matter relic density receives an additional contribution from the so-called superWIMP mechanism [24] which strongly depends on the identity of the next-to-lightest odd particle. We identify an important region of the parameter space where this contribution is always negligible and freeze-in production is dominant.
Another interesting setup we discuss occurs when the dark matter particle (the lightest odd particle) is not N 1 but H 0 , a WIMP. In that case, N 1 , which is produced via freeze-in, decays into the dark matter, increasing its relic density and allowing for new viable regions in the parameter space. We show that interesting signals from direct and indirect detection experiments are generally expected in this configuration.
The rest of the paper is organized as follows. In the next section, we introduce the model and discuss the experimental bounds it is subject to. Then in section 3 we obtain the conditions necessary to ensure that N 1 does not reach thermal equilibrium in the early Universe. Our main results are presented in sections 4 and 5. In the former, we carefully study the production of FIMP dark matter and obtain the corresponding viable parameter space. Section 5 is dedicated to the case where the dark matter particle is H 0 and receives a contribution to its relic density from FIMP decays. Finally, we present our conclusions in section 6.

The model
The model we consider is the so-called scotogenic model [14], one of the simplest models that can simultaneously explain neutrino masses and dark matter at the TeV scale. In it, the SM is extended with a second Higgs doublet H 2 ≡ (H + , H 0 2 ) and three Majorana neutrinos N j (j = 1, 2, 3), all odd under an exact Z 2 symmetry (the SM fields are instead even under it). This symmetry forbids the coupling between H 2 and the quark fields, which would give rise to flavor changing neutral currents, and it guarantees the stability of the dark matter particle.
The Lagrangian of the model contains the following new terms involving the singlet fields The above expressions tell us that, if we want to generate light neutrino masses with new physics at the TeV scale, the product λ 5 y 2 2,3 must necessarily be very small (∼ 10 −11 ). This condition can be satisfied in different ways, however. One can, for instance, set λ 5 ∼ 10 −9 so that y 2,3 ∼ 0.1. Or one could fulfill it with λ 5 = 10 −1 and y 2,3 ∼ 10 −5 . Moreover, we can use equation (5) and (7) to set a lower bound on the Yukawa couplings: Smaller values of y 2,3 would fail to reproduce the observed neutrino mass scale. Since equation (3) has the same matrix structure as the usual seesaw equation, one can adapt the Casas-Ibarra parametrization [40] to it and express the Yukawa couplings in terms of the experimental data on neutrino masses and mixing angles. For the analogous case of two-right handed neutrinos which is relevant in our scenario, this procedure introduces only one free parameter [41,42], an angle that we take to be real. We assume a normal hierarchical spectrum for the neutrinos and took their oscillation parameters from [43]. In this way, we guarantee that all the models we consider in the following are compatible with current neutrino data.
The same interactions that generate neutrino masses induce lepton flavor violating processes such as µ → e γ and τ → µ γ at the 1-loop level [15,34]. Since these processes have not been observed, one must ensure that the predicted branching ratios are below the present experimental bounds. Given that the current limits read BR(µ → e γ) < 5.7 × 10 −13 [44] and BR(τ → µ γ) < 4.4 × 10 −8 [45], the former decay typically gives a stronger bound. In this model, the branching ratio for the µ → eγ process is [15] BR(µ → e γ) = 3α em where the loop function F 2 (x) varies in the range [3 × 10 −5 , 0.14] for x = [10 4 , 0.1] and we have already taken into account the fact that y 1 is negligible. Notice, in particular, that large Yukawa couplings, y 2,3 0.1, are strongly disfavored. In our analysis, we always impose that BR(µ → eγ), computed from equation (9), be below the experimental limit. Another important bound that must be taken into account is the dark matter constraint -the requirement that the predicted relic density agrees with the observed dark matter density. In this model there are two viable dark matter candidates, the lightest neutral scalar and the lightest singlet fermion, and the predicted relic density depends on how they were produced in the early Universe. While most previous works have assumed the usual freeze-out scenario, we want to examine the possibility that N 1 does not reach thermal equilibrium in the early Universe and is instead produced via freeze-in.

Out of equilibrium conditions
The basic requirement of the FIMP (or freeze-in) mechanism is that the dark matter particle does not reach thermal equilibrium in the early Universe. In the scotogenic model, only the fermions, which are gauge singlets, can play the role of FIMPs. Equilibrium can be prevented if their Yukawa interactions are sufficiently suppressed. In this section we analyze the different processes that can produce singlets and obtain the conditions necessary for them not to reach equilibrium. In particular we show that only one of them, denoted by N 1 , can play the role of a FIMP.
Because the singlet fermions have Yukawa interactions of the form N k L H 2 , they can be produced via the two-body decay of the odd scalars. The decay rate for the production of N 1 is approximately given by Γ(H 2 → N 1 L) = M H 2 y 2 1 /(8π) -see equations (16) and (17) below. Then, the out of equilibrium condition for this decay reads which for M H 2 ∼ 100 GeV implies If y 1 were larger than this value, N 1 would be produced abundantly enough to reach thermal equilibrium. This small value of the Yukawa coupling implies that, as already anticipated in the previous section, N 1 gives a negligible contribution to neutrino masses -see equation (8). The heavier singlets, N 2,3 , can be produced either via scalar decays (H 2 → N 2,3 L) or, if they are heavier than the scalars, via the inverse decay H 2 + L → N 2,3 , both of which are in equilibrium for y 2,3 10 −8 .
Since the bound from neutrino masses requires y 2,3 10 −6 , we can conclude that N 2 and N 3 necessarily reach thermal equilibrium in the early Universe. This model, therefore, admits only one FIMP: N 1 .
N 1 can also be produced via the decay of the heavier singlets or via 2 → 2 scatterings of SM leptons or odd scalars. All these processes are however subdominant and do not modify the equilibrium condition obtained above.
Contrary to our findings, it was stated in [17] that all three singlets could be out of equilibrium while explaining neutrino masses. The reason for this erroneous conclusion is that they failed to recognize the importance of the scalar decays as a production process for the singlets. Instead, they assumed that singlets were pairproduced via the annihilation of two leptons or two odd scalars. Since the rates of these processes depend on the neutrino Yukawa couplings to the fourth power (rather than the second), the out of equilibrium condition gives the wrong (and weaker) bound y i 10 −4 , which is consistent with the limit from neutrino masses. From our discussion, it should be clear though that the decays of the odd scalars (or the inverse decays mentioned above) cannot be neglected as they are the dominant production process for the singlets. Once these decays are taken into account it follows that only one singlet can be out of equilibrium.
Summarizing, the bound from neutrino masses implies that in this minimal setup (with three singlet fermions) only one FIMP is allowed. We will next show that this FIMP can easily account for the observed dark matter density via freeze-in.

FIMP dark matter
In this section we analyze the case where the dark matter candidate -that is the lightest odd particle-is the singlet that does not reach thermal equilibrium in the early Universe, which we denote by N 1 . That is, we consider the spectrum The N 1 relic density, Ω N 1 h 2 , will therefore have two contributions: one from the freeze-in mechanism and another one from the late decay of the next-to-lightest odd particle, which we call the superWIMP contribution [24]. We have therefore Since these two contributions are entirely independent -they become relevant at different temperatures and do not depend on the same parameters-we will study them separately.

The freeze-in contribution
Let us discuss first the freeze-in contribution to N 1 production. Since N 1 has a direct coupling to leptons and to odd scalars, its production will be dominated by the decays of the scalars (H 0 , A 0 , H ± ) while they are in equilibrium with the thermal bath. The N 1 yield, Y N 1 (T ) = n N 1 (T )/s(T ), is computed by solving the following Boltzmann equation [2] s T where s is the entropy density of the Universe, H(T ) is the expansion rate of the Universe at a given temperature and γ N 1 (T ) is the thermal averaged FIMP production rate. We have that where X = H 0 , A 0 , H ± and is a SM lepton. In this equation, K 1 (x) is the Bessel function of the second kind, and g X is the number of internal degrees of freedom of particle X. Specifically, g H 0 ,A 0 ,H + ,H − = 1. The decay rates that enter into this expression are calculated as where the approximations are valid unless there is a strong mass degeneracy between N 1 and one of the scalars. It is easy to verify that the decays of the heavier singlet fermions, N 2,3 → N 1¯ , give a negligible contribution to dark matter production. In fact, the corresponding decay rate is given by which is always much smaller than (16) and (17). Other negligible processes are the production of dark matter via scatterings of two Z 2 -odd particles or two SM particles. Both are always subdominant because the corresponding cross-sections are proportional to the fourth power of the Yukawa couplings. Thus, the N 1 abundance, Y N 1 , is solely determined by the Yukawa couplings (Y ν α1 ) and by the spectrum of odd scalar particles: m H 0 , m A 0 , m H ± . As we will see, for our purposes it is often a good approximation to consider all odd scalars to be degenerate, in which case we denote their common mass by m S .
From equations (14)(15)(16)(17), taking into account that s(T ) = 2π 2 g s T 3 /45, H(T ) = 1.66 √ g ρ T 2 /M P l and K 1 (x) ∼ 1/x for x 1, we have at high temperatures T > m S : Therefore, on the one hand we have that at T > m S the yield always scales as the square of the scalar masses and of the N 1 Yukawa couplings. On the other hand, at T m S the scalar particle abundance becomes Boltzmann suppressed and the production of dark matter is no longer efficient. As a result we have We have studied quantitatively the freeze-in production of dark matter in this scenario by solving numerically the Boltzmann equation (14) with the initial condition Y N 1 = 0 for T m S . Figure 1 shows the predicted dark matter abundance as a function of the temperature for different values of y 1 . The upper line corresponds to y 1 = 10 −8 and the lower one to y 1 = 10 −12 . In this figure the common scalar mass, m S , was set to 400 GeV. As stated before, the other parameters of the model are irrelevant. Notice, from the figure, that the abundance has the typical freeze-in behavior: it increases steadily until the so-called freeze-in temperature is reached, remaining constant afterward. Since the freeze-in temperature is determined by the mass of the decaying particle, it is the same for all the lines, as observed in the figure. Finally, the abundance is seen to depend quadratically on y 1 , as expected from equations (16), (17) and (20).
The dependence of Y N 1 on m S is illustrated in figure 2, which displays the dark matter abundance as a function of the temperature for different values of m S . In this figure y 1 was set equal to 10 −10 . One can clearly see that the freeze-in temperature increases with m S , with the result that the asymptotic value of Y N 1 decreases with m S . In fact, at low temperatures Y N 1 is about ten times smaller for m S = 2 TeV than for m S = 200 GeV. Notice from figures 1 and 2 that equation (20) is actually a very good approximation for the final yield obtained through the freeze-in mechanism.
In the previous two figures we have assumed a common mass, m S , for all the odd scalars. In general, however, there will be a mass splitting between the three  different states. To demonstrate that such mass splitting does not significantly affect our results, we show in figure 3 the dark matter abundance as a function of the temperature for different mass splittings. Notice that the variation in the final abundance due to the different kind of spectra is indeed very small. It is, therefore, a very good approximation to compute the dark matter abundance assuming that all odd scalars have the same mass m S . The relic density of dark matter, Ω N 1 h 2 , is related to the asymptotic value of Y N 1 at low temperatures by where T 0 = 2.752 K is the present day CMB temperature. It is this quantity that should be compared with the observed dark matter density as measured by WMAP  [46] and PLANCK [47]. For dark matter production via the freeze-in mechanism, the N 1 relic abundance can be estimated as where we used equations (20) and (21). Notice that this expression has the expected dependence on m S , y 1 and M 1 . Figure 4 displays the N 1 relic density as a function of M 1 for m S = 400 GeV and different values of y 1 . For M 1 we considered a minimum value of 1 keV as indicated by phase space density analysis [48,49] and by the requirement of cold or warm dark matter. The maximum value was taken to be 100 GeV in agreement with the idea that all odd particles live at or below the TeV scale. The horizontal band shows the region that is compatible with current observations. Notice that as we increase the  mass a smaller value of y 1 is needed to be consistent with the data. Hence, whereas a keV particle requires y 1 ∼ 10 −8 a 100 GeV particle requires y 1 ∼ 10 −12 .
The viable parameter space for freeze-in dark matter in the scotogenic model is shown in figure 5. It displays, in the plane (M 1 ,y 1 ), the regions that are consistent with the observed dark matter density for different values of m S . The freeze-in mechanism is thus able to explain the dark matter over a wide range of masses, from the keV to the TeV scale. Notice that at a given dark matter mass, the heavier m S the larger y 1 . This figure is one of our main results, as it indicates the regions in the parameter space of the scotogenic model where the observed dark matter density can be accounted for entirely via freeze-in.
If N 1 is very light, M 1 ∼ 1 − 10 keV, the resulting dark matter is warm rather than cold, with important implications for structure formation in the early Universe. Without freeze-in it is not possible to obtain warm dark matter in the scotogenic model because N 1 would thermalize and later decouple while relativistic, yielding  a relic density about three orders of magnitude larger than observed. To make such scenario compatible with current observations would require either entropy dilution, e.g. via the decay of some other particle, after N 1 production [50] or a nonthermal production mechanism within a low reheating temperature scenario [50], both entailing significant departures from the model. Freeze-in provides instead a natural and simple way of obtaining warm dark matter in the scotogenic model. If, in addition to freeze-in, other mechanisms contribute to dark matter production, the lines in figure 5 provide an upper bound on the coupling y 1 at a given value of M 1 and m S . As mentioned at the beginning of this section, in the scotogenic model the relic density of N 1 receives also a superWIMP contribution from the decays of the next-to-lightest odd particle after it has frozen out. Let us now turn our attention to that contribution.

The superWIMP contribution
In the superWIMP mechanism, the contribution to the dark matter relic density from the late decay of the next-to-lightest odd particle (NLOP from now on) is given by where Ω f reeze−out N LOP h 2 is the relic abundance, obtained via the usual freeze-out mechanism, of the NLOP . In the scotogenic model, there are essentially two possibilities for the NLOP: N 2 or one of the scalars. Next, we will in turn consider these two options.

N 2 as the NLOP
If N 2 is the NLOP it will decay into dark matter via the scalar-mediated threebody process N 2 → N 1 . The requirement that this decay happens after the N 2 freeze-out (at T ∼ M 2 /20) implies that This condition yields an upper bound on the product y 1 y 2 : which is always satisfied in this scenario -the out-of-equilibrium condition gives a stronger bound. On the other hand, the lifetime of N 2 should be smaller than about 1 second in order to not affect the Big Bang Nucleosynthesis (BBN) epoch. This requirement implies a lower bound on the product of the Yukawa couplings, namely At high values of the dark matter mass, this condition is very restrictive. If, for instance, M 1 ∼ 100 GeV and m S , M 2 ∼ 1 TeV, it is not possible to satisfy it as we know, from figure 5, that y 1 should be no larger than about 10 −12 (to avoid dark matter overproduction) and that y 2 cannot be of order 1 due to the µ → eγ bound. For M 1 ∼ 1 keV and the same values of m S and M 2 , y 1 should be smaller than about 10 −8 and the above bound is satisfied for y 2 10 −4 . Taking y 2 ∼ 10 −2 as the upper limit on y 2 allowed by µ → eγ, the BBN constraint would exclude models with y 1 10 −10 or equivalently with M 1 100 MeV. We can also use equation (26) to set a lower bound on the mass of N 2 . Since m S 100 GeV, y 1 10 −8 and y 2 10 −1 -10 −2 , we get M 2 10 GeV. Thus, the FIMP mechanism combined with the BBN constraint above tells us that M 2 and M 3 necessarily lie around the electroweak scale.
Regarding the value of Ω f reeze−out N 2 h 2 , previous studies have already shown that N 2 -N 2 annihilations are not very efficient and usually require, to be consistent with the observed dark matter density, values of the Yukawas couplings so large that they run into conflict with the bounds from µ → eγ. Coannihilations between N 2 and the scalars significantly help to increase the total annihilation rate, reducing the relic density and alleviating the tension with the µ → eγ bound. This situation is illustrated in figure 6, which displays a scatter plot of the N 2 relic density versus M 2 .
In it we have randomly varied all the parameters of the scotogenic model over a wide range: 1 keV ≤ M 1 ≤ 100 GeV, 100 GeV ≤ M 2 ≤ 1 TeV, 1 TeV ≤ M 3 ≤ 3 TeV, M 2 ≤ m H i ≤ 3 TeV, 10 −12 ≤ |Y ν α1 | ≤ 10 −8 , 10 −3 ≤ λ L ≤ 1. All points in this figure satisfy the constraints from neutrino masses, µ → eγ, and BBN. To precisely compute the relic density we used micrOMEGAs [51], which automatically includes all the relevant processes and takes care of possible resonant or coannihilation effects. With the goal of isolating the effect of coannihilations, we have divided the sample into two sets according to the mass splitting between N 2 and the scalars. The masssplitting is small for the red points (allowing for coannihilations) and large for the blue points (excluding coannihilation effects). The horizontal band corresponds to the observed dark matter density. Notice that coannihilations are essential to obtain a relic density in agreement with the observations. If m H 0 > 1.5 M 2 (blue points), the N 2 relic density after freeze-out is always very large -at least four orders of magnitude larger than the observed dark matter density. Thus, compatibility with current data requires M 1 /M 2 10 −4 , according to equation (23). And since M 2 is at most of order TeV, M 1 necessarily lies below the GeV scale. A large hierarchy between M 1 and M 2 is thus an essential condition in this scenario. If, on the contrary, M 1 < m H 0 ≤ 1.5 M 2 (red points), the N 2 relic density can even reach values below the observations. Consequently, no strict bounds on M 1 /M 2 can be derived based on the relic density.
If the dark matter density were dominated by the superWIMP contribution, h 2 , one would need to ensure that the dark matter is non-relativistic at the onset of structure formation; otherwise it would behave as hot dark matter. This condition leads to a relation between the N 2 decay time and the ratio M 1 /M 2 . A detailed analysis of this issue can be found in [17]. They found, in particular, that one can obtain warm dark matter for 24 keV M 1 24 MeV(M 2 /100 GeV). In figure 6 we have also displayed, for three different values of M 1 (100 GeV, 100 MeV, 100 keV), the regions where the superWIMP contribution accounts for the entire dark matter density. If M 1 = 100 MeV, for example, then along the dashed line the superWIMP contribution agrees with the observed relic density. Models above that line are excluded (for that value of M 1 ) as they overproduce dark matter whereas models below that line require the freeze-in contribution to be compatible with the data. Even though the relic density constraint is satisfied along the dashed-dotted line for M 1 = 100 GeV, that value of M 1 is actually ruled out by the BBN bound, as explained before. If M 1 = 100 keV, the superWIMP contribution accounts for the dark matter along the solid line and one obtains warm dark matter.

A scalar as the NLOP
If one of the scalars is the NLOP, its direct decay into N 1 and SM leptons after decoupling from the thermal bath will give an additional contribution to the dark matter abundance. The condition that the decay takes place after the scalar freezeout but before BBN translates into where we have implicitly assumed that the decaying scalar and N 1 are not highly degenerate. These bounds are easily satisfied for the range of parameters relevant for freeze-in -see figure 5. For concreteness, in the following we assume the NLOP scalar to be H 0 but it must be kept in mind that the results for the other scalars are similar. The relic density of H 0 in this scenario is very much alike that in the inert-doublet model. A remarkable feature of this model is that if M W < m H 0 500 GeV the relic density is always too small to satisfy the dark matter constraint. The reason being that the annihilation into gauge bosons are so efficient that they deplete the dark matter abundance well below the observed value. Only for masses above 500 GeV (or below M W ) it is possible to satisfy the dark matter bound. Figure 7 shows a scatter plot of the H 0 relic density versus m H 0 obtained after varying all the parameters of the scotogenic model (m H i ≤ M 2 ≤ 1 TeV, M 2 ≤ M 3 ≤ 3 TeV, 100 GeV ≤ m H i ≤ 1 TeV and the others as before) and selecting those consistent with neutrino masses, µ → eγ, and BBN. As before, the horizontal band shows the observed dark matter density. Notice from the figure that the relic density increases with the mass and that, as expected, it only crosses the experimental value for masses above 500 GeV or so. This fact has a very important implication: if m H 0 < 500 GeV the superWIMP contribution to the relic density is negligible and the entire dark matter density has to be explained via the freeze-in mechanism. That is, in contrast to the case where N 2 is the NLOP, we can identify an important region of the parameter space, m H 0 < 500 GeV, where the freeze-in contribution is always dominant.
If m H 0 > 500 GeV, the superWIMP contribution could be the dominant one. In that case, since the H 0 relic density is never much larger than the observed dark matter density, a mild hierarchy between N 1 and H 0 is required, m H 0 /M N 1 4 (see figure 7). Even for such large values of m H 0 , however, the freeze-in contribution can dominate the N 1 relic density.

Implications
As we have seen, FIMP dark matter can indeed be realized in the scotogenic model. In general, the relic density of N 1 is the sum of a freeze-in contribution and a contribution from the decay of the NLOP. Whether one or the other dominates depends strongly on the parameters of the model. Let us now briefly discuss the implications of this scenario.
A generic prediction of FIMP models is that the NLOP, which must decay into the FIMP, is very long-lived [2], providing a possible way of testing these scenarios at colliders such as the LHC. In the scotogenic model, the most interesting signal occurs when the charged scalar is the NLOP. In that case the relic density is expected to be dominated by the freeze-in process and, from equation (22), we have that Now, let us suppose that this charged scalar, with a mass in the range [100 GeV, 1 TeV], is produced at the LHC. Its decay width is given by equation (17). Therefore, taking into account the value of y 1 derived above we obtain Thus, the H + decay length, l(H + ), is (ignoring for the moment the Lorentz boost factor) 3 meters 1 TeV Including the Lorentz boost factor amounts to multiplying this upper limit by a factor from 2 to 7. Thus, for dark matter masses in the range [10 keV, 1 MeV] the decay length is below 10 meters and H + decays inside the detector, leaving a charged lepton plus missing energy signature that could be searched for at the LHC. If the decay happens instead outside the detector, evidence for H + could be found at the LHC via searches for long-lived charged particles [52]. It is beyond the scope of the present paper, however, to determine whether these signals can actually be used to set meaningful constraints on this scenario. Another generic feature of FIMP dark matter is the absence of signals at direct or indirect detection experiments -a direct consequence of the feeble interactions that are required to prevent the dark matter from reaching thermal equilibrium in the early Universe. These experiments provide, nonetheless, an unambigous way of falsifying this scenario: as soon as a positive signal is confirmed in any dark matter detection experiment we would learn that dark matter does not consist of FIMPs and more specifically that the scenario we studied in this section is ruled out. Such signal would instead give a strong support to the WIMP paradigm of dark matter. But if the next generation of dark matter experiments, such as XENON1T [53], fails to find evidence of dark matter, the WIMP framework would be in trouble and alternative scenarios that can naturally explain the absence of such evidence would become much more appealing. In that hypothetical future the FIMP scenario could become the standard framework to account for the dark matter. Only time will tell which of these two possible outcomes regarding dark matter detection will actually be realized.

FIMP decay into dark matter
In the previous section we assumed that the singlet fermion that does not reach thermal equilibrium in the early Universe (N 1 ) was also the lightest particle odd under the Z 2 symmetry, and consequently the dark matter candidate. It may well be though that N 1 is not the lightest odd particle so that the dark matter candidate is instead one of the neutral scalars or another singlet fermion. In that case, N 1 is unstable and decays into the dark matter, increasing its relic density. Thus, N 1 modifies the regions where the dark matter constraint is satisfied, allowing for regions which in the usual freeze-out scenario are under-dense (Ω f reeze−out h 2 < 0.1) to become compatible with the observed dark matter density. Since the singlet (say N 2 ) relic density obtained via freeze-out is typically larger than the observed one, see e.g. figure 6, an additional contribution from FIMP decay is usually not welcome as it will only help in very specific cases. Much more interesting is the situation where one of the neutral scalars is the dark matter candidate, for we know that over a significant region of the parameter space its freeze-out relic density is very small -see e.g. figure 7. For definiteness, we take H 0 as the dark matter particle and assume that all the odd scalars are lighter than N 1 , m H 0 < m A 0 , m H ± < M 1 . Notice that, contrary to the discussion in the previous section, the dark matter particle in this case is a WIMP.
The H 0 relic density will receive two contributions, one from freeze-out and one from the late decay of N 1 . We can then write with Let us now proceed to calculate Ω f reeze−in N 1 h 2 in this case. The dominant freeze-in production process is the inverse decay of N 1 , X + → N 1 , where X denotes an odd scalar and is a SM lepton. The N 1 yield, Y N 1 (T ) = n N 1 (T )/s(T ), is computed by solving the same Boltzmann equation as in the previous section, equation (14), but with a different production rate where g N 1 = 2 because N 1 is a Majorana fermion. The decay width for the three decay channels of N 1 into scalars are given by Hence, the total decay rate of N 1 is where the last approximation is valid unless N 1 is highly degenerate with the scalars.
The abundance Y N 1 at certain temperature T can then be expressed as whereas the N 1 relic density is Finally, we can approximate Ω N 1 −decay Thus, a coupling of order 10 −12 is required to account for the entire dark matter density via the decay of N 1 . The above result holds provided that N 1 decays after the H 0 freeze-out, H 0 ). Since the N 1 decay rate is given by whereas one can see that this condition is easily satisfied. In order to not alter the predictions of BBN, one must also ensure that Γ N 1 1/0.3 sec −1 = 2.2 × 10 −24 GeV, which is seen to be fulfilled for the values required to obtain the correct dark matter density.
The idea then is that if Ω f reeze−out H 0 h 2 < Ω DM h 2 we can always choose a value of y 1 such that the contribution from the decay of N 1 compensates for the deficit and one gets a relic density in agreement with the observations, Ω H 0 = Ω DM . That is, the presence of the FIMP allows us to enlarge the viable parameter space of the model by rescuing those regions where freeze-out gives a too small relic density. In particular, the region m H 0 500 GeV becomes viable within this setup.
The resulting scenario is quite similar to that discussed in [37]. The difference being the mechanism that allows to increase the relic density. In [37] it was the m H 0 (GeV)  The solid line shows the current bound by the LUX experiment [54] whereas the dashed line displays the expected sensitivity of XENON1T [53].
coannihilations with the singlet fermions whereas in our case is the late decay of the FIMP.
Since the dark matter particle H 0 is a WIMP, the usual direct and indirect detection signals are expected and one must make sure that current bounds are respected. The dark matter phenomenology of H 0 is reminiscent of that in the inert doublet model. Direct detection, for instance, proceeds via a Higgs mediated diagram and is determined by the coupling λ L = (λ 3 + λ 4 + λ 5 )/2. Figure 8 shows a scatter plot of the spin-independent direct detection cross section versus the dark matter mass. The figure was obtained after randomly varying the different parameters of the model (1 TeV ≤ M 1 ≤ 3 TeV, M 1 ≤ M 2 ≤ 3 TeV, M 2 ≤ M 3 ≤ 3 TeV, 100 GeV ≤ m H i ≤ 1 TeV and the others as before) and imposing the known experimental bounds (neutrino masses, µ → eγ, etc.). The H 0 relic density is consistent with the observed dark matter density thanks to the contribution from N 1 decays. For comparison we show the current experimental bound (solid line) [54] and the m H 0 (GeV) expected sensitivity of future experiments (dashed line) [53]. Even though several models are already excluded (those above the solid line) and many more will be probed by future experiments (those above the dashed line), one can still find models with small values of σ SI over the entire range of masses we explore. Direct detection bounds therefore do not restrict the range of the dark matter mass in this scenario.
Unsurprisingly, the indirect detection bounds turn out to be more constraining. Indeed, since Ω f reeze−out H 0 h 2 Ω DM h 2 , we expect annihilation rates larger than those typically associated with WIMPs, σv H 0 σv thermal ∼ 3 × 10 −26 cm −3 /s. Figure  9 shows a scatter plot of σv versus the dark matter mass. As before, the correct relic density is obtained via N 1 decays and the experimental bounds were taken into account. In blue we show the models that are excluded by the direct detection bound on σ SI -see figure 9. In red we show instead the points that are consistent with that bound. Notice that σv can indeed be much larger than the so-called thermal value.
The solid and dashed lines show the current bounds obtained by the Fermi-LAT collaboration, for dark matter annihilation into b quarks [55] and W gauge bosons [56], respectively. They exclude all models with m H 0 300 GeV. For higher values of the dark matter mass, m H 0 300 GeV, one can easily find models compatible with both direct and indirect detection constraints.
In contrast to the scenario with FIMP dark matter discussed in the previous section, this setup, where H 0 is the dark matter particle and the decay of N 1 contributes to its relic density, can be probed by both direct and indirect detection experiments. And as we have seen, the expected signals are generally significant.

Conclusions
We have shown that in the scotogenic model -one of the simplest extensions of the SM that can account for neutrino masses and dark matter at the TeV scale-one (and only one) of the singlet fermions, N 1 , can be out of equilibrium in the early Universe and behave as a FIMP, with important implications for the phenomenology of this model. This setup predicts, for instance, that one of the light neutrinos is essentially massless. Within this framework the dark matter candidate can be a FIMP, N 1 , or a WIMP, H 0 . In the former case, the relic density of dark matter receives two contributions, one from freeze-in and another one from the late decay of the next-to-lightest odd particle -the superWIMP contribution. The freeze-in contribution was found to be dominated by the decays of the scalars and its dependence with the different parameters of the model was examined in detail. Specifically, we determined the regions in the plane (M 1 , y 1 ) where freeze-in can account for the observed dark matter density and found that they span a wide range of masses, from the keV to the TeV scale. The superWIMP contribution was also discussed and shown to strongly depend on the identity of the next-to-lightest odd particle. In the latter case, when H 0 is the dark matter particle, the relic density is not only the result of a freeze-out but receives and additional contribution from the late decays of N 1 . This second contribution allows to increase the dark matter relic density, opening up new viable regions in the parameter space of the model. Thanks to this contribution from N 1 decay, regions that within the standard scenario feature a too small relic density, such as m H 0 500 GeV, can become compatible with the observed dark matter density. We demonstrated that in this case one generally expects observable signals at direct and indirect dark matter experiments.