Dark Matter from Freeze-In via the Neutrino Portal

We investigate a minimal neutrino portal dark matter model (DM) where a right-handed neutrino both generates the observed neutrino masses and mediates between the SM and the dark sector, which consists of a fermion and a boson. In contrast to earlier work, we explore regions of the parameter space where DM is produced via freeze-in instead of freeze-out motivated by the small neutrino Yukawa couplings in case of $\mathcal{O} \left( \mathrm{TeV} \right)$ heavy neutrinos. For a non-resonant production of DM we can predict the heavy neutrino mass to be $M_N \approx 10 \, \mathrm{TeV}$ and achieve a lower bound on the DM mass of $m_\chi \gtrsim 5 \, \mathrm{TeV}$. For the resonant production of DM, we find that it can be produced via freeze-in or freeze-out even with couplings of $\mathcal{O} \left( 10^{-5} \right)$. In the latter case to produce the observed for DM density we find $m_\chi \approx 100 \, \mathrm{eV}$.


Introduction
Both Dark Matter (DM) and neutrino masses provide a strong hints for beyond standard model physics (BSM). A simple way to accommodate for neutrino masses is to introduce right-handed neutrinos which are SM singlets, thereby allowing for a Majorana mass term. This enables mass generation via the type I seesaw mechanism. Furthermore, the resulting heavy neutrino state N is massive, electrically neutral and within a certain mass range even (cosmological) stable. If the heavy neutrino is considered to be a DM candidate it has to be stable. To ensure that its mass must satisfy M N < 2m e . Therefore, the Yukawa coupling has to be very small, namely y ν 10 −6 . As a consequence, the production rate is small, thereby allowing for DM production via the freeze-in mechanism 1 [2]. In freeze-in scenarios, DM production never becomes efficient, i.e. the interaction rate Γ is always small compared to the Hubble parameter H, Γ H. Thus, DM is never in thermal equilibrium with the SM. On the contrary, freeze-out scenarios require Γ > H for a certain period of time (see figure 1). To account for the observed DM relic abundance via freeze-in of the decay h → N ν, the heavy neutrino mass should be of O (10 keV). However, the possibility of keV sterile neutrino DM [3] might be excluded by experiments soon [4], more precisely by the non observation of the decay N → νγ in gamma rays, see e.g. [5]. Moreover, in case of M N > 2m e , the heavy neutrino N is not stable and therefore not a DM candidate. But even in this case, the right-handed neutrino can act as a mediator to DM since it is a SM singlet. This possibility is referred to as neutrino portal DM (NPDM) which has recently drawn attention in the literature [6][7][8][9]. However, these works only consider DM production via the freeze-out mechanism. This paper explores a minimal NPDM model where DM is produced via the freeze-in mechanism. The paper is organized as follows: In section 2, we present the model and its general properties. In section 3, we briefly review the Boltzmann equations which are solved for the considered model analytically within some special cases in section 4. Section 5 summarizes the results of a numerical solution of the Boltzmann equations. In section 6 we conclude.

Setup
In this section, we introduce the particle content and the general properties of the model. A model with the same particle content was investigated in [7] where DM production within freeze-out scenarios was explored. In addition to the SM particle content, the model includes a right-handed neutrino ν R to accommodate for the observed neutrino masses via the type I seesaw mechanism. Moreover, the model contains a dark sector consisting of a fermion χ and a scalar φ. While they are uncharged under the SM gauge groups, they are charged under a dark symmetry, e.g. a dark U(1). Assuming the SM particles to be uncharged under the dark symmetry renders the lighter particle of χ and φ to be a stable DM candidate since the dark symmetry forbids couplings between SM and dark sector particles. In this scenario, the resulting heavy neutrino N acts as a mediator between the DM and the SM particles since the singlet ν R can couple to χ and φ via a Yukawa coupling as long as the expression χφ is a singlet under all gauge groups. The parts of the Lagrangian relevant for the neutrino mass generation and the coupling to DM are given by Since the focus of this work is the DM production mechanism we restrict our analysis to only one generation of SM and DM particles. We do not take into account any contribution to the DM relic abundance from a possible Higgs portal interaction arising from the term (φφ * ) (hh * ) in the scalar potential and furthermore assume that φ does not acquire a VEV. Moreover, effects resulting from kinetic mixing of the vector mediators of the dark symmetry with the SM gauge bosons are neglected, thus our analysis focuses on the neutrino portal to DM only. As already mentioned, the observed neutrino masses are generated by a type I seesaw mechanism leading to mass eigenstates ν and N with masses of approximately m ν = y 2 where v is the VEV of the Higgs. They are described in terms of the interaction eigenstates by: This mixing causes an interaction between the ν, N and the Higgs as well as a coupling of N to the SU (2) L gauge bosons. As presented in [10], the resulting interactions between the heavy and the light neutrinos are given by: The matrices B and C are defined as in [10] and in case of only one generation and real Yukawa couplings they yield: Thus, the couplings relevant for heavy neutrino production are given by whereas the coupling of the heavy neutrino to the dark sector is governed by: Note that the parameters y ν and M N are not independent and related by the seesaw mechanism requiring y ν = √ m ν M v −1 . Therefore, the couplings in the eq. (2.7)-(2.9) can be rewritten as:

Boltzmann Equations
Determining the relic abundance of the DM candidate requires solving the Boltzmann equations, which describe the time evolution of the particle number densities in the expanding universe. Adopting the formalism used in [11] the Boltzmann equations can be written as: Here, n i is the number density of particle species i and H is the Hubble parameter. The 3Hn N term takes the expansion of the universe into account while the right hand side governs the impact of scattering processes which occur with a certain thermal rate γ eq . The equilibrium number densities n eq i are given by the momentum integral over the distribution function f of the respective particle species which is approximated with a Boltzmann distribution in our case: For a two to two scattering involving only CP conserving interactions the quantity γ eq results in For a CP-conserving decay we obtain Here z = M N T holds and Γ N is the decay rate of the particle in its rest frame. Next, to simplify the form of the Boltzmann equations they are written in terms of the quantity Y = n s E instead of the number density where s E = 2π 2 g s eff 45 T 3 is the entropy density. This leads to Finally, for a numerical solution, it is useful write eq. (3.5) in terms of log 10 (Y N ) and log 10 (z): γ eq (N a · · · ↔ ij . . . ) Hn eq N n N n a . . . n eq N n eq a . . .

Relic Abundance: Analytic Estimates
The 2 ↔ 2 scattering processes responsible for producing DM (the lighter particle of χ and φ) can be classified into two categories: SM Particle Scattering and Heavy Neutrino Scattering.
The SM particle scattering processes involve two SM particles in the initial state, have χ and φ in the final state and are mediated by the heavy neutrino. Consequently, we have σ ∼ y 2 ν y 2 χ . The heavy neutrino scattering processes have two heavy neutrinos in the initial state and produce a pair of χ or φ. Here, we have σ ∼ y 4 χ . In addition to the different dependence on the couplings y ν and y χ , for SM particle scattering processes we can assume the SM particles to be in thermal equilibrium whereas the heavy neutrino N is only in thermal equilibrium if the production via (inverse) decays as e.g. (νh → N ) h → νN is sufficiently efficient. Consequently, even in the case of y χ y ν the heavy neutrino scattering might not be the dominant production channel since n N n eq N 1 could be the case. All contributing diagrams are displayed in figure 2.

SM Particle Scattering
For the rest of the discussion, we assume that the dark sector particles are roughly the same in mass and therefore replace m φ = m χ and furthermore assume M N > M W . As discussed in the end of chapter 2, for M N M W the coupling of the heavy neutrino to the Higgs and a light neutrino is much larger than the coupling to the SU (2) gauge bosons. Therefore, we only take the contribution of vh ↔ φχ into account. The relevant cross section is given in appendix A (A.1) and for m φ = m χ results in Here, Γ N is the total decay width of the propagating neutrino which can decay into vh for M N > m h and into χφ for M N > m χ + m φ . The decay width is also given in appendix A (A.2). Next, we use eq. (3.3) to determine the thermal rate of the process. There are two cases to be distinguished: • The resonant case with M N ≥ m χ + m φ where M 2 N ≥ s min and thus it is integrated over the resonance of the cross section σ vh↔χφ First, we discuss the non-resonant case where we neglect the contribution of the decay width Γ N . Moreover, for M N m χ the integral in eq. (3.3) is solvable analytically and results in: Therewith, eq. (3.5) for the dark sector particles results in: The number density of DM particles is given by n χ + n φ and we assume all SM particles to be in thermal equilibrium, i.e. n SM = n eq SM . Furthermore, since we are investigating a weakly interacting dark sector we assume n χ(φ) n eq χ(φ) during the time of production. Thus, The factor of 2 arises due to the assumption that the slightly heavier particle of χ and φ is long lived but might eventually decay to the stable one, hence Y DM = Y χ + Y φ . Integrating from z = 0 to z = ∞ yields: Remarkably, the result is inverse proportional to the DM mass m χ , i.e. the energy density is independent of m χ . This allows for predicting the value of the product of the Yukawa couplings y ν y χ by setting Y DM (z → ∞) = Y DM,exp , with The experimental values for Ω DM , Ω B , the density parameter for baryons, and Y B ,the baryon number density in a co-moving volume, are taken from [12] and m B , the average baryon mass, is approximated with the proton mass. Evaluating Y DM (z → ∞) = Y DM,exp results in: The implications of this result are discussed in chapter 4.3 Next, we discuss the resonant case, i.e. M N ≥ m χ + m φ . As it was pointed out in [13], in this case it is useful to approximate the Breit-Wigner peak in eq. (4.1) with: which is valid as long as b a, i.e. Γ N M N . With that we find The integration of the approximated Boltzmann equation (4.4) results in: Since the result is not proportional to m −1 χ fitting the observed DM density depends on the DM mass m χ itself. Again, we postpone the discussion of the result to chapter 4.3.

Heavy Neutrino Scattering
The cross sections for the heavy neutrino scattering processes are given in Appendix A. For the case of M N m χ ≈ m φ they result in 11) σ N N →φφ = y 4 Although it is not possible to find an analytic expression for the process N N → φφ it turns out that σ N N →χχ ≈ σ N N →φφ for all s. Therefore, we approximate its contribution to the DM production in the derivation of an analytic result by the contribution of the process N N → χχ. The Boltzmann equation for the DM density by again assuming n χ n eq χ and the factor four arises since two χ/φ particles are produced and due to γ eq (N N → φφ) ≈ γ eq (N N → χχ). Furthermore, we assumed n N = n eq N . This assumption is only valid if the processes that keeps N in thermal equilibrium, e.g. the inverse decay νh → N , are sufficiently efficient. The integration of eq. (4.14) yields: As for the SM particle scattering in the limit of M N m χ , the DM density is inverse proportional to its mass and thus allows for a prediction of y 4 χ ≈ 10 −5 . However, for a realistic result we have to consider both -SM particle scattering and heavy neutrino scattering -processes. This is discussed in the next chapter. For the case where the SM scattering processes are in the resonant regime, i.e. M N > m χ +m φ , the cross section given by eq. (A.5) consists of two terms. The first can be integrated in the limit M N m χ to obtain the thermal rate (4.16) Even within this limit it is not possible to obtain an analytic result for the second term in the cross section. However, the contribution of this term to the cross section is smaller or equal compared to the contribution of the first term for all center of mass energies s. Hence we approximate the thermal rate by two times (4.16). Here, the integration of eq. (4.14) results in: (4.17)

Discussion of the Analytic Results
In the limit of M N m χ ≈ m φ we found analytic solutions for the DM relic density for both types of processes. Combining both results yields: Therefore, the coupling y χ is required to be smaller than 10 −5 in order to not overproduce DM. In principle, the couplings y χ and y ν are unrelated. However, both describe a coupling to the right-handed neutrino and -if the heavy neutrino is lighter than O (10 15 GeV) -both couplings are required to be relatively small. This motivates the idea that they might be suppressed by the same mechanism, thus resulting in y ν ≈ y χ . 2 . Considering a model which generates y χ ≈ y ν allows for constraining the mass of the heavy neutrino since with that eq. (4.19) reads 5π + 2 2 (4.20) Thus, to fit the observed DM density (4.6), M N ≈ 10 TeV is required. However, we achieved this result by assuming that the heavy neutrino is always in thermal equilibrium, m χ M N and by only taking into account the dominant processes of the SM particle and heavy neutrino scattering each. From eq. (4.18), we see that the contribution of the heavy neutrino scattering processes only accounts for roughly ten percent of the produced DM in case of y χ = y ν . Thus, the result will not be altered significantly if the heavy neutrino is out of equilibrium. Also, taking into account the sub-dominant processes does not have a significant impact since they are suppressed by The only significant change will occur in areas of the parameter space where m χ ≈ M N , thereby violating the assumption of m χ M N . For these reasons in section 5, we solve the Boltzmann equations numerically for the case of y ν = y χ . In addition, we found an analytic solution for the DM relic density in the limit M N m χ where the SM particle scattering processes are in the resonant regime: Again, we obtain an upper bound on y χ from the measurement of the DM energy density by setting y ν = 0 which results in y χ 10 −5 M N m −1 χ 1 4 . Moreover, in case of y χ y ν we find the observed DM energy density if y χ ≈ 10 −12 M N mχ . However, if y χ y ν does not hold the approximation of n χ n eq χ we used to derive (4.18) does not apply anymore. To illustrate that we look at the case y χ = y ν , where (4.18) results in:  Using eq. (3.2) we find that Y eq DM 10 −2 . Therefore, n χ n eq χ cannot be satisfied. Hence, the freeze-in scenario does not apply here. Nevertheless, it is still possible to account for the correct amount of DM. In this case, we recover a freeze-out like scenario since due to the resonance the interaction rate becomes as large as the Hubble parameter although the system is only weakly coupled. Thus, DM comes into equilibrium with the SM and freezes out as soon as the interaction rate becomes smaller than the Hubble parameter. This occurs approximately at T = M N . 3 Consequently, the number density can be estimated by the equilibrium density at freeze-out and since we assumed M N m χ the number density (3.2) at freeze out yields: Equating this result with eq. (4.6) yields a DM mass of m χ = O (100 eV). We summarized our results for the case y χ = y ν in a schematic plot (see fig. 3).

Numerical Analysis
In this section we present the results of the numerical solution of the Boltzmann equation for the DM candidate χ. We solved (3.5) for χ and N in the non-resonant regime, including the processes νh ↔ χφ, N N ↔ χχ, N N ↔ φφ and N ↔ νh or h ↔ N ν. Furthermore, we set y ν ≈ y χ and m φ = m χ . As a consequence of our analytic approximations we expect the scenario to work for M N ≈ 10 TeV and m χ 5 TeV. However, especially if M N ≈ m χ our approximations from the analytic computation do not apply. Also N , as shown in appendix B will reach thermal equilibrium via the process N ↔ νh around T = M N . Thus, the contribution of N N ↔ χχ and N N ↔ φφ might be smaller than expected. We solved the problem numerically for M N < 2m χ while assuming that during the time of the DM production n χ = n φ and that φ eventually decays into χ so that the DM relic density is given by n χ + n φ . The initial conditions are that the SM particles are in thermal equilibrium while N , χ and φ are not present. The interpolation of the thermal rates and the numerical solution of the Boltzmann equations are performed with Mathematica. The obtained results are summarized in figure 4 where we compared the obtained DM density Y th to the experimentally observed DM density Y exp for different values of m χ and M N . As expected, for M N m χ the correct density is obtained for a constant value of M N of M N ≈ 30 TeV. Also we see that, in fact, there is a lower bound on the DM mass which results in M N ≈ 5 TeV. Lastly we comment on the case y ν = y χ : Here, the coupling y χ which is determined by the mediator mass via eq. (4.19) is larger than y ν when M N 10 TeV. In this case, most of the DM gets produced via heavy neutrino scattering. Since we overestimated its contribution to the relic density by assuming n N = n eq N a larger coupling y χ will be required and therefore relax the upper bound on y χ 10 −5 . The numerical analysis of the resonant case is not performed within this paper. As in the non-resonant region, we expect the most significant deviations from the analytic result for cases where M N ≈ m χ .

Conclusion
We have investigated a minimal neutrino portal DM model where the dark sector consists of a scalar φ and a fermion χ. In addition, the SM is extended by a right-handed neutrino which generates the neutrino masses via a type I seesaw mechanism and, furthermore, acts as a mediator between the SM and DM. Motivated by the small Yukawa couplings of the type I seesaw mechanism in case of small heavy neutrino masses of M N O (PeV) we studied DM production via the freeze-in mechanism which, in contrast to the freeze-out mechanism, requires interaction rates of Γ H. We derived analytic solutions for the resonant (M N > m χ + m φ ) and non-resonant (M N < m χ + m φ ) DM production regime. Adding the requirement that the coupling of the righthanded neutrino to the SM is of the same order of magnitude as its coupling to the dark sector increases the predictivity of the model and allows for predictions of the mediator or the DM mass respectively. In the non-resonant regime, we find M N ≈ 30 TeV and a lower bound on the DM mass m χ 5 TeV. Within the resonant regime, however, for y χ = y ν the resonant production of DM is strong enough to bring DM into equilibrium with the SM. Thus, we recover the usual freeze-out mechanism although the couplings between DM and the SM are very small. Moreover, in this scenario we can predict a DM mass of m χ ≈ 100 eV. For y χ y ν , nonetheless, DM production via freeze-in is still possible. To satisfy the observed DM energy density the coupling of the right-handed neutrino to DM has to be y χ ≈ 10 −12 M N mχ . Thus, we demonstrated that producing the observed DM energy density within this model of neutrino portal DM is viable even with small couplings between the SM and the dark sector. An interesting extension of this model would be the explanation of the small Yukawa couplings of the right-handed neutrino to the SM and DM by a suppression mechanism. One explanation could be an extra dimensional model where the heavy neutrino in contrast to the SM and DM particles propagates in an extra dimension since it is the only singlet under all considered gauge groups. Thus, all couplings to the right-handed neutrino are suppressed by a volume factor which can even generate y χ ≈ y ν . Due to the tiny couplings of the right-handed neutrino to SM and DM direct detection seems not viable. Therefore, exploring the detectability of this model might be an interesting prospect for a future work. Moreover, we only demonstrated a few working cases of this model. For example we considered only the case of m χ = m φ and M N > M W . Thus, an exploration of M N < M W which leads to different dominant interactions between the heavy neutrino N and the SM or a hierarchic dark sector might be interesting. Here, Γ N is the total decay width of the propagating neutrino which can decay into vh for M N > m h and into χφ for M N > m χ + m φ . The decay width is given by: The decay rate for the process N → νh for M N > m h is given by:

A Cross Sections
With that we find the interaction rate: Comparing this result to the Hubble Parameter with M N m h yields: This quantity is greater than one for Temperatures of roughly T ≤ M N . Consequently we can expect the heavy neutrino to be in thermal equilibrium with the SM for T M N .