Thermal effects in freeze-in neutrino dark mater production

We present a detailed study of the production of dark matter in the form of a sterile neutrino via freeze-in from decays of heavy right-handed neutrinos. Our treatment accounts for thermal effects in the effective couplings, generated via neutrino mixing, of the new heavy neutrinos with the Standard Model gauge and Higgs bosons and can be applied to several low-energy fermion seesaw scenarios featuring heavy neutrinos in thermal equilibrium with the primordial plasma. We find that the production of dark matter is not as suppressed as to what is found when considering only Standard Model gauge interactions. Our study shows that the freeze-in dark matter production could be efficient.


Introduction
The most significant experimental and observational evidences for physics beyond the Standard Model (SM) are the lack of a dark matter (DM) candidate and the unknown origin of neutrino masses (and of the lepton mixing).The existence of sterile neutrinos allows to accommodate both these puzzles in a rather straightforward way.Heavy neutral leptons (like right-handed neutrinos), singlets under the Standard Model gauge group, are the most popular option for the generation of the light neutrino masses via different variants of the seesaw mechanism.In this work, we will focus on low-scale seesaw scenarios, i.e. we will assume that the heavy neutral leptons (HNL) responsible for neutrino mass genertion have masses at a scale not exceeding few TeVs.Within this framework it is, moreover, possible to accommodate the masses and couplings of some of these new states to make them stable on cosmological scales and, consequently, be potential DM candidates.Such requirement impacts, in particular, the mass scale of the sterile neutrino DM, favoured to be O(1 − 100) keV.
These HNLs mix with the active neutrinos via Yukawa couplings to the SM Higgs boson.In the case of the DM, such mixing is responsible of the production process dubbed Dodelson-Widrow mechanism (DW) [1].Although this mechanism works for a wide range of masses, it is severely constrained by a variety of observational evidences.First of all, X-ray searches for the decay of the sterile neutrino to an active one and a photon constrain, as function of the mass, the coupling of the sterile neutrinos to the active ones and provide an upper bound, of the order of tens of keV, to the mass of the DM candidate [2].Very light masses, up to few keV, are instead ruled out by the so called Tremaine-Gunn bound and by structure formation consideration as, for example, phase-space density [3,4] and the number of satellites of the Milky-Way [5][6][7].For masses above 2 keV, another constraint holds from Lyman-α forest data analysis (from which it is possible to infer the spectrum of matter density fluctuations, determined by the DM properties).However, this constraint is strongly model dependent and the bounds are related to the warm dark matter (WDM) production mechanism, and to which extent this mechanism contributes to the total DM abundance.It was shown in [8] that, once all these constraints are accounted for, the DW mechanism can only account for at most ∼ 50% of the total DM abundance.A possible way out would be guaranteed by the presence of a lepton asymmetry (in such a case, the DM production mechanism is dubbed Shi-Fuller [9]).This would lead anyway to a very constrained parameter space [10,11] unless further ingredients are added to the theory.
Another economical production mechanism for DM would be provided by the decay, à-la-freezein [12][13][14][15][16], of heavy right-handed neutrinos.The effectiveness of such mechanism was studied in [8], in the framework of a neutrino model based on the inverse seesaw (ISS) [17], as well as within the Type-I seesaw mechanism [18].In the ISS framework the spectrum of the SM in enlarged by the inclusion of right-handed (i.e.featuring in the interaction basis Yukawa coupling with the SM Higgs) and sterile (i.e.no Yukawa couplings with the Higgs present in the interaction basis) heavy neutrino states.The mass spectrum of the theory features the three light massive neutrinos, and a series of pairs of heavy pseudo-Dirac neutrinos.In the case in which the number of sterile neutrino fields is greater than the number of right-handed ones, additional physical states, with masses lying at intermediate values between the light neutrinos and heavy pseudo-Dirac pairs, are present [17].The most minimal model of this type, complying with neutrino oscillation data, is the so called ISS(2,3) (2 right-handed neutrinos and 3 sterile neutrinos).In this case, besides the two pairs of pseudo-Dirac heavy states, an intermediate state (with mass m DM ∼ O(keV)) arises, above the light neutrino masses, providing a good DM candidate [17].
In this work we will revisit the studies in Ref. [8,18] on the possibility that a sterile neutrino produced from the decay of heavier HNL accounts, at least partially, for the dark matter component of the Universe.In the ISS mechanism, the Yukawa couplings are generally large enough such that the HNL states thermalise in the early Universe, while the light sterile state does not; this enables the freeze-in production of DM via the decay of a heavy HNL state into it and the Higgs boson.The freeze-in mechanism ensures an efficient DM production for a sufficiently suppressed mixing with the active states to evade observational constraints [18,19].However, as recently pointed out in [20], following a seminal work presented in Ref. [21], the mixing angle between active and sterile neutrinos is subject, in the early Universe, to potentially large thermal corrections.Thermal masses are indeed relevant at the time of DM production through freeze-in, and dominate over the keV mass of the neutrino DM considered here, such that the effective light neutrinos propagating in the medium are found to be approximately the flavour ones.This renders, for example, freeze-in production via W and Z decays not capable of accounting for the correct amount of DM [21] 1 .
We extend the study conducted in [21] to realistic seesaw scenarios, including, in addition to the DM state, additional heavier RH neutrinos (with masses of the order of a few hundreds GeV) with large enough Yukawa interactions to be in thermal equilibrium.These models can generally account for the massive nature of active neutrinos and lepton flavour mixing.These RH neutrinos can then decay into the keV DM and a Higgs boson, typically at temperatures of the order of their masses.Our present work also differs form the one in Ref. [21] by the fact that all the sterile neutrinos, including the DM state, are of Majorana nature.This scenario will ease the embedding of a lowscale fermionic seesaw rendering the framework more realistic and ultimately UV complete.In our analysis, we will consider both the type-I and inverse seesaw mechanisms as benchmark scenarios to illustrate our findings.We focus on production through active-heavy neutrino mixing, thus studying temperatures below the electroweak (EW) phase transition, T SSB ∼ 160 GeV [23], which we take in the following to be instantaneous.Moreover, in order to be able to make a comparison with previous estimations of the keV neutrino production rate, we only take into account the production through two-body decays.Note however that potentially important contributions to the production rate can arise from 2 → 2 scatterings as well as from the Landau-Pomeranchuk-Migdal (LPM) effect at high temperatures [24][25][26], but we leave the inclusion of such effects, and thus the study of the production rate above the EW crossover, for future work.
This work is organised as follows: the framework and the dark matter production are discussed in Section 2. Section 3 details the formalism we use and the strategy to evaluate thermal corrections taking into account the Higgs and gauge boson contributions to the neutrino self energy for a general case.We show in Section 3.5 how we compare in the case with one light neutrino and one sterile state with the results of earlier work in [21], finding qualitatively similar conclusions.The analysis for dark matter production below the EW crossover is thoroughly conducted in Section 4, first without HNL, and then in the realistic case with two heavy neutrinos in addition to the dark matter candidate.Our findings are summarised in Section 5. Further technical details of the computations are collected in Appendix A and B, while we devote Appendix C to show DM production in an example with a larger dark sector including a new scalar, not suffering from the aformentioned mixing suppression, for completeness.
2 Common origin for neutrino masses and DM

Neutrino masses from approximate lepton number symmetry
Among the several ways to explain the lightness of neutrino masses, the type-I seesaw mechanism [27][28][29][30][31][32][33] provides a simple one by introducing a number of RH neutrinos2 .On general grounds, given that these RH neutrinos are complete singlets under the SM gauge group, their introduction allows not only for the usual Yukawa coupling to the SM Higgs boson, as for the other fermions of the SM, but also for a Majorana mass term, M N , which breaks lepton number by two units.The SM Lagrangian would thus be enlarged with the following lagrangian where Φ ( Φ = iσ 2 Φ * ) and L L are the Higgs and lepton doublets, respectively, and N R the RH neutrinos.Both the Majorana mass M N and the Yukawa couplings Y ν are matrices in flavour space.
After spontaneous symmetry breaking (SSB), the Higgs develops a vacuum expectation value, v H , and generates a Dirac mass term between SM (active) neutrinos and the RH ones, We can work, without loss of generality, in a basis in which the matrix M N is real and diagonal.In the seesaw limit, where the scale of the Majorana mass terms is much larger than the Dirac ones (so that it is possible to define the expansion parameter Θ = m D M −1 N ), the light neutrino mass matrix results at leading order from the expression In this case, the active-heavy neutrino mixing, which would change leptonic weak interactions, is given at leading order by θ ≃ m * D M −1 N and is in general suppressed due to the relative high-scale of M N .The phenomenological impact of the RH neutrinos can be important if they are not excessively heavy, and have sizeable mixings θ.In the so-called vanilla seesaw, where no special relation is present in the submatrices M N and Y ν , the active-sterile mixing is strongly suppressed by the ratio of mass scales between active and sterile neutrinos, |θ| ≲ m ν /M N ≲ 10 −5 GeV/M N resulting, for Majorana masses above the GeV scale, to mixings that are way beyond the reach of current experiments.On the other hand it is important to notice that, being Eq.( 2) a matrix relation, it is perfectly reasonable to make it reproduce current data on neutrino masses and mixing while having an active-sterile mixing that does not follow the vanilla scaling.This may be due to accidental cancellations in the matrix equation or, more interestingly, to an underlying symmetry of the theory.In particular, there exists a one-to-one correspondence between a global lepton number symmetry and massless active neutrinos in the SM extended with fermionic singlets [35] implying that the mass scale of active neutrinos is suppressed (and protected from large radiative corrections) if an approximate lepton number symmetry is present in the theory, even in the presence of sizeable Yukawa couplings.Such fermionic singlet extensions of the SM can feature new physics signals strong enough to be testable in current experiments, while correctly reproducing neutrino oscillation data [36].
Within the type-I seesaw, a convenient parametrisation to effectively explore the parameter space of the model is the Casas-Ibarra (CI) one [37], in which the Dirac mass matrix depends on the lowenergy neutrino oscillation parameters as where U PMNS is the lepton mixing matrix, m a diagonal matrix containing the masses of the light active neutrinos, M an analogous matrix for the heavy neutrino masses and R is an orthogonal matrix which can be parametrised by three complex angles3 ω ij , R = V 23 V 13 V 12 , where The CI parametrisation in Eq. (3) guarantees to recover the input values of U PMNS and m from the diagonalisation of the full Lagrangian in Eq. ( 2), as long as the seesaw hypothesis, Θ = m D M −1 N ≪ 1, is satisfied.This means that the imaginary parts of the angles ω ij in Eq. ( 4) cannot be too large, as the Yukawa couplings grow exponentially with them, although it is evident that the specific upper bound is also a function of the matrix M .In the framework of the CI parametrisation, a lepton-number conserving scenario is recovered in the limit ω 12 , ω 23 → 0, ω 23 → i∞, M 1 → 0 and M 2 → M 3 : in this case the mass spectrum of the theory features two massless particles and two massive Dirac states, thus mediating only lepton number-conserving interactions.Such exact limit clearly does not reproduce the observed neutrino data, while a perturbed scenario with an approximate lepton number symmetry (like in for instance the case where ) is phenomenologically viable, while featuring a relatively low new-physics scale and sizeable active-sterile mixings (in this case the mass spectrum contains the three light active neutrinos, a light Majorana state with mass M 1 and with a feeble mixing with the active neutrinos, and a heavy pseudo-Dirac pair with mass M 2,3 ).Other lepton number symmetric-inspired constructions of the type-I seesaw are possible, cf.e.g.[38].
In some variants of the type-I seesaw realised at low scale, other singlet (from SM gauge interactions) neutrinos (ν S ) are considered, as in the case for the ν-MSM [39,40], the Inverse [32,41,42] (ISS) and Linear [43,44] (LSS) seesaw mechanisms; these variants allow to have large neutrino Yukawa couplings with a comparatively low seesaw scale.We will thus focus on mechanisms that rely on an approximate lepton number symmetry rather than on a hierarchy of scales to explain the light active neutrino masses.Beside their phenomenological impact, these simple extensions of the SM may also accommodate the observed DM relic abundance, as in most successful realisations, the lightest sterile state can account for the DM, while the other heavier states are employed to generate the masses of the active neutrinos at tree-level.
In a minimal setup4 , we introduce two sets of RH neutrinos, N R and ν S , respectively, which only differ in that they have opposite lepton number assignments, L = 1 and L = −1.The most general Lagrangian (regarding mass and interaction with the sterile neutrinos) which respects lepton number conservation is then where Y ν and M are Yukawa and Dirac mass matrices, respectively.After SSB the mass matrix can be arranged by blocks as5 At this level, light neutrino masses are exactly zero at all orders, thanks to the lepton number symmetry, while the heavy states form Dirac pairs with masses of order m N ∼ O M 2 + m 2 D .The active-heavy neutrino mixing, θ = m * D M −1 , is not necessarily suppressed (considering M at a scale between tens and hundreds of GeV).Only once small lepton number violating terms are introduced in the Lagrangian, light neutrinos will acquire a mass.The possible LNV terms that can be considered to complete L LC are given below (all these terms break lepton number in two units, ∆L = 2) which translate into the following neutrino mass matrix: Under the approximate lepton number symmetry, such that µ, µ R , ϵm , we can safely neglect the effect of µ R , which contributes to light neutrino masses at loop level 6 .Again, while in principle there is no reason for lepton number to be approximately conserved, the inclusion of these small LNV terms is still natural in the 't Hooft sense, as lepton number protects them from receiving large loop corrections.Depending on which sources of LNV are included in the Lagrangian, one can distinguish among different scenarios, namely the linear seesaw (when only ∆M LSS is included), the inverse seesaw (including only the small Majorana mass, µ, for ν S ), or, what we name the LISS scenario, when both µ and ϵm ′ D contributions are included.Upon the inclusion of these LNV terms, light neutrinos acquire masses proportional to the sources of LNV, while heavy neutrinos now form pseudo-Dirac pairs whose mass splitting is also proportional to the LNV terms.In the case of the LSS one finds for light neutrino masses while for the ISS case, we have where we have defined the mass splitting of the pseudo-Dirac pairs as ∆m heavy .
In the context of the ISS, it was found that at least 2 pairs of N R and ν S were necessary to explain oscillation data [17], while the addition of an extra ν S field to the particle content would additionally provide a good DM candidate at the scale O(µ).The resulting model was hence labelled as ISS(2,3) [8].Indeed, for heavy states at M ∼ TeV, µ has to lie at the keV-scale to explain light neutrino masses, and thus the decoupled ν S field becomes a good warm DM candidate in the few-keV range whose interactions with the SM are suppressed by powers of O (m ν /µ, µ/M ) with respect to those of the SM neutrinos.On the other hand, within a type-I seesaw construction, only two N R are necessary to explain neutrino oscillations, while a third, lighter N R with suppressed Yukawa couplings can be added and identified as the DM candidate.
In the following, we will start by considering the "single family" case, in which we study only one SM neutrino, the DM candidate, and one Dirac pair of heavy neutrinos.This will already capture some of the most important features we will consider for the DM production, namely the existence of a heavy neutrino species capable of decaying into the SM bosons and the DM neutrino, as well as the relevant gauge interactions with the plasma of each neutrino species.This extends over previous studies where only the DM candidate and a single SM neutrino were considered [21,22], opening the window to new production channels for the DM.Later on, we will consider full models that completely account for oscillation data as well as including the keV neutrino DM candidate, and point out regimes that cannot be described in the single family approximation.

Minimal framework
As previously discussed, the main ingredients of the framework under scrutiny are represented by the presence of at least one heavy decaying neutrino, n h , capable of existing in thermal equilibrium in the early stages of the history of the Universe, as well as a sterile neutrino DM candidate, n DM .In order to ensure its cosmological stability, its coupling with the SM states should be suppressed, which motivates the assumption that it never existed in thermal equilibrium in the early Universe and its abundance is negligible prior to the freeze-in process.
In the following, although we develop the formalism to compute and understand DM production from decays on full generality, accounting for any possible neutrino spectra and mixing between the different species, we will focus on two possibilities, motivated by Refs.[8,18].For very low-scale seesaw realisations (such as the so called ν-MSM, where the new states lie at the GeV scale), one generally expects production rates from decays to be suppressed because of the large mass difference between the parent particle and the Higgs boson, making the decay possible only for particles at the very end of the thermal distribution 7 , and thus we will consider the heavy neutrinos to lie around the EW scale.
We will first consider a simplified scenario motivated by the ISS(2,3) model proposed and studied in Ref. [8], which in the single family limit presents a simplified neutrino spectrum with just one light neutrino, n l , the keV DM, n DM , and two pseudo-Dirac heavy neutrinos, n h with h = 1, 2. Neglecting contributions from µ R in Eq. ( 8), the mass matrix in the basis (ν L , N c R , ν c S , ν c DM ) T is given by where all elements are now just mass parameters.Upon diagonalisation, one finds that the active-DM mixing, U α4 , and the active-heavy mixing, U α5( 6) are given by From here it is clear that, in order to comply with the strong X-ray bounds on the active-DM mixing, a hierarchy must always be present either between m D and M to make θ ≪ 1, between µ SS ′ and m DM , or both.This in turn would translate into a suppression of DM production rates.Finally, the mass gap between the pseudo-Dirac pair would be of order µ SS .Note that this suppression of the couplings are generally to be expected in any neutrino mass model when working in the single family approximation.
Second, we will work and compare with full neutrino mass models accounting for oscillation data.In particular we will work with the type-I seesaw, where only two RH neutrinos need to be introduced in the minimal case to explain oscillations, on top of the DM candidate.This in turn translates into having one massless active neutrino.To consistently explain oscillation data when studying the parameter space, we adopt the Casas-Ibarra parametrisation introduced previously in Eq. ( 3), working in the approximate lepton number symmetric limit.In contrast to the single family approximation, it is possible to have large couplings while complying with experimental and astrophysical bounds thanks to cancellations in the matrix equations.
Both in the single family approximation as well as with complete frameworks, we will in general be able to write the relation between the flavour and mass eigenstates as where U is the total neutrino mixing matrix whose dimensionality depends on the number of states.
We have separated the mass eigenstates by blocks, with light neutrinos, n l , the DM candidate with mass m DM , n DM , and the HNL, n h , whose masses are m N8 .The ISS can be readily accommodated by changing

DM production
The most popular option for the production of keV sterile neutrino DM is represented by oscillation processes in which active neutrinos are converted into DM in the early Universe.We recall that this kind of production mechanism, in its minimal incarnation, is the Dodelson-Widrow mechanism [1].DM production via oscillation gets enhanced when occurring in an environment with a sizable lepton asymmetry, case known as the Shi-Fuller (SF) production mechanism [9].In both scenarios, most of the DM is produced at relatively low temperatures, close to the QCD phase transition.Although thermal corrections are relevant also for these mechanisms, the computations discussed in this work do not strictly apply since they are aimed for much higher temperature regimes.Furthermore, the DW and SF mechanisms are nowadays very strongly constrained by both DM Indirect Detection (see e.g.[5] for a review) and structure formation [45,46].Nonetheless, in the following, even if the DW mechanism cannot account for all the observed DM, it will inevitably contribute to the final abundance, as it represents an unavoidable source of DM at lower temperatures.On the other hand, we will not take into account the SF production as it relies on a lepton asymmetry, which we will assume to be negligible here.
Our reference scenario will be, instead, represented by the case of DM production at earlier stages of the history of the Universe (i.e. at higher temperature regimes) via decays of heavy SM states (like the W/Z bosons) or BSM states, namely heavy sterile neutrinos.A very natural option, along this line, was proposed in [8].There, keV scale DM was embedded into the ISS(2,3) mechanism for the generation of the SM neutrino masses.This framework allows for the existence of O(0.1 − 1) TeV heavy neutrinos with sufficiently large Yukawa couplings to be in thermal equilibrium in the early Universe.These HNL can efficiently produce DM à-la freeze-in via their decay into a SM Higgs and a sterile neutrino.The DM relic density can be estimated as [13-16, 47, 48] where an eventual sum over multiple heavy neutrino states is implicitly assumed in the first term while the second term represents the sum for decaying SM bosons, B = H, W, and Z, into DM.A back of the envelope estimation from Eq. ( 14) leads to the fact that, neglecting production through SM boson decays, for a DM mass m DM ≃ 10 keV (typical scale for the ISS(2,3)) and m N ≃ 150 GeV for the HNL mass, a decay rate Γ(n h → n DM + SM) ∼ O(10 −16 ) GeV is needed to generate the correct amount of DM.
In the scenario considered in Ref. [8], the DM production process considered was the decay of heavy pseudo-Dirac neutrinos into the DM itself and a Higgs boson, finding that one could indeed simultaneously explain the DM relic density with a keV sterile neutrino and the neutrino oscillation phenomenon.This conclusion was then further studied in the context of the type-I seesaw in Ref. [18], while the production through weak-gauge boson decays into DM was investigated in Ref. [22], however neglecting possible thermal effects that arise at temperatures around the EW scale [21].
Other processes, relevant at higher temperatures, would be for example 2 → 2 scatterings and the so-called LPM effect.The former arises, for example, through scatterings of gauge bosons with leptons producing a Higgs and the DM candidate in the final state.The rate can be estimated, at high temperatures, as , with i the matrix element corresponding to the DM candidate.The latter arises from the resummation of multiple soft scatterings which become relevant (and can even be the dominant effect) at temperatures above the EW symmetry breaking.Thus, we expect its effect to be mild on the temperature ranges under consideration.Nonetheless, we stress that a proper inclusion of these effects is necessary in order to study the complete DM production at temperatures above T SSB ∼ 160 GeV.

Decay channels for DM production
The spectrum arising in the neutrino sector from simultaneously having the dark matter candidate and explaining oscillation data through heavier states, here assumed to lie around the EW scale, translates into a rich array of possibilities for DM production.
First, the mixing between active neutrinos and DM allows for the production through weak boson decays into a DM candidate and a lepton (either a charged lepton for W decays or a mostly active neutrino for Z decays9 ).This contribution was studied in Ref. [22] neglecting thermal effects, while the inclusion of thermal effects was first done in Ref. [21] in the context of Dirac neutrinos.In the following we will assume that neutrinos are Majorana particles, which quantitatively change the picture.Additionally, one could consider the Higgs boson decays into DM and a light neutrino, but given their rather light masses, this contribution is negligible as the coupling goes through the mass.
When heavy neutrinos, n h , have large enough mixings with the active ones, they can be present in the plasma at electroweak temperatures, opening up the possibility for various new decay channels not present when considering only the DM candidate.Assuming they are heavier than the weak gauge bosons and the Higgs, DM can be produced through the decay of n h into a Z boson and DM.Moreover, the rather large mixings necessary for the heavy neutrinos to be in thermal equilibrium translate into non-negligible couplings with the Higgs boson, making the production channel n h → H +n DM relevant.What is more, from the statistical distribution of the final states, the decay into the Higgs boson could potentially show an enhancement when compared to the decay of a boson into a couple of fermions.This is thanks to the fact that the decay to a scalar can have Bose enhancement, while a pair of fermions in the final state tends to suppress the rate with respect to the zero temperature case due to Pauli blocking [49].
Given the richness in the possible decays to DM, in the following we will study several scenarios which introduce successively a new contribution to the production rate.We thus consider the following cases, that are summarized in Table 1 with their production channels: • Considering solely one species of active neutrinos and the DM candidate, without the introduction of HNL.This is the simplest scenario, which would correspond as well to the case studied originally in the proposal of the Dodelson-Widrow mechanism [1].In this case the only relevant parameters are the active-DM mixing and the DM mass.For the production through decays, only weak-gauge bosons could produce a non-negligible amount of DM, given that the Higgs couples to the rather light masses.Although previously studied in Ref. [21] in the context of Dirac neutrinos, here we will instead consider them to be Majorana.When referring to this simple case, we will call it the "toy 2 × 2" scenario, given that it could describe DM but not account for neutrino masses and mixings.
• Introducing HNL in an ISS-like construction within the single family approximation.This is inspired by Ref. [8] where it was found that the ISS(2,3) could simultaneously explain oscillation data and the DM abundance through a keV-scale candidate.In our simplified one-family ISS, we will thus have one light neutrino, the DM candidate, and a pair of almost degenerate HNL, and thus refer to it as "ISS-4 × 4".We choose the parameters to comply with current constraints on mixing between the different species both for the DM and the heavy species, but point out that this case does not explicitly reproduce oscillation data, although there is in principle enough freedom when going to the full case to fit data.In any case it proves useful because it already introduces the non-negligible coupling between the Higgs boson and the heavy neutrinos, as well as the possible decay of heavy neutrinos into a gauge boson and the DM neutrino.
• Minimal type-I seesaw, fully reproducing oscillation data through the Casas-Ibarra parametrisation [50] from Eq. ( 3), with two almost degenerate HNL and the DM candidate.We choose the type-I as the coupling between DM, HNL and the Higgs is parametrically different in terms of Yukawa couplings to the one of the ISS.Additionally, we will first neglect the Higgs contribution to fully study how the addition of HNL changes the production compared to the "toy 2 × 2" scenario, and finally introduce also the Higgs contribution, which is non-negligible, and compare with the results from the "ISS-4 × 4".

Scenario Decay channels Toy
Table 1: Summary of the possible DM production channels through two-body decays available in the different cases we consider.Note that we assume here that the HNL (n h ) are heavier than the Higgs.In the opposite case, one could still have the decay H → n h + n DM .The difference between the "ISS 4 × 4" and the full type-I seesaw is that in the latter the matrix structure of the couplings allows for regimes in which cancellations between elements can allow for large couplings.
3 Production rates in thermal field theory

Formalism and strategy
We will make use of the real-time formalism of thermal QFT [51,52] in order to compute decay rates and dispersion relations for the collective propagating modes in the plasma [49,53,54].Our aim is to generalise the results of Ref. [21] to any neutrino mass model which includes SM-singlet fermions, without any assumption on the size of the mixing matrix elements, U αi , or the neutrino masses.This will allow to understand whether the presence of heavier neutrinos mixing with the DM can change the picture described in Ref. [21] where only a light neutrino and the DM were considered.
In the real-time formalism we have a doubling of the degrees of freedom, such as the free propagator becomes a 2 × 2 matrix.For a particle of mass m, taking the symmetric-Keldysh contour and in the absence of a chemical potential, its components can be written as [53] where f (k 0 ) is the equilibrium distribution function given by with η = ±1 for bosons and fermions, respectively, and β is the inverse of the temperature of the plasma.The factor d αβ encodes all the Lorentz structure of the propagator, and is given by for fermions, −η µν + kµkν m 2 , for vector bosons.(17) Even if one has to deal with a matrix structure for the propagator, there are some relations between the components of the propagators which can help simplify computations.In particular, the following relation can be found between the retarded self-energy of a fermion, σ R , and the (+−)-component of the self-energy in the real-time formalism [53]: It is indeed the imaginary part of the retarded self-energy that will eventually be related to the dumping rates of quasi-particle excitations in the plasma [53].Regarding the real part of σ R , it can be obtained through the following dispersion relation [21] Reσ As will be apparent in the next section, in general we will be able to write the self-energy contributions for massive neutrinos as where A ij are some combination of coupling constants appearing at the Lagrangian level, and the index k runs over mass eigenstates.Thus, we can still make use of the dispersion relation for σ k in order to find the real and imaginary parts of Σ.Given that the couplings A ij do not depend on momenta, in the following we can use the relation from Eq. ( 18) including them in Eq. (20).When A ij ∈ R, this quantity will be directly ImΣ ij , but this will not necessarily be the case.Thus, in the following we will dub the product which is the final quantity we are interested in.
On general grounds, in the following, we will decompose the self-energy as (suppressing generation indices to ease the notation) with p an unit vector in the direction of the momentum ⃗ p. Thanks to this decomposition, we will be able to project Σ and work with the scalar functions Σ (I) instead, obtained as When necessary, we will also take into account the chirality projectors P L(R) .As shown in Ref. [55], once the one-loop self-energy corrections are computed, the rate at which a species in the plasma reaches equilibrium will be given by the imaginary part of the self-energy.This is precisely the quantity we are interested in within the freeze-in framework.The advantage of using the one-loop corrections is that it directly includes all possible interactions with the plasma, as well as the correct dependence with temperature and momenta for the production rates.
Figure 1: Feynman diagrams contributing to the neutrino self-energy.In the SM extension with singlet fermions, there is the Higgs contribution (left) and the gauge (W and Z bosons) contribution from the W and Z bosons (right).The imaginary part of the self-energy will be related to the rate at which each species would reach equilibrium.

Higgs contribution
We study here the thermal self-energy contribution from the Higgs boson coupling, which after SSB can be written as [56] where C ij ≡ α=e,µ,τ U † iα U αj and m i(j) are the different neutrino masses and the second line uses the fact that neutrinos are Majorana.Using the propagators from Eq. ( 15), we obtain iΣ where r µ ≡ (p − q) µ and f B(F ) corresponds to the distribution function for bosons (fermions).To ease the notation, we have defined . From this expression we readily obtain which can be furthered simplified by decomposing the Dirac delta functions as with 10 ω k ≡ ⃗ q 2 − m 2 k , leading to Notice that the latter expression agrees with Ref. [57] in the limit m i(j) → 0.
Projecting the self-energy matrix as prescribed in Eq. ( 22) we obtain where the sub-index L(R) corresponds to the different chiralities and the "prefactors" ξ while one can notice that ξ . Thanks to this decomposition, we can obtain I Σ (I) ij in terms of some constant pre-factors and the following three master integrals 11 ) It will prove useful to change variables as where the variable ω H has limits ω ± H = (p ± q) 2 + M 2 H . Upon this change of variables, the latter integrals become where ) is an energy variable introduced in order to better identify the regions of support of the δ functions in Eq. ( 30), as detailed in Appendix A (and not to be confused with the LNV scale introduced in Section 2.2).

Gauge (Z) contribution
The interaction between massive neutrinos and the Z boson is given by 11 From now on to ease notation we write p ≡ |⃗ p|.
with g the SU (2) L coupling constant and c W ≡ cos θ W the cosine of the weak mixing angle.Just like in the Higgs case, the second equality holds when neutrinos are Majorana.Following the same procedure as in Section 3.2, we arrive at with ω Z ≡ (⃗ p − ⃗ q) 2 + M 2 Z and γ µ the Dirac gamma matrices.Using the relations from Eq. ( 22) we find where the "prefactors" P (I) L(R) are given by with the functions A 0 , A 1 and B given below The prefactors P L .Thus, we can define the following master integrals for the Z gauge boson case12

W -boson contribution
The interaction for the W -boson when including HNL is given by Even if the W -boson also contributes to the neutrino self-energy corrections, the structure of the Lagrangian is such that these terms can be readily taken into account from the Z-boson ones.Indeed, by making the changes g/(2c and noticing that in the W case we get P (I) R = 0 and P (2) L = 0 identically, one can obtain the neutrino self-energy corrections from W -boson exchange using the same master integrals as for the Z case.

Propagating states in the medium
Given the self-energy corrections to the propagator, written on full generality as14 we are interested in finding the propagating states in the medium, which follow Dirac's equation It proves useful to decompose the wave function Ψ in terms of helicity and chirality components, as [21] where the second term in Eq. ( 42) is the helicity operator, with the possible eigenvalues h = ±1.By using this decomposition of Ψ, together with Eq. ( 40), we arrive at the following set of coupled equations From the first line in Eq. ( 43) we can solve for ζ h as and substitute in the second line in order to find the equation of motion for φ h as We can identify the terms in brackets in Eq. ( 45) as the perturbed inverse propagator for left-handed fields in the medium, S −1 .Indeed, in the absence of self-energy corrections we recover the usual dispersion relation in vacuum: In the following, what we will need to do is to precisely find the dispersion relations for the neutrinos in the hot plasma, given by the complex zeroes of the inverse propagator from Eq. ( 45).

Simplified 2 neutrino case
It is instructive to consider an example where we only have the active neutrinos and the sterile one which will become the DM candidate.This allows to simplify expressions and obtain an analytical understanding of the different terms contributing to the production rates.In this "toy 2 × 2" scenario, the neutrino mass matrix, in the interaction basis, can be written as with tan 2θ ≡ 2m D /m DM ≪ 1, where θ is the active-heavy mixing angle.
The fact that we are interested in the production at temperatures around the EW scale allows to neglect the comparably lighter neutrino masses in the self-energy corrections, greatly simplifying Eq. ( 45).Additionally, we can work in the original flavour basis, in which the self-energy contributions become 15 where σ L(R) is an appropriate combination of the W and Z boson contributions to the retarded self-energies 16 .Thus, neglecting the contributions from Σ (2) , the inverse propagator, S −1 , becomes where we have defined the following combinations of self-energy corrections, that depend on the helicity eigenvalue (h = ±1): In order to find the dispersion relations in the medium for the propagating states, we need to find the complex zeroes of S −1 .These are given, to leading order in the mixing angle θ and neglecting light neutrino masses, by: From the imaginary part of p 0 we can directly obtain the relaxation rate of the DM, which is given by where ξ ≡ m 2 DM /(2p) and the rate depends on the helicity of the states, h, through the dependence of α h and Ω h ≡ 2p ∆ h + iγ h on it (see Eq. ( 50)).Although we are not interested in contributions to the DM production rate beyond those involving a SM boson, for completeness we show in Appendix C the resulting production rate within this formalism when including a direct coupling between the singlet neutrinos and a new scalar, which has been extensively studied in the literature [20,48,58,59].The expression from Eq. ( 52) agrees with that of Ref. [21] for Dirac neutrinos, which would translate into α h → 1, and in which the so-called "effective mixing angle", θ h ef f , is introduced as We show in Fig. 2 the DM production rates through weak boson decays for positive helicity (left panel) and negative helicity (right panel) for Dirac neutrinos (dashed lines) and Majorana neutrinos (solid lines).The active-DM mixing is set to |U α4 | ∼ 10 −6 and m DM = 10 keV.For RH helicity neutrinos, the term proportional to α h does not contribute considerably to the production rate as it is suppressed through the momenta of the neutrino, which is of the order of the temperature.Only very low momenta DM would have a non-negligible contribution, but the probability of producing such a neutrino is suppressed by the distribution function.Thus, most of the production of RH helicity DM goes through the second term in Eq. ( 52), which depends on Ω h .This second term is the only one present for Dirac neutrinos, and thus both rates are qualitatively similar in the left panel.On the other hand, LH helicity neutrinos present a completely different production rate depending on whether they are of Majorana or Dirac Nature.The outstanding difference arises because, when neutrinos are Dirac particles, only the second term in Eq. ( 52) is present.In the SM, weak interactions only couple to LH neutrinos, which receive large self-energy corrections, as can also be noted from the first line of Eq. ( 51), where Ω h appears at leading order playing the role of a "thermal mass".This in turn translates into a suppression of the effective mixing angle defined in Eq. ( 53) which itself suppresses the production rate (dashed orange line in the right panel of Fig. 2) 17 .On the contrary, if neutrinos are Majorana particles, then the term α h is present thanks to the Z boson couplings to RH neutrinos, and the production rate for LH helicity DM is greatly enhanced, and qualitatively similar to the one for RH helicity DM (although not equal because of the different W -boson contribution).

Decay rates for different contributions
We start our discussion by studying the DM production rate, Γ h s , as a function of the temperature through τ ≡ M W /T , and subsequently as a function of momenta through y ≡ p/T , variables which will prove useful later on to estimate the total DM abundance that can be generated through this Figure 3: Production rates in the "toy 2 × 2" scenario without HNL (dashed lines), the "ISS-like" scenario (dotted lines) and the type-I seesaw (solid lines) for RH helicity neutrinos (left panel with cold hues) and LH helicity neutrinos (right panel with warm hues).Additionally, dash-dotted lines represent the results from the type-I seesaw without the contribution from the Higgs.The vertical black dotted line represents the temperature at which SSB takes place and above which the results do not apply.The solid black line represents the Hubble expansion rate.The momenta is fixed such that y ≡ p/T = 1/10, and the mass of the HNLs, when present, is approximately given by m N .The active neutrino-DM mixing is set to |U α4 | ∼ 10 −6 and the DM mass to m DM ∼ 10 keV.
mechanism.Throughout all our discussion, we set the mass of the DM candidate to m DM = 10 keV, its mixing with active neutrinos to |U α4 | ∼ 10 −6 , and the mass of the heavy neutrinos, when present, to m N ∼ 150 GeV, unless otherwise specified.Additionally, whenever heavy neutrinos are present, we set their mixings with active neutrinos such as, at most, they saturate current bounds on them [60].
In Fig. 3 we show the production rates for both RH helicity (h = +1, left panel using cold hues) and LH helicity (h = −1, right panel on warm hues) DM as a function of the temperature, for the different cases we have motivated in Section 2.4 and summarized in Table 1.These were: • "Toy 2 × 2" scenario where only weak gauge boson decays to DM are relevant (dashed lines labelled "w/o HNL"), given that the Higgs boson couples through neutrino masses, which are negligible in comparison with the other scales of the problem lying around the EW scale.
• Type-I seesaw, but only taking into account weak gauge boson contributions (dash-dotted lines identified as "Type-I w/o H"), to assess the impact of the possible heavy neutrino decays into a gauge boson and DM.
• "ISS-like" scenario where heavy neutrinos are present (dotted lines called "ISS-4 × 4") and thus the Higgs contribution needs to be taken into account for consistency.We work in the single family approximation, which as discussed in Section 2.2, translates into a suppression of the couplings.
• Type-I including all contributions, notably the Higgs decay channel (solid lines labelled "Type-I").
Whenever we compare the "ISS-like" construction with the type-I, we set the heavy neutrino scale and mixings to similar values in order to be able to extract meaningful conclusions.The black solid line represents the Hubble expansion rate as a function of τ , which serves to compare and verify that, indeed, Γ h s (τ ) ≪ H(τ ) and thus considering freeze-in production of DM through decays is well justified.Finally, the gray shaded area for small values of τ (larger temperatures) represents the temperature above the EW crossover at which the Higgs vev becomes zero.For temperatures above T SSB ∼ 160 GeV, we are in the symmetric phase where the results of the present study cannot be applied.In this regime production would rather go through 2 → 2 scatterings, but its contribution is beyond the scope of this work.Moreover, we assume the transition from the symmetric to the broken phase to be instantaneous.Thus, we consider for the Higgs and gauge bosons their T = 0 masses, leaving the inclusion of the full temperature dependence of the Higgs vev and thermal masses in the propagators for future work.Indeed, the temperature dependence of the Higgs vev would translate into a suppression of the gauge bosons and Higgs masses for temperature regimes close to T SSB .In this case, the inclusion of the Higgs and gauge boson thermal masses in the high temperature limit (HTL) would be necessary in order to have a consistent description of the production rate.The thermal mass of the Higgs, given by becomes particularly relevant in this regime, being the dominant contribution to the scalar mass and regularizing IR divergences that would appear otherwise [25,26,61].We have nonetheless checked that, under the assumption of the instantaneous phase transition, the inclusion of the Higgs thermal mass does not qualitatively modify our results.Moreover, we could consider as well the inclusion of thermal masses for the HNL, which can be estimated as On the other hand, there are numerous experimental constraints on the size of the mixing between active and heavy neutrinos, Θ, which can be translated into bounds on the size of the Yukawa couplings as Using the result from Eq. ( 56) in the expression of the mass, we find an upper bound on its size When compared with M 2 i , for the temperatures and masses we have investigated, and given the bounds Θ † Θ exp ii ≲ 10 −1 − 10 −2 , we find that thermal masses for the HNL are usually smaller than the actual vacuum mass.
The first feature we can observe by comparing both panels is that, in all instances, the production rates for both helicities are qualitatively similar (although not quantitatively).This stems from the fact that we have Majorana neutrinos, and thus, both Z and Higgs contribute to both helicities, while interactions with the W boson only contribute to the LH helicity neutrinos.
Second, we note that in all cases, except when including the Higgs contribution in the type-I seesaw like construction, the rates are qualitatively very similar, peaking around τ = M W /T = 1 with values O 10 −24 GeV.It is thus clear that, without the Higgs contribution, the presence of the heavy neutral leptons does not alter significantly the dynamics in the early Universe, and thus the amount of DM generated through decays involving gauge bosons, to be quantified in the next section, is expected to be negligible regardless of the neutrino mass model details and spectrum beyond the presence of DM and its mixing with active neutrinos.
On the other hand, focusing now on the solid lines in which we consider a type-I scenario and include all relevant contributions, the production rate gets considerably enhanced compared to the others, while still satisfying Γ h s (τ ) ≪ H(τ ).The most remarkable difference here is with the "ISS-like" construction, which has a very similar spectrum, but differs in the structure of the Yukawa couplings between singlet neutrinos and the Higgs, thus producing a negligible contribution.The reason for this suppression is that, in order to satisfy the bounds on the DM-active neutrinos mixing angle, the full ISS(2,3) model requires a hierarchy of values in the mass terms coupling right-handed neutrinos and sterile singlets [19]; it is clear that such a structure cannot be achieved in a simplified model, where the only way to comply with experimental constraints is to suppress the overall Yukawa scale, thus also suppressing the HNL-Higgs couplings.We stress that the "ISS-like" construction analysed here is a simplified model, and we leave the study of the full (more complex) (2,3) ISS for a future work.In particular, we find that the only way to reproduce a similar spectrum as in the full type-I with our single family approximation for the ISS is in general through a suppression of either the active-DM coupling or the active-heavy mixing, if not both, and thus the rates do not differ notably from the case without HNL.On the contrary, with the full type-I structure, cancellations can lead to spectra satisfying all constraints while featuring large couplings.We conclude then that the single family approximation, applied in this case for the ISS, cannot completely describe all the possibilities full models describing oscillations offer, as argued in previous sections, and it is therefore necessary to work with full neutrino mass models.Note that for the full type-I seesaw scenario the rate tends to increase with temperature, to be expected from the possible Bose enhancement factor in the final state of the heavy neutrino decay n h → H + n DM .We emphasise however that these results are only applicable for temperatures below the EW crossover T SSB .
We show in Fig. 4 the production rates, using the same color code for the different helicities as before, for several values of the ratio of momenta over temperature, y, for the type-I construction.For large temperatures all the rates saturate to similar values, while larger values of y translate into a larger production rate at smaller temperatures.
When estimating the DM abundance produced through this mechanism in the next section, we will indeed need to take into account the rate for any y as well as for τ ≥ τ SSB .We show in Fig. 5 a comparison of the production rates as a function of y between the different cases, for τ ≡ M W /T = 1 fixed.It is again clear that the one family approximation for the ISS does not contribute further to the production than the "toy 2 × 2" case.The peak of the production appears at p ∼ T /20, even if in the full type-I scenario, the production rate is somewhat broader given the various decays contributing.
From the production rates in Fig. 5, and more clearly in Fig. 6, it is obvious that the probability to produce both a large or low momenta DM neutrino decreases exponentially, while around values p ∼ T /5 the maximum of the production is achieved.Another interesting feature is the reduction of the rate as the Universe expands (larger values of τ ).Indeed, the decays can only happen when the heavy mother particles (mainly the heavy neutrinos and Higgs) are present in the thermal bath.At temperatures much below their masses, their abundance starts to be Boltzmann suppressed, and thus the production rate becomes negligible when the temperature of the Universe is below T ∼ m N /15, the usual freeze-out temperature of the mother particles.
In the next section, we will estimate the final DM abundance in the two limiting cases found here, namely the "toy 2×2" scenario where only light neutrinos and DM are relevant, and with a benchmark point within the full type-I description which possibly renders a much larger DM abundance beyond the Dodelson-Widrow contribution through oscillations and collisions at lower temperatures [1].4) for the four benchmark points studied in Fig. 7, identified by their DM production neglecting thermal effects.In all cases the HNL masses are m N ∼ 150 GeV up to small corrections, while in the last column we additionally have non-zero Majorana phases in the PMNS mixing matrix, α 1 ∼ 1.06π and α 2 ∼ 0.38π.
production occurs when the energy budget of the Universe is dominated by its radiation content 19 , hence a(t) = T 0 /T (t) with T 0 being the present day temperature.By introducing the new variables, y = p(t)/T (t) and τ (t) = M W /T (t), the Boltzmann equation can be rewritten as: where we have conveniently defined Γ h prod ≡ Γ h s f DM,eq , while the Hubble parameter can be written as: in a radiation dominated Universe.The DM number density can be computed by integrating over both τ and y: The fraction of DM generated through this mechanism, normalized to the observed DM relic density, can thus be written as [62]: where g DM and g d ∼ 100 are the number of degrees of freedom for DM and the number of relativistic degrees of freedom at decoupling, respectively.As a first estimate for the amount of DM that can be produced through freeze-in decays of heavier particles, we compute the relic density in the "toy 2 × 2" scenario, as it allows to efficiently integrate over both p and T .On general grounds, even if the Majorana nature of the neutrinos translate into qualitatively different results to those found in Ref. [21], we find good agreement with the estimates done there, finding that the produced DM abundance from weak gauge boson decays is orders of magnitude below the observed one, given the observational constraints on the neutrino mixing.Particularly, the fraction of DM that can be accounted for through freeze-in is f DM ≡ Ω DM /Ω exp DM ∼ 10 −10 .On the other hand, as extensively discussed in Section 4.1, full neutrino mass models contribute to the production through decays of the HNL to DM and the Higgs.This can be further observed in Fig. 7 where we show the rates for different points in parameter space for the type-I seesaw with a HNL scale O(150) GeV.The solid lines correspond to m DM = 5 keV, while the rest have masses m DM = 10 keV, with different complex angles in the CI parametrisation, specified in Table 2. From here we note that, while complying with all bounds 20 we can find rates ranging from a fully thermalised species (solid   2, translating into slightly different rates.The DM mass is set to 5 keV for the solid lines while m DM = 10 keV for the rest.The solid black line represents the Hubble expansion rate.lines) which would freeze-out (Γ h s > H) to production fully in the freeze-in regime satisfying Γ h s ≪ H.The red band roughly represents the interphase in which freeze-in stops being a complete description of the problem and neglecting the build-up of the DM abundance is no longer justified.This is given by the conservative assumption that, for rates satisfying Γ h s ≤ H/10, freeze-in is a complete description.Note that these rates change with the momenta p, so that, depending on this, one might be in the freeze-in or freeze-out regime for different regions of the distribution function at different temperatures.Nonetheless taking all this into account is beyond the scope of this work, and thus we estimate the final DM abundance within freeze-in for points of parameter spaces that satisfy Γ h s < H, leaving a full scan of parameter space and complete treatment of the Boltzmann equations for future work.Upon integrating Eq. ( 59) we find in general at most a 2-order of magnitude reduction of the final abundance.For instance, for the CI parameters corresponding to the dotted lines in Fig. 7 giving a DM fraction of around 14% when neglecting thermal effects, we find that f DM ≳ 0.4%21 , while the dash-dotted line produces, when including the thermal effects here studied, f DM ≳ 1%.Given that there are regions of parameter space ranging from freeze-in to rates as large as to thermalise, we conclude that it is interesting to thoroughly study these cases, but leave the full and precise solution of the Boltzmann equations, as well as the scan of the parameter space for future work.

Conclusions
Despite all experimental and observational efforts to find and understand the nature of dark matter in the last decades, we still lack a compelling signal in any one of the current searches other than gravitational evidence.It is therefore interesting to consider DM candidates beyond the WIMP paradigm, which might not interact with our direct detection experiments.
In this context, it proves extremely appealing to consider a joint origin for the DM as well as for light neutrino masses.Indeed, under the general assumption that neutrino masses are generated through the introduction of singlet RH neutrinos, whose Majorana mass could lie at any scale as it is not set by EW symmetry, a mostly singlet neutrino, at or around the keV-scale, would represent a good DM candidate.Given its small couplings due to mixing suppression, the production of such a DM candidate would in general be non-thermal, helping also to produce a colder spectrum to comply with structure formation bounds.
The most minimal scenario, accounting just for the DM and its mixing with active neutrinos, produces DM through neutrino oscillations and collisions which build up a large enough DM population, while keeping it always out of equilibrium in the early Universe.This is the well-known Doddelson-Widrow mechanism, which however cannot account for all the DM given the bounds from X-rays.It is therefore necessary to account for DM through another production mechanism, for example that of freeze-in through particle decays.
While freeze-in production has been extensively studied before through the introduction, for example, of scalars directly coupled to the DM, here we have instead studied DM freeze-in through SM boson and HNL decays to DM.As soon as one introduces the DM candidate mixing with active neutrinos, weak gauge bosons are subject to decay to DM and one SM lepton.Moreover, the presence of heavier singlet neutrinos is necessary once neutrino oscillations are to be explained, automatically making Higgs and the HNL decays to DM themselves relevant for the DM production.Compared to the introduction of new scalars coupled in the dark sector, these channels represent an irreducible contribution to the final DM abundance on top of the DW mechanism, which becomes active at much colder temperatures, and thus it is necessary to understand whether enough DM can be produced or not through decays to DM involving a SM boson.For this purpose, there has been a number of works studying some of these channels separately.
On the one hand, Refs.[8,18] focused on DM production through HNL decays to the Higgs boson and the DM.The advantage of this channel is that it is not directly related to the heavily constrained active-DM mixing, thus potentially decoupling the DM production from the DM decay ones allowing to search for it.Contrary, the case for weak-boson decays to DM and either a charged lepton or a mostly active neutrino was studied in Ref. [22].While in both instances the observed DM abundance could be successfully generated, important thermal corrections were neglected [21].
While the weak boson decays to DM could not generate a large enough DM abundance due to the large thermal corrections that SM neutrinos receive at EW temperatures, it was necessary to study how the presence of the heavy neutrinos, and consequently the introduction of the Higgs boson mediated channel, might alter this picture.Just from naive arguments, the possible decay of the HNL to a boson and the DM can potentially take advantage of Bose enhancement in the final state, in contrast to the expected Fermi blocking in the final states from scalar decays into a pair of fermions.
In this work, we have studied, for the first time, the production of a keV-scale neutrino DM species taking into account the dominant thermal effects when including HNL as well, extending over the work of Ref. [21].While we find indeed that weak boson decays do not contribute further to the DW mechanism, given their very large thermal corrections suppressing active-DM mixing at large temperatures, the Higgs channel still represents an interesting avenue which needs further exploration.
Among the different scenarios we consider, increasing in complexity from a scenario with only DM and an active neutrino to describing neutrino oscillations and DM within the type-I seesaw, we find that simplified scenarios working in the single family approximation, where only one of each type of neutrino is introduced (namely an active neutrino, the DM one, and a HNL) do not produce a large DM abundance beyond weak boson decays, even in the presence of HNL which could decay to the Higgs boson.This is because in these cases, given the bounds on active-heavy neutrino mixing and on active-DM mixing from X-rays searches, couplings always need to be suppressed and thus production rates are negligible.On the contrary, working in full models completely accounting for oscillation data opens the door to a non-negligible production of DM through freeze-in.In these cases, one can still find solutions complying with all bounds while having not necessarily suppressed couplings.This is naturally realised if there is an underlying approximate lepton number symmetry, in which case light neutrino masses are also stable under radiative corrections.We explicitly verified this possibility in the type-I seesaw, finding solutions where the DM candidate even thermalises, although our focus was on regions of parameter space within the freeze-in regime.We studied several benchmark points, and found that the final DM abundance is, at most, suppressed by about three orders of magnitude with respect to the expectation when neglecting thermal effects.Given that there are regions of parameter space where an overabundance of DM can be found without thermal effects, we deem it promising to make a detailed study of this production mechanism with the thermal effects.
We leave a full scan of the parameter space for future work, given the expensive computational cost to estimate the final abundance for a single point, but motivated by the results we found here, where non-negligible rates for DM production are possible thanks to the presence of HNL and the Higgs boson in the plasma at EW temperatures.Among possible improvements that can also be taken into account, we highlight the introduction of thermal masses for the SM bosons (although we do not expect a qualitative change of the conclusions here drawn), as well as the Higgs vacuum expectation value evolution which would open the window to considering DM production beyond the SSB temperature.In this context, it is also necessary to include other production channels, namely 2 → 2 scatterings, as well as to consistently describe the LPM effect for temperatures above the EW crossover, in order to have a complete description of the production of the DM at all temperatures.As already stressed, figuring out the contribution to DM from these channels, taking into account the dominant thermal effects below the EW crossover, is of relevance given that whenever considering keV neutrino DM and the origin of neutrino masses in realistic models, both the DW contribution as well as this contribution from HNL decays to DM and the Higgs would contribute and generate the minimal DM within the model.Whether these two mechanisms working at different temperatures can explain the observed DM abundance or not would thus translate into the necessity or not of new physics in the dark sector on top of the singlet neutrinos when considering this species of DM.

Acknowledgement
We are extremely grateful to E. Fernandez-Martinez for his valuable comments and discussions which were at the origin of this work.This project has received funding /support from the European Union's Horizon 2020 research and innovation programme under the Marie Sk lodowska -Curie grant agreement No 860881-HIDDeN and under the Marie Sk lodowska-Curie Staff Exchange grant agreement No 101086085 -ASYMMETRY.ML is funded by the European Union under the Horizon Europe's Marie Sklodowska-Curie project 101068791 -NuBridge.

A Regions of support for massive particles running in the loop
In general, upong integration over the loop-momentum, q, we find four different combinations of Dirac delta functions, which define the regions of support of the integrals, over which their contribution is non-vanishing.For a boson of mass M B22 running in the loop we have in general which can be divided in the following different cases (taking only the T -dependent part): The second Dirac delta forces p 0 > 0 thus describing the decay n i ⇆ n k B and its inverse process.
Let us obtain the regions of support for q in detail in this case.We know that ω B ∈ ω − B , ω + B , which translates into By squaring on both sides we arrive at the following equality which can only be satisfied if The solution to the inequalities from Eq. ( 64) is then found to be q± = This is only valid if p 0 > q2 ± + m 2 k and Eq. ( 66) are satisfied, and defines the region of support for this case.One can check that in the limit of massless fermions running in the loop the result from Ref. [21] is recovered.
We need p 0 < 0 in order for the second Dirac delta to contribute, given that ω k (ω B ) are positive definite, and describe a similar process as before changing neutrinos for antineutrinos.The regions of support in this case are q± = 3. [f B (−p 0 + ω k ) + f F (ω k )] δ(q 0 − ω k )δ(p 0 − ω k + ω B ): For p 0 > 0 this term would describe the process n k ⇆ n i B, while for p 0 < 0 one would have B ⇆ n k n i .We have the following for the integration limits over q

B Useful integrals with distribution functions
Throughout the computations we need to perform some integrals over the distribution functions of bosons and fermions.In particular, we make extensive use of the following ones 23 : C Example of production rate with new scalars While we have extensively discussed SM and HNL contributions to the DM production, for completeness, we consider here the case in which the dark sector is comprised of other species, for example a singlet scalar, which could be at the source of the DM mass.This could be motivated by the assumption that lepton number is not explicitly broken in a seesaw scenario but rather dynamically when a singlet scalar develops a vev.On general grounds, we can have the following Lagrangian24 where ϕ is a singlet scalar with lepton number L ϕ = −2.The scalar potential, V , is such that, upon spontaneous symmetry breaking, both the SM Higgs boson Φ and the scalar singlet ϕ develop a vev, v H and v ϕ respectively.In the mass basis, we obtain a new interaction term between the mostly singlet scalar, φ, and the neutrino mass eigenstates One can directly compute the contribution of such an interaction to the dispersion relations for the neutrinos, just by exchanging v H → v ϕ and A L ij → m i − C ij m i − C * ij m j in the Higgs contribution computed in Section 3. To leading order in the mixings, and neglecting the contribution from light neutrino masses, the self-energy correction arising from the mostly-singlet scalar 25 , Ω h ϕ ≡ 2p(∆ h ϕ +iγ h ϕ ), contributes to the DM production rate as It is easy to notice that, in the absence of this extra scalar interaction, we recover the original result which depends on the neutrino mixing, see Eq. ( 52).More importantly, there is now a contribution to the DM production which is not suppressed by the neutrino mixing θ, and proportional to the Yukawa coupling in the dark sector, Y N , thus allowing for an efficient DM production as expected [48,58,59,89].

R
can be obtained by changing C ij → C * ij in the expressions of P (I)

Figure 2 :
Figure 2: DM production rates for RH helicity (left panel) and LH helicity (right panel) neutrinos in the "toy 2 × 2" scenario.Solid lines are for Majorana neutrinos while dashed ones for Dirac neutrinos.The momenta has been fixed such that y ≡ p/T = 1/10.The solid black line represents the Hubble expansion rate as a function of τ ≡ M W /T , while the gray shaded area represents temperatures above the EW crossover for which the results do not apply.For the neutrino parameters we fix m DM = 10 keV and |U α4 | ∼ 10 −6 .

Figure 4 :
Figure 4: Production rates in the type-I seesaw for different values of y ≡ p/T as a function of the temperature.The active neutrino-DM mixing is set to |U α4 | ∼ 10 −6 and the DM mass to m DM ∼ 10 keV, while the HNL mass is given by m N .Left panel with cold hues shows the production for positive helicity DM while the right panel with warm hues for negative helicity.In both, the solid black line represents the Hubble expansion rate.

Figure 7 :
Figure 7: Production rates for positive (left panel in cold hues) and negative helicity (right panel with warm hues) DM, for the type-I seesaw, as a function of τ for different points in parameter space.The HNL mass is always set around m N ∼ 150 GeV, while the CI parameters change as showin in Table2, translating into slightly different rates.The DM mass is set to 5 keV for the solid lines while m DM = 10 keV for the rest.The solid black line represents the Hubble expansion rate.

Table 2 :
DM mass and CI parameters from Eq. (