Neutrino dark matter and the Higgs portal: improved freeze-in analysis

Sterile neutrinos are one of the leading dark matter candidates. Their masses may originate from a vacuum expectation value of a scalar field. If the sterile neutrino couplings are very small and their direct coupling to the inflaton is forbidden by the lepton number symmetry, the leading dark matter production mechanism is the freeze-in scenario. We study this possibility in the neutrino mass range up to 1 GeV, taking into account relativistic production rates based on the Bose-Einstein statistics, thermal masses and phase transition effects. The specifics of the production mechanism and the dominant mode depend on the relation between the scalar and sterile neutrino masses as well as on whether or not the scalar is thermalized. We find that the observed dark matter abundance can be produced in all of the cases considered. We also revisit the freeze-in production of a Higgs portal scalar, pointing out the importance of a fusion mode, as well as the thermalization constraints.


Introduction
One of the outstanding mysteries of the Universe is the nature of dark matter (DM). An attractive minimal option is provided by sterile neutrinos, whose existence is strongly suggested by the observed neutrino masses. The smallness of the latter can elegantly be explained by the seesaw mechanism [1,2,3,4,5,6]. When the active-sterile neutrino mixing is sufficiently small, the lightest sterile neutrino can be very long-lived and play the role of DM. In the simplest scenario proposed by Dodelson and Widrow [7], such neutrinos can be produced via mixing with the active neutrinos in a thermal bath of the Standard Model (SM) particles, although the sterile neutrinos do not thermalize themselves. This minimal option now appears to be in conflict with a number of observations [8,9,10,11,12,13,14,15,16,17] (see e.g. [18] for a review).
An interesting possibility, which we explore in detail, is to generate the observed relic abundance of the sterile neutrino DM through the freeze-in mechanism [41,42]. This scenario requires a tiny coupling and a negligible initial DM abundance. The correct relic density is then gradually built up via this feeble coupling along the evolution of the Universe. A successful realization of this mechanism involves an extra scalar field, whose decay produces sterile neutrinos [29,30]. The vacuum expectation value (VEV) of such a field can be responsible for the Majorana neutrino masses [43,44]. For a small enough coupling, the freeze-in mechanism is at work and the correct DM density can be produced.
In our work, we perform an in-depth analysis of the freeze-in production of sterile neutrinos in the mass range up to 1 GeV, taking into account different possible production regimes, relativistic reaction rates with the Bose-Einstein distribution function, thermal masses and main effects of the phase transitions. Previous studies have mainly focused on the keV mass range [29,30,45,34,46], in which case the active-sterile mixing angle is required to be below 10 −5 or so. In our case, the requisite mixing must be much smaller calling for a symmetry justification. The required symmetry can be identified with the neutrino parity which acts on the lightest sterile neutrino only.
We find that the neutrino production often takes place in the relativistic regime, in which case the Bose-Einstein distribution should be used for the initial state scalars. This differs from many previous studies which have resorted to the non-relativistic Maxwell-Boltzmann approximation. In order to take quantum statistics into account, we follow the approach of [47], [48] and extend it to asymmetric reactions. The resulting rate enhancement depends strongly on the thermal masses and ranges from O(1) to two orders of magnitude in the vicinity of the 2d order phase transition.
We also take into account the most important effects of the phase transitions. First of all, these affect the presence or absence of certain couplings which depend on scalar VEVs. In addition, the mass change at the phase transition can facilitate particle production. The DM production modes depend on whether or not the scalar is thermalized. It couples to the SM via the Higgs portal [49,50,51]. Then, its thermalization depends on the Higgs portal coupling and the maximal temperature. To this end, we revisit the scalar production through the Higgs portal couplings and the consequent thermalization constraints. We find, in particular, that the 2 → 1 reaction (fusion) plays an important role in this analysis.
In this work, we are mainly interested in reproducing the correct DM relic abundance. To this end, we solve the Boltzmann equation for the number density rather than the momentum distribution function (unlike e.g. [34]). We reserve the latter for future work.
The paper is organized as follows. In Section 2, we introduce our model, discuss the leading thermal corrections and compile the current constraints on sterile neutrino DM. In Section 3, we generalize the relativistic reaction rates of [47], [48] to asymmetric reactions. Thermalization constraints are discussed in Section 4. Our main results are presented in Sections 5 and 6. We conclude in Section 7.

The model
In this work, we focus on a simple set-up: the SM is extended by a real singlet S and some number of right-handed (sterile) neutrinos ν R i . 5 The lightest of them is assumed to constitute long-lived dark matter.
We assume that the Majorana masses are produced entirely by the singlet scalar VEV. This can be implemented through a lepton-number discrete symmetry forbidding the bare mass: The relevant Lagrangian terms are then The above symmetry results in 2 useful properties: • no inflaton coupling to ν R i ν R j is allowed (assuming that the inflaton carries no lepton charge). Otherwise, inflaton decay would normally dominate the neutrino production.
• diagonalizing the neutrino mass matrix diagonalizes the S coupling to neutrinos (neglecting the Dirac contributions). As a result, there is no flavor change and a heavier ν cannot produce a lighter ν via its decay with S-emission. Thus, we can focus on direct freeze-in production of the lightest ν.
Let us denote the lightest mostly-sterile neutrino ν and its coupling to S λ: 5 Their number can be significantly larger than 3 as motivated by string theory [52]. Here, a sterile neutrino is defined as a fermion that has a Yukawa coupling with the SM neutrinos as well as a Majorana mass term.
Its mass is then M = λ S neglecting the Dirac mass contribution. Throughout this paper we assume that the relevant Yukawa couplings are very small such that the usual Higgs decay does not produce a significant amount of dark matter. The resulting active-sterile mixing angle Θ ∼ y H /λ S 10 −5 is also very small.

The scalar sector
The scalar sector of the model includes the Higgs field H and the real scalar S. The potential invariant under the S → −S symmetry is given by Here we use the unitary gauge H = (0, h/ √ 2) T . Both H and S must develop non-zero VEVs v and u, respectively. These are given by The mass matrix at this point is Since the couplings are real and we require v 2 > 0, u 2 > 0, the mass matrix M 2 is positive definite if and only if M 2 can be diagonalized by the orthogonal transformation where O = cos θ sin θ − sin θ cos θ (9) and the angle θ satisfies The mass squared eigenvalues are given by The above equation implies sign(m 2 2 − m 2 1 ) = sign(cos 2θ) sign(λ s u 2 − λ h v 2 ). We will primarily be interested in the small mixing case. (E.g., for a light singlet, the meson decay and LEP constraints on the mixing angle are very strong.) Thus, it is convenient to employ the small angle approximation sin θ 1 and neglect the θ 2 -terms. In this case, the eigenvalues can be relabelled according to the state composition and satisfy The mixing angle can then be expressed as This form is convenient since stability considerations bound the first factor by 1. Clearly, for m s close to m h our approximation fails. When m h and m s are substantially different, the mixing angle is bounded by |θ| < m s /m h , m h /m s . In most of our parameter space, the mixing angle is indeed very small.
In what follows, the sign of θ is unimportant, so will denote by θ the magnitude of the mixing angle.

Thermal corrections
At high temperature, the scalar potential gets modified by the thermal corrections. The main effect is captured by the thermal masses which amount to the replacements where c h 3 16 g 2 + 1 16 Here g 1,2 are the SM gauge couplings and y t is the top-quark Yukawa coupling. At high T , the minimum of the potential is at v = u = 0. The transition to non-zero VEVs takes place at the critical temperatures: v = 0 → v = 0 at T v c and u = 0 → u = 0 at T u c . The dynamics of the transition is somewhat complicated and proceeds in steps: at the first stage, one of the VEVs stays zero and the other becomes non-zero, while at the second stage both of the fields attain non-zero VEVs. On the other hand, we find that the neutrino DM production depends on the critical temperature rather weakly. Therefore, it suffices for our purposes to approximate the critical temperature by that of the first stage, T 2 c = |µ 2 i |/c i . The critical temperatures can then be expressed in terms of the physical masses and the couplings for sin θ 1: It is important to include the thermal masses (14) in the calculation of the reaction rates. This is dictated by their correct high temperature behaviour.
The neutrino thermal masses, on the other hand, can be neglected since these are suppressed by λ 2 .

Constraints on sterile neutrino dark matter
In Fig. 1, we collect the most stringent limits on the active-sterile mixing sin 2 Θ as a function of the sterile neutrino mass M . We assume the abundance of sterile neutrinos to be equal to the dark matter density measured by Planck [53,54]. Since our dark matter candidate decays into active neutrinos and other SM states, there are various strong constraints on this scenario. The most relevant ν decay modes are [55,56,57,18]: where ν a indicates an active neutrino. For a heavier ν, further decay modes become relevant, e.g. those involving muons. Here we are assuming that the mixing with the electron neutrino dominates.
In the dark grey region, the sterile neutrino lifetime is shorter than the age of the Universe. Sterile neutrinos are always produced in a thermal bath via the sterile-active mixing. This leads to the "overproduction" constraint indicated by the dashed purple line, above which the sterile neutrino abundance exceeds that of dark matter.
The sterile neutrino radiative decay is particularly relevant for X-ray and gamma-ray line searches. For sterile neutrino masses M 50 keV, searches of decaying dark matter signals have been carried out using a wide range of X-ray telescopes like XMM-Newton [8,58], Suzaku [59], HEAO-1 [8], INTEGRAL [14,13], Swift [60] and CHANDRA [61,62]. We collect most of them in the dark blue shaded area. 6 Among them, the CHANDRA satellite provides the strongest limits [17]. Most recent bounds from the X-ray microcalorimeter NuSTAR [64], looking at the Galactic Bulge, are displayed in dark cyan. The limits from searches for sterile neutrino decay lines using the Gamma-ray Burst Monitor onboard the Fermi Gamma-Ray Space Telescope (Fermi-GBM) [65] are shown in red. The green region is further constrained by INTEGRAL [13] searching for spectral lines from dark matter with a mass up to 14 MeV, decaying in the Milky Way halo. Gamma-ray lines searches further constrain our sterile neutrino dark matter parameter space at higher masses: we show the bounds from COMPTEL [66,57] (magenta), EGRET [67,57] (orange) and Fermi Large Area Telescope (Fermi-LAT) [16] (red).
Finally, measurements of the cosmic microwave background (CMB) allow us to constrain sterile neutrino decays leading to early energy injections [68,69,70,71,72]. The relevant decay modes are ν a e + e − and to ν a γ. Using the bounds on the corresponding decay rates from Ref. [73] with appropriate photon flux rescaling, we obtain the CMB bounds shown by the yellow dashed lines.
We see that the resulting constraints on the mixing angle are very strong. For example, for M close to 1 GeV, the bound on Θ is of order 10 −18 . Such small values appear unnatural. Within our simple model, the angles are input parameters, while in various extensions their small values can be justified by flavor-dependent symmetry. Indeed, in addition to the lepton number Z 4 , one may impose a Z 2 symmetry which acts on the lightest Majorana neutrino ν R 1 : This forbids the corresponding Yukawa couplings and sets Θ = 0. Assuming that this Z 2 is broken at some scale, the effective Yukawa couplings can be generated by higher dimensional operators. As a result, very small mixing angles can be generated. Since in the limit Θ → 0 the system becomes more symmetric, small mixing angles are natural according to the t'Hooft criterion [74]. Longevity of the lightest sterile neutrino can be achieved at small masses and/or small mixings. While most research efforts have focused on the first option, here we are considering the second possibility in more detail. We also see that, given the vast (Θ, M ) parameter space, dark matter decay may be observed, e.g. via monochromatic X-or gamma rays.

Relativistic rates for asymmetric reactions
Neutrino dark matter can be produced through a number of reactions. These include both scattering and decay which take place in the relativistic regime, i.e. when the temperature exceeds the particle masses. Since there are bosons in the initial state, relativistic Bose-Einstein enhancement can be very significant and the reaction rates must take it into account. The relevant results for symmetric reactions, that is, involving particles with the same mass in the initial state, have been obtained in [47], [48]. In our case, some reactions can be asymmetric, e.g. H + S → X, and these results must be generalized to particles of different masses.
In this section, we generalize the relativistic reaction rates based on the Bose-Einstein statistics [47], [48] to processes involving particles with different masses. The a → b reaction rate per unit volume is given by the general expression Here M a→b is the QFT transition amplitude, in which we also absorb the initial and final state symmetry factors, and f (p) is the momentum distribution function. For the freeze-in scenario, the density of the final state particles is small so that the enhancement factors 1 + f (p j ) can be set to one. On the other hand, it is important to keep the full Bose-Einstein distribution functions f (p i ) for the initial state and their replacement by the Maxwell-Boltzmann ones can lead to a rate underestimate by orders of magnitude. We are particularly interested in the 2 → 2 reactions. The reaction rate can be expressed in terms of the cross-section, with The cross section is defined by where M is the QFT transition amplitude. Here we absorb the symmetry factors for the initial state directly into σ(p 1 , p 2 ). The calculation is most easily performed in the center-of-mass (CM) frame, so let us convert the integral into that frame. The CM frame for each pair p 1 , p 2 is the frame where p 1 + p 2 has only zero spacial components. Let us consider how the integration measure transforms as we go to the CM frame. The Lorentz invariant measure is such that Any time-like vector p can be Lorentz-transformed to the form with the explicit parametrization in terms of rapidity η and angular coordinates θ, φ being In other words, in the convention p = (p 0 , p 3 , p 2 , p 1 ) T , we have The p-vector in the form (E, 0, 0, 0) T corresponds to the CM frame and E > 0 is half the CM energy. The p-integration measure becomes where Ω p is the solid angle in p-space. Now apply the same Lorentz transformation Λ(p) to the vector k, where we have used the fact that Λ(p) is a constant Lorentz transform with respect to the variable k so that the measure remains invariant. We drop the prime for convenience, remembering that k now is in the CM frame. Ω k denotes the corresponding solid angle in that frame.
Let us now integrate the delta functions. We can explicitly integrate over k 0 and |k|. In the CM frame, the delta functions become This enforces We then have which allows us to rewrite the integration measure as where in the integrand one must set k 0 and |k| to their values given by Eq. (32). Note that E is half the CM energy.
Since the cross section in the CM frame is a function of E only, the angular dependence comes entirely from the distribution functions. We have Here we have used k 3 = |k| cos θ k . Integration over Ω p gives 4π and the integral over φ k is 2π. Let us now integrate over θ k . The integral can be reduced to (This result can most easily be obtained by the change of variables to y = e bx .) We thus get with k 0 and |k| given by (32). This expression reduces to that of [47], [48] for equal masses, It is important to note that the masses here must include thermal corrections (14). This is necessary for the correct high temperature behaviour: only when a thermal correction to m has been included. The latter also regularizes the infrared divergence in the massless limit.

Thermalization constraints
In this work, we focus on freeze-in production of sterile neutrinos. Freeze-in calculations are reliable only if the produced particles do not thermalize. This requires the coupling between the thermal bath and the frozen-in particles as well as self-interaction of the latter to be sufficiently small. In this section, we delineate parameter space consistent with these conditions. We use relativistic formulas for the reaction rates, taking into account quantum statistics for the initial state.
Let us consider the regime where H and S develop VEVs v and u, respectively. We can parametrize them in the unitary gauge as In terms of the 4-component Majorana neutrino ν, the relevant to our calculation terms in the Lagrangian are where M = λu and we have neglected the scalar mixing. The Majorana notation has the advantage that the νν final state includes all combinations of 2-component neutrinos and antineutrinos.

Sterile neutrino thermalization
We show in Sec. 5 that the main production channel for sterile neutrinos is the decay s → νν.
Here we assume that s is in thermal equilibrium and m s M . For a sufficiently large λ, the decay is efficient and the neutrino number density n ν (T ) approaches its equilibrium value at a given temperature, n eq ν (T ). In this case, the reverse process νν → s becomes important and the neutrinos tend to equilibrate with the thermal bath of s. Thus, we use the non-thermalization criterion n ν (T ) < n eq ν (T ) (40) for any T down to temperatures around M/3. (At lower T , n eq ν (T ) is exponentially suppressed.) The number density n ν is calculated via the Boltzmann equatioṅ where H is the Hubble rate, with g * being the number of active SM degrees of freedom. Γ 12 (s → νν) is the reaction rate per unit volume (see Section 5 for an explicit expression). It is calculated with the Bose-Einstein distribution for s, while neglecting the final state Pauli blocking due to the low density of ν, as is usual in freeze-in computations. Since the issue of bona fide thermalization is quite complicated in any case, this approximation is adequate for our purposes.  In the shaded region, n ν n eq ν such that the reverse process νν → s is important.
The solution to the Boltzmann equation for fixed λ, m s and zero initial n ν is then compared to the equilibrium neutrino density at a given T . If inequality (44) is satisfied for any T above M/3, the freeze-in approximation is adequate. Our numerical results for m s M are presented in Fig. 2. In the shaded region, the neutrino density equals or exceeds its equilibrium value. The kink at roughly 1 GeV appears due to the significant change in g * at the QCD phase transition. We see that only quite small couplings, e.g. below 10 −8 at m s ∼ 1 GeV, are consistent with the freeze-in approximation. The bound can be approximated by The qualitative behaviour of λ(m s ) can be understood from λ 2 -dependence of the rate and the fact that the main contribution to n ν comes from temperatures of order m s . We discuss this in more detail in Sec. 5.2.
In the vicinity of the shaded region, the neutrino density is significant such that the final state quantum statistics (Pauli blocking) can have a tangible impact on the reaction rate. This effect would reduce the rate, hence our bound is somewhat more restrictive than the true one and can be viewed as conservative.
Let us note that other possible "thermalization" conditions appear in the literature. For example, one can compare the neutrino production rate to the Universe expansion rate. If the former dominates, one expects the neutrino sector to be quickly populated. In our case, this corresponds to n −1 s Γ 12 (s → νν) 3H. While such a condition often leads to similar results, there are notable exceptions. In particular, the above inequality is always satisfied at low enough temperatures regardless of the coupling. This, however, does not mean that the neutrino sector gets populated. Indeed, when s is non-relativistic, n −1 s Γ 12 (s → νν) is approximately constant, while H decreases as T 2 . As a result, all the s-quanta available at the corresponding temperature get converted into ν pairs. Yet, since for relativistic neutrinos n eq s (T ) n eq ν (T ), the neutrino density increase is negligible and ν's do not thermalize. Another exception is the situation in which the production is intense but short in duration, e.g. around a phase transition. In this case, the accumulated density can still be small.

Thermalization of s
In this work, we assume that the dominant source of s-quanta is the Higgs thermal bath. It is important to understand under what circumstances the processes h ↔ ss, hh ↔ ss and hh ↔ s lead to thermalization of s. As in the previous section, we use the criterion for any T > ∼ m s to ensure non-thermalization of s. 7 The number density n s is calculated numerically via the Boltzmann equatioṅ where Γ i are the reaction rates Γ 12 (h → ss), The resulting bounds on λ hs are shown in Fig. 3. For a light s, the decay mode h → ss dominates, while for a heavy scalar the production is typically dominated by the fusion mode hh → s. The latter is sensitive to the s-VEV u = m s / √ 2λ s , so additional input such as the coupling λ s is required. This VEV grows very large at small λ s which results in a large reaction rate. Note that in the vicinity of the shaded region, the final state Bose-Einstein enhancement factor can be non-negligible, so our procedure overestimates somewhat the upper bound on the coupling.
The bound on λ hs at m s m h is independent of m s , This is because Γ 12 (h → ss) is independent of m s in this regime and the production stops around T ∼ m h /5. At larger m s , the scattering reaction hh → ss becomes significant. The rate Γ 22 (hh → ss) scales as T 4 in the relativistic regime and the resulting n s (T ) ∝ T 2 . The yield is dominated by low temperatures consistent with the relativistic scaling, that is, T ∼ m s . We thus obtain The fusion channel hh → s is more complicated. For m s 2m h , it becomes active at temperatures below T ∼ m s , that is, when the Higgses still have enough energy to produce s and the Higgs thermal mass is not too large for the process to be kinematically allowed. The fusion becomes inefficient below T ∼ m s /6. In this regime, the reaction rate does not follow any simple scaling law and numerically we obtain The appearance of λ s can be understood from the reaction rate scaling as λ 2 hs /λ s for a fixed m s .

Self-thermalization due to λ s
Even if λ hs is small, the s-sector can thermalize due to self-interaction λ s s 4 . This happens when the number changing processes ss ↔ ssss become efficient and the density n s starts being sensitive to λ s . The specifics of self-thermalization are computationally involved. In the symmetric phase u = 0 at large s-densities close to equilibrium, the (necessary) thermalization condition on λ s has been derived in [47]. Here we are interested in the broken phase u = 0 at low s-densities and thus have to resort to simple estimates. We assume that the initial n s is created via the Higgs thermal bath and study which values of λ s do not affect its evolution.  In the shaded regions, the ss → ssss process is efficient. The bound depends on n s and thus is sensitive to λ hs and the s-production mode. The displayed constraints correspond to two benchmark values of λ hs : 10 −12 and 10 −9 . Here the scalar mixing and the EW transition effects have been neglected.
The 2 → 4 reaction rate at low s-density can be written as where σ 24 is the corresponding QFT cross section and v rel is the relative (Møller) velocity. We are interested mostly in the relativistic regime, in which case σ 24 (ŝ) ∼ 10 −4 λ 4 s ln 2 (ŝ/2m 2 s )/ŝ, whereŝ 4m 2 s is the Mandelstam variable. This result can be verified with CalcHEP [75]. For fixed λ hs and λ s , the density n s (T ) n eq s (T ) is calculated via the Boltzmann equation in the previous subsection.
Although the momentum distribution of s is non-thermal, the characteristic energy of the s-quanta can be approximated by T . This is because n s (T ) is dominated by the late time production in the Higgs thermal bath at temperature T . In the relativistic regime, we may takê s ∼ 4T 2 to calculate the average cross section and v rel ∼ 2. The number changing interactions are efficient if 2Γ 24 > 3n s H, so to ensure non-thermalization we require where we have taken into account the fact that this ratio is maximized at the lowest temperature consistent with the relativistic scaling.
Our numerical results are shown in Fig. 4. Equation (50) makes it clear that the bound on λ s increases with m s . Other qualitative features can be understood from the discussion in the previous subsection. We see that the upper bounds on λ s are significantly above those in [47] (cf. Fig. 2). This is expected since the number density in our case is significantly below its equilibrium value.
Let us emphasise that the above bounds have been obtained under a number of simplifying assumptions. First of all, we have neglected the small scalar mixing, which is not expected to affect the results significantly. We have also assumed that the density of produced particles is low enough such that the final state quantum statistics is unimportant. Finally, we have ignored EW phase transition effects. These can have a non-trivial impact on the bounds. In particular, as we discuss in Sec. 6.2, the hh → s mode can be active even at light m s due to the Higgs mass reduction close to the critical temperature. In this sense, the presented constraints can be viewed as conservative.
5 Sterile neutrino production I: thermalized s

Reaction rates
In the thermal bath of h and s, there are a few channels for ν production, see Fig. 5. The reactions s → νν and ss → νν take place at both high and low temperatures, while hh → νν and hs → νν require the presence of scalar trilinear vertices which only appear below the corresponding critical temperatures.
The relevant interactions are given by Eq. 39. Note that the field VEVs and the degrees of freedom depend on the temperature. At high temperatures, the VEVs vanish, u, v = 0, and the single Higgs d.o.f. is replaced by 4 massive Higgs scalars h i . In this work, we neglect the gauge boson contributions suppressed by an extra power of the gauge coupling (see e.g. [76]). We also neglect the scalar mixing θ 1 apart from the reaction h → νν, which is absent at leading order in θ.  In what follows, we neglect the SM-like Yukawa coupling of the lightest sterile neutrino. As mentioned before, its tiny value can be justified by the neutrino parity.
Below we summarize our results for the reaction cross sections which are to be inserted into Eq. 36 or its equal-mass analog. The masses that appear in the rates are meant to be the thermally corrected masses.

hh → νν
The calculation is easiest performed in the CM frame. The amplitude for the νν final state is Here the combinatorial factor 2 × 2 coming from two identical particles in the initial and final states is included;ŝ = 4E 2 is the Mandelstam variable. The neutrino 4-momenta are denoted by p, p and u, v are 4-spinors.
The spin sum and phase space integration yield where in our convention we include both the initial and final state phase space symmetry factors of 1/2 in the cross section.

hs → νν
The corresponding amplitude is The resulting cross section is As before,ŝ ≡ 4E 2 , although h and s have different energies in the CM frame.

ss → νν
The process ss → νν can proceed both through the s-channel and the t, u-channels at 2d order in λ. The amplitude is wheret,û are the Mandelstam variables. The 4-momenta of the initial state particles are denoted by p 1 and p 2 . The resulting cross section is where the symmetry factors of 1/2 for the initial and final states have been included directly in the cross section. To get 4F σ CM ss→νν , one uses which holds for ss → X processes. We find good numerical agreement with the corresponding CalcHEP [75] result.

s → νν
This process is allowed when m s > 2M . The calculation of the decay rate s → νν is straightforward with the result The consequent reaction rate is We note that the propagators in the above formulas should include the relevant (small) widths to avoid singularities, while the masses include the thermal corrections.

Dark matter abundance: m s > 2M
In this subsection, we solve the Boltzmann equation for the neutrino number density and find parameter regions consistent with the observed DM abundance. Here we assume that m s > 2M such that the decay mode s → νν is available. Note that the thermal correction to M is suppressed by λ 2 and can therefore be neglected.

Qualitative behaviour of the Boltzmann equation solution
Consider freeze-in production of N particles in the reaction i → N . In the relativistic regime, the reaction rate scales as T l , where l depends on the interaction type. Using entropy conservation g * s a 3 T 3 = const with a being the scale factor, one can trade the time variable for T . The resulting Boltzmann equation can be written as where and we have taken the number of d.o.f. to be constant in the range of interest. Assuming that the initial density is zero at temperature T 0 , the solution reads while for l = 5 it is n(T ) = cT 3 ln T 0 T . For renormalizable interactions, l ≤ 4 and the result at late times is insensitive to T 0 : On the contrary, non-renormalizable interactions lead to the "UV freeze-in", i.e. the density dominated by the early time production at T 0 , while for l = 5 n(T ) = cT 3 ln T 0 T . In our work, the relevant reactions are of the type 1 → 2, 2 → 2 and 2 → 1. Their temperature scaling will be discussed later.

Results
The Boltzmann equation describing evolution of the ν number density readṡ HereΓ The theta-functions make sure that the processes involving scalar trilinear vertices are switched off above the critical temperatures. Further, they take care of the different number of Higgs d.o.f. before and after electroweak phase transition. The rates Γ 22 and Γ 12 are calculated according to (36) and (59) using the results of the previous subsections with non-zero v and u. The Higgs decay rate Γ 12 (h → νν) is given by sin 2 θ Γ 12 (s → νν). The Boltzmann equation in the relativistic regime has a simple solution. We find that the most important contribution comes from Γ 12 . Since λ hs and λ s are small, we may neglect the s-thermal mass at late times, in which case Eq. 59 yields while at very high temperatures it scales as m 4 s (T ) ∝ T 4 . In this regime, where the constant is proportional to m 2 s . The dark matter yield is conveniently expressed in terms of Y , where g * s is the number of d.o.f. contributing to the entropy. It is proportional to the total number of the DM quanta. The observed DM density requires The solution (72) is valid in the relativistic regime, that is, down to temperatures of order m s . Thus, the resulting Y ∝ 1/m s . Our numerical results for the total DM relic abundance and the full reaction rates are shown in Fig. 6. We find that the DM yield is dominated by the decay s → νν at temperatures T ∼ m s and the required coupling is This applies to the regime m s M . In this case, the DM yield Y due to the s → νν decay is independent of M and proportional to λ 2 . Thus, in order to get the right relic abundance, the relation λ ∝ 1/ √ M is enforced (while smaller M for the same λ lead to under-abundance). We find that these conclusions apply quite generally, beyond the parameter choices of Fig. 6.
Given the correct relic abundance, small and large values of M are excluded by perturbativity and the Higgs mixing or the presence of a tachyonic scalar. Indeed, since λu = M and m 2 s = 2λ s u 2 , For a fixed relic density and other parameters, λ s ∝ 1/M 3 so that at low M it blows up while for large M it violates 4λ h λ s > λ 2 hs . Since our focus is on freeze-in production of neutrino DM, we exclude significant values of λ. These lead to efficient ν production such that n ν is close to its equilibrium value. In this case, the reverse process νν → s becomes important and the system tends to equilibrate. Although such a possibility is not excluded by observations, it does not correspond to freeze-in neutrino production.
The approximation θ 1 applies in all of the allowed parameter space: θ ranges from 10 −1 in the lower right corner to 10 −5 in the upper left corner of the plots. Close to the tachyonic region however, θ ∼ m h /m s or m s /m h such that the relations (12) receive non-negligible corrections. Therefore the tachyonic region border is only approximate.
The stability condition 4λ s λ h > λ 2 hs combined with the right DM yield Y impose a lower bound on m s , m s > 10 8 λ To get the limit of 1 MeV, we have used λ hs > 4 × 10 −8 required for thermalization and the warm DM bound M > 1 keV (taking the number of SM degrees of freedom at m s to be 10). 8 The main DM production channel is s → νν. We find that the relativistic effects in this reaction are tangible. Fig. 7 shows that replacing the Bose-Einstein distribution with the Maxwell-Boltzmann one can lead to up to a 65% error in the reaction rate. The Bose-Einstein enhancement is sensitive to the thermal mass: for lower couplings the effect is more pronounced. This is natural since the distribution peaks at low energies while the thermal mass provides a lower bound on how low the energy can be.

Light s: m s < 2M
In this case, the main production channel s → νν becomes less significant. The process is kinematically allowed at very high temperatures, when u = 0 and the Majorana neutrino mass vanishes. It stops after the phase transition to u = 0. The produced number density is diluted by the subsequent Universe expansion. As a result, the processes like ss → νν and h → νν become equally important or even take over the leading role.
In case of a very light s, there are a number of non-trivial constraints to be observed. In particular, one must make sure that s decays before BBN. Since s cannot decay into neutrinos, the decay proceeds through the mixing with the Higgs. The decay modes and widths are discussed in Appendix B. We impose the constraint τ < 1sec, which ensures that s does not contribute to the relativistic degrees of freedom at BBN and does not destroy light nuclei. Furthermore, for m s < 2m µ , there is a strong constraint on the mixing angle with the Higgs. Rare Kaon decays require θ 10 −4 [80]. For heavier s, the bound relaxes to 10 −3 or so [81].
Finally, since we are assuming a thermal abundance for s, the Higgs portal coupling must be large enough to ensure thermalization via h ↔ ss, λ hs > 4 × 10 −8 . Although the available parameter space is quite limited, we find that it is still possible to obtain the right DM relic density. Two examples are shown in Fig. 8. In this case, the strongest constraints are imposed by τ s < 1 sec and the absence of tachyons, 4λ h λ s > λ 2 hs . The latter is significant since a light s requires a small λ s . In the allowed parameter space, the bound θ 10 −4 is then satisfied.
As seen in the plots, different reactions dominate at different times. At high temperatures, s → νν dominates but the resulting DM density gets diluted. At later times, h → νν and ss → νν become important. The plateau regions producing the correct DM relic density (Fig. 8, upper row) are associated with ss → νν as the leading (or next-to-leading) production mode. The corresponding rate scales as T 4 down to temperatures of order M . Thus, the resulting yield Since the required Y ∞ also scales as 1/M , the PLANCK line corresponds to a plateau in the (M, λ) plane. The DM yield associated with the different reactions is shown in the lower row of Fig. 8. The left panel confirms that more than 50% of the yield in the plateau region is indeed provided by ss → νν. We also observe that s → νν makes a significant contribution and tilts the PLANCK line in analogy with Fig. 6. At somewhat larger masses, the Higgs decay h → νν becomes more important. The amplitude for this process is proportional to λ θ which is approximately constant for a fixed M : Thus, the resulting PLANCK region is almost vertical in the (M, λ) plane. The lower right panel of Fig. 8 shows that the dominant DM yield is produced at electroweak temperatures via h → νν. To the left of the PLANCK curve, our DM is under-abundant. The neutrino thermalization constraint of Fig. 2 is not directly applicable here since the channel s → νν is not available. We find that, in the allowed parameter region, n ν is below its equilibrium value so the neutrinos can be treated as non-thermal.
6 Sterile neutrino production II: non-thermal s It is possible that s never reaches thermal equilibrium either due to its large mass or due to its small couplings. In general, there is a variety of non-thermal s-production mechanisms in the Early Universe. Its direct coupling to an inflaton would lead to perturbative and/or non-perturbative production, e.g. via parametric resonance [82]. Furthermore, light scalar field fluctuations during inflation generate an s-condensate which then decays into s-quanta. However, these mechanisms are sensitive to further details of the complete UV model, for instance, to the Hubble rate during inflation [83]. In particular, for small Hubble rates such contributions are suppressed. In what follows, we focus on s-production from a Standard Model thermal bath and assume that the other sources are subdominant.

Heavy s
If s is very heavy while the temperature is not high enough, the singlet does not thermalize and can be integrated out. DM production proceeds through Higgs annihilation hh → νν and decay h → νν due to the Higgs-singlet mixing. We find that the decay mode dominates for the parameter values of interest.
It is instructive to consider the channel hh → νν separately. When this mode dominates, one recovers the so-called "UV freeze-in" scenario. In this case, the DM abundance is sensitive to the maximal temperature T 0 < T u c . The Boltzmann equation at high T reads where the factor of 8 takes into account 4 Higgs d.o.f. above the EW transition scale. Since As a result, the DM yield Y = n/s SM ∝ T 0 is determined by the UV end of the evolution. This is unlike the usual freeze-in scenario where the IR behaviour is more important.
Although the decay channel h → νν opens up only below the EW breaking scale, numerically it turns out to be more important and the sensitivity of the DM abundance to T 0 is weak. Our numerical results are presented in Fig. 9 which shows the regions with the right relic abundance. The DM production amplitude is proportional to the combination λθ which is fixed for a fixed M . This makes the production rate independent of λ and the PLANCK region vertical in the (M, λ) plane. As before, our DM is under-abundant to the left of the PLANCK line.
for g * 107. We have verified that the neutrino thermalization constraint is insignificant and n ν is below its equilibrium value.
In the allowed parameter space, the mixing angle ranges from 10 −2 to 10 −6 . As before, the θ 2 corrections become significant close to the tachyonic region border.

Small couplings: freeze-in production of s
Here we consider the possibility that the λ hs and λ s couplings are so small that s never reaches thermal equilibrium (see, e.g. [84] for early work). Assuming that the initial abundance of s is zero or negligibly small, the s quanta are produced by the Higgs thermal bath via the usual freeze-in mechanism. Subsequently, they decay into sterile neutrinos leading to the required DM abundance. Due to the s − h mixing, s decays also produce SM particles, yet this gives only a small correction to the entropy since the density of s is far below its equilibrium value.
There are a few s-production channels: hh → ss, h → ss and hh → s, where the last two reactions are possible only below the corresponding critical temperatures. hh → s is a new reaction type, not considered before. Hence, it is instructive to consider it in more detail.

hh → s rate
The general expression for the reaction rate reads Here |M 2→1 | 2 includes 1/2 from the phase space symmetry of the initial state. Performing the angular integrals as before and using as well as |M 2→1 | 2 = 1/2 × λ 2 hs u 2 , we find This expression is valid for a single Higgs d.o.f.

h → ss and hh → ss rates
These reaction rates have been computed in [48]. For a single Higgs d.o.f., the results read , where E is half the CM energy and we have factored out the symmetry factor 1/2!2! stemming from 2 identical particles in the initial and final states.

Results
The number density of the s-quanta is calculated according tȯ Here the θ-functions account for the EW phase transition and the change in the number of the Higgs d.o.f. We neglect the dependence on T u c since s is not thermalized and λ hs is very small. Since there is no significant back reaction of the produced s quanta on the thermal bath nor substantial entropy production via s-decay, the total DM yield can then be computed as the s-yield times the branching ration for the s decay into dark matter, The s decay width into the SM particles is given in Appendix B.  Figure 10: λ vs M generating the correct relic DM relic density ("PLANCK") for a nonthermal s. Here s is produced by the freeze-in mechanism via the Higgs thermal bath. In the excluded regions, the freeze-in calculations become unreliable due to efficient hh ↔ s or ss → ssss processes leading to thermalization.
We compute the number density of s via freeze-in calculations. Thus, it is important to observe the non-thermalization constraints. For a given λ, an increase in M implies an increase in u, which makes s production via hh → s more efficient and can lead to thermalization. A competitive constraint, which becomes stronger for light s, is imposed by vacuum stability, 4λ s λ h > λ 2 hs . Furthermore, in regions with a substantial λ s , the process ss → ssss becomes efficient and can lead to self-thermalization. We exclude these as well. An additional Kaon physics constraint θ < 10 −4 applies for a very light s, m s < 200 MeV. We find, however, that it is satisfied automatically.
Our numerical results are presented in Fig. 10. The behaviour of the PLANCK curve can be understood as follows. The factors that determine the neutrino abundance are Y s and the decay branching fraction for s → νν. Consider first the regime m s m h . In this case, Y s is determined by the fusion process hh → s, whose rate is proportional to This behaviour is seen in the right panel of the figure. In both panels, DM is under-abundant to the left of (or below) the PLANCK curve. We see that quite large values of λ up to 10 −3 are consistent with all of the constraints. One may worry that the neutrinos would thermalize via s ↔ νν at such a large coupling. However, the density of s is much lower than its equilibrium value and this reaction does not increase the number of s-quanta, while ν → νs is not allowed kinematically and νν → ss is suppressed. Thus, the system is not expected to thermalize.
For a very light s, the BBN constraint on the lifetime of s becomes significant: at small λ, it decays mostly into the photons and electrons which affect the abundance of light elements unless τ s < 1 sec.
Finally, we find that the mixing angle is very small in all the cases considered and its effects can be neglected.

On electroweak phase transition effects
The EW phase transition can have an important impact on the DM abundance. The Higgs mass reduction close to the transition opens up the fusion channel hh → s even if this process is forbidden kinematically at other temperatures. It is operative if 2m h (T T v c ) < m s , while its efficiency depends on the nature of the transition. (An analogous effect in a different setting was considered in [85].) In this work, we are interested in small couplings. Then, the electroweak phase transition corresponds either to a second order phase transition or a crossover. In the former case, the Higgs becomes massless at the critical temperature, while at the crossover it remains massive. Perturbative analysis is insufficient to distinguish the two: what appears as a second order transition typically corresponds to a crossover, as established by lattice simulations. The full analysis of the singlet scalar extension is not yet available, although for a heavy singlet or EW triplet, the nature of the transition has been determined in [86,87,88]. The second order transition is found to occur in special cases, while a crossover is very common at weak coupling. This is to be contrasted with perturbative calculations (see e.g. [89]). Similar results are expected to apply in the light singlet or triplet case. 9 Although the Higgs does not turn massless at the crossover, its mass gets significantly reduced. In the SM, this reduction reaches an order of magnitude at the (pseudo-)critical temperature [90] (see also earlier work [91,92]). Since we are mostly interested in very small Higgs portal couplings, the presence of the singlet is not expected to change the nature of the transition. Thus, we may assume m h (T v c ) ∼ 10 GeV as in the SM. To estimate the efficiency of the fusion mode, let us consider a simplified case of zero s − h 9 We thank Lauri Niemi for sharing some of his results.
mixing and employ a simple parametrization 10 where T c ≡ T v c is the EW critical temperature and c is a constant fixed by requiring m h (0) = 125 GeV. Taking a simple perturbative estimate for T c , one can then calculate the fusion rate. The resulting hh → s rate for a representative parameter set is shown in Fig. 11, left panel.
We find that this effect alone can account for all of the observed dark matter. Although short, the fusion is intense enough to produce numerous s-quanta which subsequently decay into sterile neutrinos. As one gets closer to the 2d order transition (at larger λ hs ), the Bose-Einstein enhancement becomes more pronounced. This is illustrated in Fig. 11, right panel. When both m h (T c ) and m s are far smaller than the temperature, the Bose-Einstein enhancement factor can reach orders of magnitude.
The fusion mode can be more efficient than the decay h → ss. Indeed, the fusion rate grows as u 2 which can be very large, while the decay rate remains constant for a fixed m s . Thus, the thermalization constraints in Figs. 3,4 due to the fusion mode extend to m s < 2m h as well and can be more stringent then those due the decay, depending on u. However, in view of the uncertainties, we have not included these to be conservative.
We note that our approximation breaks down at m s ∼ m h , i.e. when the mixing angle becomes significant. As pointed out in [76], the resonantly enhanced s − h mixing leads to additional scalar production. With present tools, it is however difficult to estimate its efficiency and we leave it for future work. We stress that the fusion mechanism considered here is intrinsically different and operative for small (and zero) mixing as long as 2m h (T ) < m s .

Conclusion
The lightest sterile neutrino is an attractive dark matter candidate. Although it is not stable, its longevity is guaranteed by its small mass and a small sterile-active mixing angle. In this work, we explore the mass range up to 1 GeV. In this case, tiny mixing angles are necessary which one can justify by a flavor-dependent (neutrino parity) symmetry.
We have focused on the scenario where the Majorana masses are entirely due to a VEV of a real scalar. This is enforced by a discrete lepton number symmetry, which is broken spontaneously by the scalar VEV. The scalar is then only allowed to couple to the SM quadratically through the Higgs portal.
Since the neutrinos can be very weakly coupled, the natural (but generally not unique) dark matter production mechanism is the freeze-in. We have analyzed freeze-in production of sterile neutrinos (ν) from the Higgs and singlet scalar (s) thermal bath in different regimes. These are summarized in the following In all of these cases, the observed DM relic density can be obtained. For the sterile neutrino mass range (1 keV, 1 GeV), we find that the requisite scalar-neutrino coupling varies between 10 −10 and 10 −3 . Our analysis takes into account the relativistic reaction rates with the Bose-Einstein distribution function, thermal masses and main effects of the phase transitions. All of these factors make an important impact on the final results. As byproducts, we have derived relativistic rates for asymmetric reactions as well as non-thermalization constraints on sterile neutrinos and the Higgs portal scalar.
We find a number of interesting effects which deserve further study. In particular, a light scalar can be copiously produced close to the EW phase transition/crossover through the fusion mode hh → s. Subsequent decay of the scalar into sterile neutrinos can account for all of the dark matter. However, the specifics of this mechanism require understanding non-perturbative dynamics close to the critical temperature.
The dark matter candidate studied here is long-lived. Its production mechanism is independent of the sterile-active mixing Θ, hence there is vast parameter space (Θ, M ) where dark matter decay can lead to an observable signal, e.g. in the form of monochromatic X-or gamma rays. The intensity of the signal is correlated with the dark matter density.
Acknowledgements. OL is indebted to Mark Hindmarsh, Lauri Niemi and Aleksi Vuorinen for invaluable discussions. VDR acknowledges financial support by the SEJI/2018/033 grant, funded by Generalitat Valenciana and partial support by the Spanish grants FPA2017-85216-P and FPA2017-90566-REDC (Red Consolider MultiDark). DK is supported by the National Science Centre, Poland, research grant No. 2015/18/A/ST2/00748. This work was made possible by Institut Pascal at Université Paris-Saclay with the support of the P2I and SPU research departments and the P2IO Laboratory of Excellence (program "Investissements d'avenir" ANR-11-IDEX-0003-01 Paris-Saclay and ANR-10-LABX-0038), as well as the IPhT. TT acknowledges funding from the Natural Sciences and Engineering Research Council of Canada (NSERC). Numerical computation in this work was carried out at the Yukawa Institute Computer Facility.

A Leading thermal corrections
In this Appendix, we summarize the most important thermal corrections to the effective potential in our model.
The tree-level effective scalar potential, written in terms of the vevs v, u reads where the sum runs over fermions and W inside the loop. In this expression, N f (W ) c = 3(1), Q i is the charge and G F is the Fermi coupling constant. We define the following mass ratio and the loop functions A W (τ ) = −(2τ 2 + 3τ + 3(2τ − 1)f (τ ))/τ 2 , with for τ > 1. For heavier m s , the scalar will also decay into other SM particles. Besides the SM channels, s has another important decay mode s → νν. The corresponding decay width reads Fig. 12 shows the total SM decay width and Γ(s → νν) as a function of m s with other parameters fixed at some representative values. While the neutrino width grows with m s , the SM decays get suppressed due to the decrease in the mixing angle θ ∝ 1/m 2 s . The spike in Γ s around m h m s is due to the sharp increase in sin θ. In this region, our approximations are unreliable.