Dark matter from freeze-in via the neutrino portal

We investigate a minimal neutrino portal dark matter (DM) model 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 OTeV\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {O} \left( \mathrm {TeV} \right) $$\end{document} heavy neutrinos. For a non-resonant production of DM, its energy density is independent of the DM mass. Assuming a democratic coupling structure we find MN≈10TeV\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M_N \approx 10 \, \mathrm {TeV}$$\end{document}. For the resonant production of DM, we find that it can be produced via freeze-in or freeze-out even with couplings of O(10-5)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {O} (10^{-5} )$$\end{document}. However, the measurement of the Lyman-α\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha $$\end{document} forest rules out the feeble coupled freeze-out case completely, while the resonant freeze-in production is only viable for mDM≳3keV˚\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m_{DM} \gtrsim 3 \, \mathring{keV}$$\end{document}.


Introduction
Both dark matter (DM) and neutrino masses provide strong hints for beyond standard model physics (BSM). A way to accommodate neutrino masses is to introduce right-handed neutrinos as SM singlets, thereby allowing for mass generation via the type I seesaw mechanism.
Furthermore, the resulting heavy neutrino state N is massive and electrically neutral. If it is considered to be a DM candidate it must be stable. Thus, its mass must satisfy M N < 2m e . Therefore, the Yukawa coupling has to be very small, namely y ν 10 −6 . Consequently, the production rate is small, allowing for DM production via the freeze-in mechanism 1 [2,3].
In freeze-in scenarios, DM production never becomes efficient, i.e. the interaction rate is always small compared to the Hubble parameter H , H (see Fig. 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 via freeze-in within a minimal setup, the Dodelson-Widrow mechanism [4], is already excluded by the experiment, more precisely by the non-observation of the decay N → νγ [5,6] and Lyman-α measurements [7]. However, the idea of sterile neutrino dark matter via different production mechanisms continues to be widely discussed [8].
In case of M N > 2m e , the heavy neutrino N is obviously 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, a possibility which is referred to as neutrino portal DM (NPDM) [9][10][11][12].
Within these works, the small neutrino masses are generated by the type I seesaw mechanism and DM is produced via the freeze-out mechanism. In contrast, this work explores a minimal NPDM model where DM is produced via the freezein mechanism.
In Sect. 2, we introduce the particle content and the coupling structure of the model. In Sect. 3 the method for deriving the analytic results for the DM number density while assuming a thermal shape of the distribution function is introduced. Although those analytic results, which are discussed in Sect. 4, are not exact they allow for studying the parametrics for DM production. Following in Sect. 5 we numerically solve the Boltzmann equations at the level of momentum distribution functions taking the non-thermal form of the momentum distribution into account. Section 6 summarizes the relevant constraints on the model from direct detection, lepton flavour violation and structure formation. After that we conclude.
Within the appendices, the relevant reduced cross sections are given and the method for solving the boltzmann equations at the level of momentum distribution functions is discussed in more detail. between the freeze-out case (green) and the freeze-in case (red) results from the much smaller coupling in the freeze-in case. The right panel shows the corresponding number densities compared to the equilibrium density in a co-moving volume 2 Setup A model with similar particle content was investigated in [10], where DM production within freeze-out scenarios was explored. In addition to the SM particle content, the model includes three right-handed neutrinos ν R i to accommodate the observed neutrino masses. The dark sector consists 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) or a Z 2 . 2 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 neutrinos N i mediate between the DM and the SM particles since the singlets ν R i 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 Neutrino mass generation − y χ φχν R i DM coupling +h.c..

(2.1)
Here, we assumed a universal coupling of DM to the three right-handed neutrinos. Furthermore, 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 additionally assume that φ does not acquire a VEV. 3 Moreover, effects resulting from kinetic mixing of possible vector mediators of the dark symmetry with the SM gauge bosons are neglected. Thus, our analysis focuses on the neutrino portal to DM only. After electroweak symmetry breaking the observed light neutrino masses are generated via the type I seesaw mechanism. To ensure that the observed neutrino masses and mixing angles are reproduced we utilize the following parametrization of the Yukawa coupling matrix Y ν [15]: where we assumed the Majorana mass matrix M M to be diagonal with degenerated eigenvalues, i.e. M M = diag (M N , M N , M N ). U PMNS is the PMNS matrix, v is the vacuum expectation value of the Higgs field, √ m ν is a diagonal matrix with the square root of the neutrino masses as eigenvalues, R is an orthogonal complex 3 × 3 matrix and m ν is the square root of the large mass squared difference m ν = m 2 ν . The mass-and interaction eigenstates are transformed into each other in leading order in the small parameter y ν v M −1 N by the matrix U : 3 In fact, the validity of this assumption as well as the vacuum stability of this model will be investigated in a future work since due to a fermion loop consisting of a ν R and a χ the φ mass term receives a negative contribution. In case the fermions in the loop are heavy compared to the boson those radiative corrections might lead to a negative m 2 φ and thereby break the symmetry that stabilizes DM. Similar effects have been investigated for the scotogenic model [13,14] where those effects constrain the parameter space significantly.
The mixing between the left and right handed neutrinos causes an interaction between ν, N and the Higgs as well as a coupling of N to the SU (2) L gauge bosons. As presented in [16], the resulting interactions between the heavy and the light neutrinos are given by: The matrices B and C are defined as in [16] and in case of real Yukawa couplings, as we will assume no CP violation from now on, 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 N v −1 . Therefore, the couplings in Eqs. (2.8)-(2.10) excluding the flavor dependent part can be rewritten as: Thus, for M N ≥ M W , the coupling g hν N can be expected to be dominant and the hν N vertex is the most relevant one for DM production. Whereas for M N ≤ M W , the Wl N and Z ν N vertices are expected to contribute the most to DM production as long as M N m ν .

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. In principle, the boltzmann equations have to be solved at the level of momentum distribution functions, which then are integrated to obtain the number density. For a freeze-out production of DM however those distribution functions can be safely assumed to be proportional to a Boltzmann distribution, which allows for solving the Boltzmann equations at the level of number densities directly. Although this assumption can lead to less precise results in case of freeze-in production we will still use this formalism to obtain analytic expressions for the relic density in Sect. 4. Later on in Sect. 5, a numerical solution of the Boltzmann equation is given at the level of momentum distribution functions.
Here, we review the formalism for solving the Boltzmann equation for number densities, while the one for distribution functions is discussed in Appendix A.
Adopting the formalism used in [17], the Boltzmann equations can be written aṡ Here, n i is the number density of particle species i. 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 eq i 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 with z = M N T . For the special case of freeze-in production via a twoto-two scattering process b 1 b 2 → i j the solution of this equation can be written in a compact form. Here, b 1/2 are particles in thermal equilibrium with the SM, whereas the number densities of i and j satisfy n i/j n eq i/j . Then, the Boltzmann equation for the particle species i is given by: Inserting γ eq (3.3) and integrating the equation from very large temperatures, i.e. z → 0, up to today, i.e. z → ∞, yields: (3.6) Here we use K = Hs E T −5 and z = m i T −1 . After performing the z integration with the initial condition Y i (z = 0) = 0 we are left with 4 The 2 ↔ 2 scattering processes responsible for producing DM 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 χ . All contributing diagrams are displayed in Fig. 2. The following discussion assumes only one SM and right-handed neutrino generation. However, these results can easily be translated into a three generation setup due to the assumption of degenerated heavy neutrino masses, i.e. M N i = M N and the universal coupling of the dark sector to the right-handed neutrinos. For the heavy neutrino scattering, the one generation result has to be multiplied by a factor of nine. For the dominant SM particle scattering process ν i h → χφ via a N j the one generation contribution with a neutrino Yukawa coupling of y ν = √ m ν M N v −1 has to be multiplied by where θ is a vector containing the in our case three real angles parametrizing the orthogonal matrix R. Choosing the standard parametrization for an orthogonal three by three matrix we find 10 −16 f 1 (θ ) ≤ 3. Since the Z ν i N j vertex has the same flavor structure as the hν i N j vertex the one generation result for the Z ν initial state is multiplied by the same factor as the hν initial state.
Only for the Wl initial the factor differs and results in Here, we find 10 −18 f 2 (θ ) 7.65. Scanning both f 1 and f 2 for randomly chosen values for the angles θ shows that on average f 2 ≈ 2.5 f 1 . Nevertheless, excluding the cases where f 1 is close to its lower bound, the contribution of the hν i initial state is still the dominant one due to the following reason: The production via the scattering of the gauge bosons is only viable for temperatures below the critical temperature where the SU (2) L × U (1) Y symmetry of the SM gets broken. Hence, the time of production is small compared to the Higgs neu-trino scattering. Therefore, we consider only the production via hν i → χφ for the analytic estimates, while all production channels are taken into account in the numerical solution.

SM particle scattering
For the rest of the discussion, we assume that the dark sector particles have roughly the same mass and replace m φ = m χ . The reduced cross section for the dominant production channel is given by: .
Here, N is the total decay width of the propagating neutrino. There are two cases to be distinguished: First, we discuss the non-resonant case. If we neglect the contribution of the Higgs mass, i.e. m h m χ , we can use Eq. (3.7) to determine the relic density directly: eff are the number of effective relativistic (entropy) degrees of freedom which are both assumed to be constant during this calculation with g (s) eff = 106.75. 5 Note that for obtaining this result the reduced cross section was multiplied by an additional factor of four arising from the four degrees of freedom of the Higgs doublet before the electroweak phase transition.
Remarkably in case of a heavy DM mass m χ compared to the mediator mass M N , the result is inversely proportional to the DM mass, 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 The experimental values for DM , the density parameter for baryons B , and the baryon number density in a co-moving volume Y B , are taken from [18] and m B , the average baryon mass, is approximated with the proton mass. Evaluating Y DM = Y DM,exp results in: The implications of this result are discussed in Sect. 4.3 Next, we discuss the resonant case, i.e. M N ≥ 2m χ . As it was pointed out in [19], 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 . Then, the integration of Eq. (3.7) results in: where we already used M N m χ to simplify the result. Again, we postpone the discussion of the result to Sect. 4.3.

Heavy neutrino scattering
The cross sections for the heavy neutrino scattering for the case of M N m χ result in (4.9) By again employing Eq. (3.7) we find: As for the SM particle scattering in the limit of M N m χ , the DM density is inversely proportional to its mass.
For the case where the SM scattering processes are in the resonant regime, i.e. M N > 2m χ , in the limit M N m χ we cannot find an analytic estimate for the DM relic density beside Although the factor of proportionality is unknown we expect this to be much smaller compared to the contribution of the SM particle scattering. This is due to the resonance contributing to the production via SM particle scattering. Hence, we neglect this contribution for the discussion of the analytic results.

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: (4.12) By comparing this expression with the observed DM density (4.4) one obtains Since the coupling y ν is only a function of M N the coupling y χ is fixed by the heavy neutrino mass M N . Moreover, we find y χ 10 −5 in order not to overproduce DM. In principle, the couplings y χ and y ν are otherwise unrelated. However, both describe a coupling to the righthanded 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, resulting in y ν ≈ y χ . 6 Considering a model which generates y χ ≈ y ν allows for constraining the mass of the heavy neutrino since then Eq. (4.13) reads (4.14) Thus, to fit the observed DM density (4.4), M N ≈ 10 TeV is required. Since we are investigating the non-resonant regime we have M N < 2m DM . Therefore, we find a lower bound on the DM mass of m DM 5 TeV if we naively assume the behaviour for large DM masses to be also correct for parameters close to the transition of the non-resonant to resonant regime. We achieved this result by assuming n N = n eq N , m χ M N and by only taking into account the dominant processes of the SM particle scattering. From Eq. (4.12), we see that the contribution of the heavy neutrino scattering processes accounts for roughly eighty percent of the produced DM in case of y χ = y ν . Thus, the result will be altered significantly if the heavy neutrinos are out of equilibrium during the time where the production via heavy neutrino scattering is efficient. Also, we expect a significant change in areas of the parameter space where m χ ≈ M N , whereas taking into account the sub-dominant processes does not have a significant impact since they are suppressed by and only accessible after electroweak symmetry breaking. For these reasons, we solve the Boltzmann equations numerically for various coupling structures in Sect. 5.
Additionally, 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: 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.12) does not apply anymore. To illustrate that we look at the case y χ = y ν , where (4.12) 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 feebly 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 . 7 Consequently, the number density can be estimated by the equilibrium density at freeze-out:  restricts fermionic DM to have a mass of at least roughly a keV [22]. Therefore, DM must be bosonic in this case. However, this case also is in tension with observations of the Lyman-α forest which allows to probe structures of the size 10 0−2 h −1 Mpc [23]. This issue is treated in more detail within Sect. 6. We summarized our results for the case y χ = y ν in a schematic plot (see Fig. 3).
Since we investigate a feebly coupled sector, the back reactions in the DM production processes can be neglected. Only for the processes N ↔ νh responsible for producing the mediator N the back reactions are relevant, since for most of the parameter space N equilbrates with the SM. Therefore, we solve the Boltzmann equation in two steps: 1. The N production via N ↔ νh is solved at the level of the momentum distribution function, thereby taking into account the non-thermal shape of the distribution. The details of solving the Boltzmann equations at the level of momentum distribution functions are given in Appendix A and the collision term for the process in (T ). 8 We take n N n eq N (T → ∞) = 0 as initial condition. 2. This quantity is used to solve the Boltzmann equations for DM production via heavy neutrino and SM particle scattering employing the formalism described in Sect. 3. We take vanishing number densities for the DM particles as initial conditions. The SM particles are assumed to follow their equilibrium densities throughout the production process. The final result is then given by Y DM = Y χ + Y φ for T → 0. Note that the independent solution of the Boltzmann equations for the dark sector particles and the heavy neutrino is only possible due to the tiny interaction rate, which allows to neglect the back reactions from DM production via heavy neutrino scattering.
The results are summarized within Fig. 4. From our earlier considerations in Sect. 4.3 we expect the setup to work for a constant mediator mass M N as long as m χ M N . This constant value can be obtained by solving Eq. (4.13) for a given coupling structure. Consider e.g. the case y ν = y χ , where Eq. (4.13) results in M N ≈ 10 TeV. This case is illustrated by the solid blue line in Fig. 4. For 10 TeV ≤ m χ ≤ 10 4 TeV the prediction is met by the numerical solution. For larger DM masses, however, a larger mediator mass is required to accommodate the observed relic density. This is due to the following reason: The freeze-in mechanism produces DM efficiently down to temperatures around the heaviest mass  [24], corresponding to the lower (upper) bound on the sum of the three active neutrino masses. Since the production rate of N is proportional to the sum of the active neutrino masses this choice shows the earliest and latest point of equilibration. The heavy neutrinos reach equilibrium for T ≈ 10 3 M N in both cases. In Fig. 4, the lower bound on the sum of the neutrino masses was used, but the results are negligibly different when considering the upper bound involved in the production process. For the non-resonant regime this mass is given by the DM mass itself. Therefore, DM production is efficient for T m χ . The mediator mass and, thereby the neutrino Yukawa y ν , start to increase as soon as n N n eq N (T ) 1 for T m χ , since this suppresses DM production via heavy neutrino scattering. In case of y ν = y χ heavy neutrino scattering accounts for 35 41 of the DM production if the heavy neutrinos are following their equilibrium density during the time of production. If this contribution is missing, it has to be compensated by a larger neutrino Yukawa which results in a larger mediator mass.
The heavy neutrinos reach thermal equilibrium with the SM via the inverse decay ν i h → N j for T ∼ cM N . The factor c is independent of the neutrino Yukawa y ν and the parameters θ which encode the flavor structure of the neutrino Yukawas. However, it depends on the sum of the three active neutrino masses since the decay rate is proportional to this sum. The evolution of the heavy neutrino number density is shown in Fig. 5. Here, the heavy neutrinos reach equilibrium for T ≈ 10 −3 M N . Therefore, the lines in Fig. 4  For f 1 = 0.1 the contribution of SM particle scattering is suppressed by a factor of 10 since the contribution of the SM particle scattering is proportional to f 1 . Thus, a larger coupling compared to f 1 = 1 is required. This effect can be seen in Fig. 4 where all dotted lines lie above the solid line of the same color.
The different couplings structures result in larger (smaller) mediator masses for a small (large) dark Yukawa coupling compared to the neutrino Yukawa. Additionally, the effect of a small f 1 differs for a small (large) dark Yukawa. While the increase with a larger DM mass becomes less significant for a small dark Yukawa, the absolute difference between the small and large f 1 cases becomes stronger. This is due to the different contributions from heavy neutrino and SM particle scattering for the different coupling structures.
For smaller DM masses close to the transition to the resonant regime, the correct DM relic density is obtained for values of M N very close to M N = 2m χ . Although not visible within Fig. 4, all lines follow the black line down to small DM masses until the enhancement close to the resonance is not strong enough anymore to generate a sufficient amount of DM. However, the numerical solution is not trustworthy in this area due to numerical instabilities and therefore not presented here. We estimate the lower bound on m χ by evaluating Eq. (4.2) in the limit M N → 2m χ . In the case of y χ = αy ν we obtain m χ α − 4 3 MeV.

Constraints
In this section we discuss different constraints on the model. At first we discuss constraints from structure formation which pose strong limits in the resonant regime. Afterwards we investigate the impact of direct detection bounds on our parameter space and briefly discuss charged lepton flavor violation and indirect detection.

Structure formation
Since DM particles only interact weakly with the SM they can escape from gravitational wells formed in the early universe, thereby delaying structure formation below their freestreaming scale. Given the redshift at the production time z prod the free-streaming scale is given by where v (z) is the DM velocity at a given redshift z.
The observation of absorption lines in the spectra of distant quasars mostly induced by hydrogen clouds, the so called Lyman-α forest, allows for probing structures on the scale of roughly 10 0−2 h −1 Mpc [23].
Following the lines of [25], we estimate the free-streaming scale for the case of DM in equilibrium with the SM up to a certain freeze-out temperature and for the case of resonantly produced DM still in the freeze-in regime. Within this model, the first case applies to the resonant production with a coupling structure of y ν y χ whereas the latter is present in the resonant production regime for y χ y ν . The non-resonant production regime is not investigated here due to the much larger DM masses that are required to generate the observed relic density. Therefore, we do not expect this case to be in tension with the Lyman-α forest.
As it was pointed out in [26], the free-streaming scale should only be understood as an order-of-magnitude estimator in the case of non-thermal DM momentum distribution and may differ up to O (1) factors from results obtained with dedicated tools like the CLASS-code which computes the linear matter power spectrum.
For the purposes of this work, the estimation of the freestreaming length suffices, firstly because the non-thermal momentum distribution produced by the resonant freeze-in process (Eq. (A.21)) is close to a thermal shape and secondly because the resonantly produced DM for the freeze-out case will be excluded by this method by roughly two orders-ofmagnitude.
We approximate the velocity in Eq. (6.1) by the average velocity at the production time z prod which is only redshifted and p prod = dp p 3 f ( p, T prod ) dp p 2 f ( p, T prod ) . (6.4) Moreover, the Hubble parameter is given by For the numerical evaluation, we use the cosmological parameters of [27]. Lastly, we use the relation between the temperature and the redshift T = T 0 (1 + z) give T prod in terms of the redshift. The temperature T 0 refers to the temperature today. Inserting these expressions into (6.1) allows for calculating λ f s in terms of the production time z prod and the average momentum at this time p prod . Then, the result is compared to the upper bound on the freestreaming scale of λ f s 0.1 Mpc which was derived in [25] assuming that the particle species in question accounts for all of the observed DM relic density.
In case of resonant production with y ν y χ we can assume DM to have a Boltzmann like momentum distribution, i.e. f ( p, T ) = exp(−E p T −1 ). We take the time of production to be the freeze-out temperature since the interactions of DM with the SM cease to be efficient from this point on. For this distribution the average momentum results in By comparing the interaction rate of the process vh → χφ in the resonant regime to the Hubble parameter we find that T prod ≈ M N . For mediator masses M N MeV the freestreaming scale becomes insensitive to the mediator mass itself beside the change induced by the different g eff S (T prod ). In this case, we find a lower bound on the DM mass of m χ 10 keV. However, we found in Sect. 4.3 that a DM mass of 0.1 keV is required in order not to overproduce DM within this scenario. This lies two orders of magnitude below the estimated lower bound. Therefore, the resonant production regime with y ν y χ is excluded by the Lyman-α measurement.
If, on the other hand, y χ y ν , DM does not equilbrate with the SM even in the resonant production regime. Therefore the spectrum is non-thermal and given by Eq. (A.21). We take z prod T prod as the temperature where the derivative of the total particle number with respect to the time is maximized. Therewith, we find T prod = 3.36M N which results in p prod = 0.4T prod . Here, we also find that for M N m χ the free-streaming scale is insensitive to the mediator mass and the lower bound on the mass results in m χ 3 keV.
To summarize, the Lyman-α measurement strongly constraints the resonant production regime of this model. While the case where the resonant enhancement of the production cross section is strong enough to equilibrate DM with the SM is completely ruled out, the freeze-in regime is only allowed for couplings y χ 10 −12 M N keV with m χ 3 keV.

Direct detection
Direct detection experiments search for interactions of DM with nuclei. In this model, a coupling of DM to the Z boson is generated at one loop. The corresponding Feynman diagram is shown in Fig. 6. The coupling to the Z is then given by L ⊃ g Z χχχ γ μ P L χ Z μ with [28] g Z χχ = − y 2 where we have used the best fit values of [29] for the parameters of the PMNS matrix in case of normal ordering, which Therewith, DM interacts with quarks via Z exchange. Since this process happens at energies far below the Z mass, the heavy mediator is integrated out leading to where g qv and g qa are the couplings of the SM quarks to the Z . At low energies only the vector-vector and axial-axial interactions are not suppressed by powers of the relative velocity or momentum transfer, thereby leading to a spinindependent and a spin-dependent DM-nuclei cross section, respectively [30,31]. For the spin-independent cross section we obtain [31] with μ χ N = m χ M Xe m χ +M Xe , g uv = g w 1 4 cos θ w − 2 sin 2 θ w 3 cos θ w and g dv = g w − 1 4 cos θ w + sin 2 θ w 3 cos θ w . This cross section is constrained by the XENON experiment, as shown in Fig. 7. Therefore, the freeze-in setup cannot be constrained by this measurement. There are scenarios considered in the literature which allow for having a large direct detection signature even in a freeze-in scenario [32]. Fig. 7 The expected direct detection signals for the coupling structures investigated within Sect. 5 are compared to the current bounds from XENON1T [33] (dashed black curve). The dip in the red curve is due to a cancellation appearing in the loop function In [32], the cross section is enhanced due to a very light mediator. Since in our model the interaction is mediated by a Z boson this does not apply here.

Indirect detection and HEP phenomenology
Prospects for indirect detection of DM such as the observations of γ -rays from the galactic center or the precise measurement of the CMB all rely on the efficient annihilation of DM into SM particles. In the case of neutrino portal DM this usually happens subsequently by DM first annihilating into heavy neutrinos which then decay or annihilate into SM particles. Several prospects for indirect detection were investigated in [34] for the case of freeze-out production of DM where before DM freezes out its annihilation is efficient. This, however, is not the case for the freeze-in scenario investigated in this work. Here, the process is efficient only in the direction of DM production. This leads to a suppression of the annihilation cross section σ v which enters all observables of indirect detection considered in [34] since the couplings y ν and y χ are required to be feeble. Moreover, the annihilation rate is suppressed by a factor n DM n eq DM compared to the freezeout case. For this reason we do not study indirect detection observables within this work.
The minimal version of the type I seesaw mechanism employed here induces couplings of the SM gauge bosons and the Higgs to the heavy neutrino states. This can modify electroweak precision observables and induce charged lepton flavor violation (LFV) as well as additional Higgs decay channels in case of a light heavy neutrino [35,36]. The strongest constraints come from the decay μ → eγ with B(μ → eγ ) ≤ 4.2 · 10 −13 [24]. Within this setup the decay is mediated at one loop level by a W boson and a neutrino. The branching ratio of this process is then given by [37]: where F (x k ) is a loop function with x k = m 2 k M −2 W . Since we assumed the heavy neutrinos to be mass-degenerate and the light neutrino mass is tiny compared to m W we split the sum in the numerator into two parts with F (0) = 10 3 and F M 2 N M 2 W . Additionally we neglect the small deviation from one in the diagonal elements of U PMNS U † P M N S in the denominator. Since the mixing matrix U is unitary we find Taking the best fit values from [29] we find 12. Thus, we can give the branching ratio as a function of the heavy neutrino mass only since the free parameters of the orthogonal matrix R cancel within this setup [15]. This expression is maximized for M N = 1.36M W and results in which is far below the experimental limit. For this reason, we also expect other LFV and electroweak precision observables not to significantly constrain the scenario.
Another imprint of this model could be found in additional decay channels of the Higgs if M N < m h . In this case the decays h → ν i N j and h → N i N j are kinematically accessible. As pointed out in [10,38] the dominant contribution comes from the decay into a heavy and a light neutrino. However, branching ratios of this process larger than 10 −2 are already ruled out and are typically much smaller due to the tiny Yukawa coupling [38]. Therefore, the contribution is negligible.

Conclusion
We have investigated a minimal neutrino portal DM model. The SM is extended by three right-handed neutrinos which generate the neutrino masses via a type I seesaw mechanism and, furthermore, act as mediator between the SM and DM. The dark sector consists of a boson φ and fermion χ coupled to the right handed neutrino via a Yukawa coupling. 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.
We derived analytic solutions for the number density in the resonant (M N > m χ + m φ ) and non-resonant (M N < m χ + m φ ) DM production regime. Adding the requirement that the coupling of the right-handed neutrino to the SM is of the same order of magnitude as its coupling to the dark sector allows for the prediction of the mediator or the DM mass respectively. In the non-resonant regime, we find M N ≈ 10 TeV. The non-resonant regime is studied in more detail numerically, as seen in Fig. 4.
Within the resonant regime, however, for y χ y ν the resonant production of DM is strong enough to bring it into equilibrium with the SM. Thus, the freeze-out mechanism is revovered although the couplings between DM and the SM are feeble. 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 is required to be y χ ≈ 10 −12 M N m χ . The resonant scenario is strongly constrained by the measurement of the Lyman-α forest. The freeze-out case can be excluded completely, while freeze-in with y χ ≈ 10 −12 M N m χ is only viable for m χ 3 keV. Charged lepton flavor violation, Higgs decays, indirect detection and direct detection have little impact on our parameter space due to the feeble coupling of the SM to the dark sector. Thus, producing the observed DM energy density within this model of neutrino portal DM is possible even with small couplings between the SM and the dark sector.
Although within this work CP violation in the PMNS matrix was assumed to be absent, it could be included in the analysis to explore its phenomenological imprints and its impact on leptogenesis. gration by assuming that if a particle distribution deviates from its equilibrium density it differs only by a momentumindependent factor, i.e. f i = α i f th i with ∂α i ∂ p i = 0. Furthermore, the equilibrium densities of bosons and fermions are approximated by a Boltzmann distribution.
Following the lines of [26,39] we solve the Boltzmann equations at the level of momentum distribution functions. This has the advantage of a more accurate solution and the exact shape of the momentum distribution allows for more insights into the process of structure formation. Throughout the calculation we approximate the equilibrium densities of any particle species by a Boltzmann distribution. The Boltzmann equation is given by: Here t is the time, H the Hubble parameter, f is the momentum distribution function of the particle species whose evolution is described by this Boltzmann equation, p is their momentum and C( p, T ) is the collision term which describes the impact of interactions. For the integration of this equation it is convenient to perform a coordinate transformation (t, p) → (r, x) such that the differential operator on the left hand side contains a derivative with respect to one of the new variables only. If r only depends on t and The condition (A.2) is fulfilled if For the last equality we used the conservation of entropy s(T 0 )a(T 0 ) = s(T )a(T ) = const. and g s are the entropy degrees of freedom. The conservation of entropy also allows us to relate the temperature T to the time t: Since T is only a function of t and not of p we can choose with m 0 being am arbitrary mass scale. Combining all this the Boltzmann equation results in Since in this work DM production is mainly governed by 2 ↔ 2 scattering processes we will discuss the collision term for these type of processes in more detail. For a A + B → C + DM scattering the collision term for the evolution of the momentum distribution function of DM is given by: Here, E i = p 2 i + m 2 i , M is the matrix element for the process A + B → C + DM which is the same in both directions since we are assuming CP invariant interactions and f i is the distribution function of particle species i. We assume that f C f DM f A f B which is justified since the paper explores the freeze in production of DM. Furthermore, we take f A/B = f th A/B assuming the interactions of A and B are efficient enough to keep them in thermal equilibrium. Moreover, taking f th A/B to be a Boltzmann distribution, shifting the integration over p c to p C +p DM = P and multiplying the equation The equation above can be simplified by rewriting it in terms of the reduced cross section [40]: Moreover, we change the variables of integration from d 4 P to an integration over the zero component of the center of mass momentum vector P 0 , the center of mass energy s and the angle θ between center of mass momentum P and the momentum of the DM candidate p DM , d 4 P = 2π P 2 d P 0 dPd cos(θ ) = 2π P 2 0 − sd P 0 dsd cos(θ ). To eliminate the remaining δ function we express the argument in terms of cos(θ ): where cos(θ 0 ) is the value required for cos(θ ) for a vanishing argument of the δ function. Therewith, Eq. (A.10) results in The last integral basically restricts the boundaries of either P 0 or s in the sense that if is fulfilled | cos(θ 0 )| ≤ 1 must hold. This requirement yields the inequality In case of m C = m DM 9 this results in a lower (relative minus sign) and upper bound (relative plus sign) of the P 0 integration of The last equality is given to showcase that in case of m DM = 0 only a lower bound exists, as was shown in [39], while for finite DM masses there is also an upper bound. Thus, we have The s integral and the following integration of the differential equation for an arbitrary cross section cannot be performed analytically. However, in case of a very light DM candidate (m DM ≈ 0) and a resonant production process with mediator M mediator the integral can be evaluated analytically. Moreover, this case is of special interest for this work since for resonant production the DM mass turns out to be below keV. Therefore, the exact shape of the momentum distribution is required to quantify the impact of DM on structure formation. In this case we have P + 0 → ∞ and In the last step, we assumed that the temperature where we observe the DM density is much smaller than the mass of the resonant particle. As mentioned above, to derive this analytic result we took the effective entropy degrees of freedom to be a constant. Hence the above formula is only a good approximation as long as we take T large enough to stay at a constant value of g s (T ) ≈ 100. Of course, we observe the universe at a smaller temperature. However, the above result remains a good approximation if the main part of the production has been finished before g s (T ) starts to vary significantly since for a collisionless particle species the quantity Y = n s is a constant.
By comparing the number of produced DM particles at temperature T to the number of particles for T → 0, n(T )T 3 lim T →0 n(T )T 3 , with an unapproximated n (T ) we find that for T ≈ M N 4 already over 0.99 of DM particle have been produced. Thus, as long as M N ≥ 100 GeV the result (A.22) serves as a good estimate.
Beside collision terms for 2 ↔ 2 scattering processes, the collision term for the (inverse) decay N ↔ νh is required. The procedure for performing the integration over the particle momenta follows the same lines as for the 2 ↔ 2 scattering. Thus, we only give the result for the collision term resulting from the decay that appears in the Boltzmann equation for the heavy neutrino N : (A.23)

Appendix B: Cross sections
Here, we give the relevant reduced cross sections for the case m φ = m χ . Since CP conservation is assumed the reduced cross sections for a process and its time reserved process are the same. 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 > 2m χ . The decay width is given by: