Freeze-In of radiative keV-scale neutrino dark matter from a new $\text{U}(1)_\text{B-L}$

We extend the Dirac Scotogenic model with the aim of realizing neutrino masses together with the mass of a keV-scale dark matter (DM) candidate via the same one-loop topology. Two of the Standard Model (SM) neutrinos become massive Dirac fermions while the third one remains massless. Our particle content is motivated by an anomaly free $\text{U}(1)_\text{B-L}$ gauge symmetry with exotic irrational charges and we need to enforce an additional $\mathcal{Z}_5$ symmetry. The dark matter candidate does not mix with the active neutrinos and does not have any decay modes to SM particles. DM is produced together with dark radiation in the form of right handed neutrinos via out of equilibrium annihilations of the SM fermions mediated by the heavy B-L gauge boson. In order to avoid DM over-production from Higgs decays and to comply with Lyman-$\alpha$ bounds we work in a low temperature reheating scenario with $4\;\text{MeV}\lesssim T_\text{RH}\lesssim 5\;\text{GeV}$. Our setup predicts a contribution to $\Delta N_\text{eff.}$ that decreases for larger DM masses and is below the sensitivity of upcoming precision measurements such as CMB-S4. A future observation of a signal with $\Delta N_\text{eff.}\gtrsim 0.012$ would exclude our scenario. We further sketch how inflation, reheating and Affleck-Dine baryogenesis can also be potentially realized in this unified framework.


Introduction
In most models of the Scotogenic variety one uses the lightest stable particle from the loop diagram for the neutrino masses as a DM candidate. For the original scotogenic model [1][2][3][4] this is either the lightest neutral component of the inert scalar doublet η or the lightest sterile neutrino produced as thermal WIMPS. In the Dirac version of the model [5,6] the DM can either the lightest neutral component of η or the singlet σ or the lightest vector-like neutrino. Later it was realized that keV scale FIMP DM is also possible in the scotogenic picture [7], but no mechanism was proposed for why this particular sterile neutrino is so much lighter than the other two. Reference [8] analyzed a model based on the DFSZ axion scenario [9,10], where a one loop diagram with vector-like fermions generates the keV-scale Majorana masses for a DM candidate. The authors of [11][12][13] showed that it is possible to construct models in which both the active and the sterile neutrino masses are obtained from loop diagrams. Recently a loop based extension of the seesaw scenario [14][15][16] with keV to GeV scale Majorana dark matter was put forth in [17], where two different scalar couplings were responsible for the mass generation and production from out of equilibrium Higgs decays. Unlike previous constructions we focus on Dirac neutrinos. We choose an abelian gauge symmetry as the guiding principle for building our model. After reviewing the Dirac Scotogenic model in section 2.1 we introduce our mechanism for generating the DM mass via a dimension five operator that resembles the Weinberg operator [18] in 2.2. In 3.2 we find that producing such a dark matter from out of equilibrium Higgs decays is not compatible with Lyman-α bounds on the DM mass. Section 3.4 demonstrates that the gauge symmetry is crucial for producing the correct amount of DM in the freeze-in scenario. We compute the minuscule amount of dark radiation produced by a similar freeze-in process in section 4. The necessary cosmic history can be realized in an inflationary context as explained in section 5. We close by illustrating how our set-up can potentially realize Affleck-Dine baryogenesis [19] in section 6.
2 The model

The Dirac Scotogenic model
Let us begin by reviewing the most salient features of the Scotogenic Model for Dirac neutrinos [5,6]. The goal is to generate the first diagram in figure 1. We follow the treatment of [20], where a U(1) symmetry is imposed on the fermionic sector that gets softly broken by the following trilinear term in the scalar potential where κ is a dimensionful parameter of mass dimension one. Here H is the SM Higgs and η, σ are inert doublet and singlet scalars with charges under the new symmetry. All particles and charges can be found in table 1. We start from U(1) B-L and assign conventional B-L field SU(2) L U(1) Y U(1) B-L Z 5 generations Table 1. Charges and representations for all particles participating in the neutrino or dark matter mass generation. The integers n in the fifth column are an abbreviation for ω n , where ω = e charges -1 to L and e R , whereas the right handed neutrinos ν R have the charge Q 1 = −1 so that the tree level mass term L H † ν R is forbidden by the symmetry. Here = iσ 2 denotes the anti-symmetric tensor in two dimensions. To generate this operator at loop level requires a soft U(1) B-L breaking by 1 + Q 1 units. Since we assume that H is uncharged under the new group, this means that the term κ η † Hσ has to have the same total charge Q η − Q σ = 1 + Q 1 . This soft breaking can be UV completed by considering the vev κ = λ IV v φ of another singlet scalar φ with charge −1−Q 1 , as will be shown in section 2.3. On the fermionic side we introduce two generations of vector-like pairs of SM singlets (N L , N R ) with B-L charge Q N In order to forbid a Dirac mass with L and ν R we have to require that Q N = ±1, ±Q 1 . We also need to forbid the following operators [20]: • H † η H † η with 2Q η together with N c R , N R would create ν c L ν L at loop level [1] • σσ with 2Q σ together with N c L , N L would create ν c R ν R at loop level [17] All of the above combinations of charges need to be non-zero and not divisible by |1+Q 1 |. If they were divisible by the only source of soft breaking, then an integer number of insertions of the trilinear scalar coupling in some loop diagram can generate the unwanted mass term. Once we know Q 1 we can fix all the other charges of the model. We will use the criterion of anomaly freedom to determine the rest of the particle spectrum and to find Q 1 in the next section. Before we do let us continue with our short review of the Dirac Scotogenic model: The active Dirac neutrino mass arises due to the first diagram in 1 and depends on the mass mixing in the scalar sector: After we expand all the fields into their components and in the absence of CP -violation there is no mass mixing between the CP -even (subscript R) and odd bosons (subscript I). We set m 2 η , m 2 σ > 0 in order to have an inert doublet and singlet. The real and imaginary components only mix among each other. The mass matrix after EWSB reads (2.8) and the same holds for the CP -odd fields, wherẽ m 2 η ≡ m 2 η + (λ Hη 1 + λ Hη 2 ) v 2 H , andm 2 σ ≡ m 2 σ + λ Hσ v 2 H . (2.9) We find two mass eigenstates in each case with the masses 10) and the mass eigenstates read The mixing angle is given in terms of the model parameters as Four diagrams contribute to the active neutrino masses: one mediated by each of the scalars S 1,2 and A 1,2 . Since S 1 and A 1 (S 2 and A 2 ) are mass degenerate there are only two distinct types of diagrams: two for heavier scalars of mass m 1 and two for the ones with m 2 . Due to the mixing there will be a relative sign between these two "generations" of scalars. This difference cancels out the divergent part leaving us with a finite mass matrix [21] ( N is the mass of the k-th heavy neutrino. To get a more insightful expression we work in the radiative seesaw limit [1] (2.14) After substituting in the mixing angle from (2.12) we find where the dependence on the soft symmetry breaking coupling κ is explicit and the scaling 1/M N is reminiscent of the familiar tree level Seesaw mechanism. To get a feeling for the involved scales let us estimate the neutrino mass in the single generation limit where in the above we used m 0 = O (1 TeV). Constraints on this scenario from lepton flavour violation and collider searches can be found in [22]. Note that since we will investigate a different implementation of Dark matter compared to the usual Scotogenic idea, we can push the masses of the scalars and N to values (far) above the electroweak scale, avoiding all laboratory constraints.

Extension for radiative DM mass
We proceed by introducing four Weyl fermions which are chiral under U(1) B-L . Usually one charges three right handed neutrinos with Q B-L = 1 so they form a vector-like pair with the ν L from the leptonic doublet ( the e L form a vector-like pair with e R ). However there are other anomaly free choices such as two right handed neutrinos with Q B-L = −4 accompanied by another one with Q B-L = 5. The idea of having chiral charges was originally put forth in [23] and applied to dark matter in [24][25][26][27][28]. Here we propose a new realization of this idea: Two Weyl fermions will be right handed and of equal charge Q 1 in order to form two massive Dirac fermions with ν L . Therefore our model predicts that the third SM neutrino remains exactly massless. The remaining two fermions will be the right handed χ 3 and the left handed χ 4 , which combine to form a Dirac fermion, that will be identified with the dark matter candidate. Since we have a gauge symmetry in mind, we need to find an anomaly free set of charges. As we only consider SM singlets with chiral charges there are only two conditions for cancelling the gravitational and U(1) 3 B-L anomalies from the Standard Model: Here the signs reflect the fact that only χ 4 is left handed. The system of equations is underdetermined and has infinitely many solutions. In order for the same one-loop topology and soft breaking to generate the dark matter mass term Without the absolute value we find no solutions. For 1 + Q 1 = −(Q 3 − Q 4 ) we find two possible solutions with irrational charges and One can see that both sets of solutions are related by exchanging Q 3 ↔ −Q 4 . The only solution possible for 3 copies of ν R with the same charge Q 1 would be Q 1 = −1 and Q 3 = Q 4 , which would allow for a term L H † ν R at tree level and hence will not be investigated further. This is why our model predicts only two massive SM neutrinos. Note that formal quantum-gravitational conjectures [29] seem to exclude abelian gauge theories with irrational charges in curved space-time. We do not consider this line of reasoning further for our purely phenomenological study. Let us emphasize that for this particle content we need a soft breaking by |1 + Q 1 | = 1 unit. However in that case any of the previously mentioned unwanted mass terms could arise at the loop level via some number of insertions of the trilinear term. Furthermore since we break the gauge symmetry by only one unit, there will be no residual Z N symmetry that also stabilizes the dark matter. To remedy both shortcoming we resort to imposing an ad-hoc Z 5 symmetry as well. The choice of an odd N was motivated by the need to forbid bilinear terms. All the charges and representations to realize the original Dirac Scotogenic model [5,6] with our exotic choice of U(1) B-L charges can be found in the table 1. Let us focus on the dark matter mass now: Motivated by Zee's model [30] for neutrino masses we consider the second topology depicted in figure 1. We add a pair of vectorlike doublets (D L , D R ) with Y = −1/2 together with another pair of vector-like singlets (S L , S R ) with Y = 0: Here we coupled the fermion χ to η, σ * instead of η † , σ, which was the case for the active neutrinos, because we need a soft breaking by plus one unit of B-L, whereas the active neutrinos needed a breaking by −1. Since both components of χ are SU(2) L singlets unlike for the SM leptons, we do not only need a chirality flip on the internal fermion line but an insertion of the Higgs doublet as well: These couplings are the reason why all D, S have the common B-L charge 1 + Q 3 see table  1. B-L forbids all Majorana masses S c L S L , S c R S R , χ c L χ L , χ c R χ R as they need a breaking by 2(1+Q 3 ), 2Q 3 , 2Q 4 units. Since Q 3,4 are irrational numbers, no loop graph with an arbitrary number of soft symmetry breaking insertions by one unit can ever accidentally produce these terms. Hence we will leave χ L,R uncharged under the Z 5 . This automatically forbids any mass mixing between χ and the ν L,R as well as the kinematically allowed radiative decay χ → νγ or the three-body decay χ → ννν. Consequently the DM candidate will be absolutely stable. All other mass terms of the schematic form LD, DH † e R , SN , Sν R , Sχ, N χ are each forbidden by at least one of the symmetries or both. The dark matter mass term from figure 1 depends on the mass mixing in the scalar sector as well as on the mixing between the D and S. Their mass matrix reads and we find the following eigenvalues The diagonalization simplifies in the limit Y SD = −Y DS and we arrive at with a mixing angle The dark matter mass arises due to eight loop diagrams in the mass basis. Since S i and A i are mass degenerate there will be only four distinct kinds of diagrams. For a fixed intermediate F j there are two diagrams depending on S 1 (A 1 ) and S 2 (A 2 ) again with a relative sign. Consequently all divergences will cancel in the sum and the resulting DM mass is finite. For a fixed intermediate S i (A i ) there are two possible diagrams involving F 1 and F 2 , both with a relative sign due to the fermionic mass mixing. This explains the structure of the expression for the DM mass: with C 2 = −C 1 = 1. By working in the radiative seesaw limit and invoking the definition of the mixing angles (2.12) and (2.27) we finally obtain Note that since we generate the dark matter mass via a dimension five operator compared to the active neutrinos (see figure 1), whose mass is an effective dimension four operator, there is another inverse power of the heavy suppression scale M F when compared to (2.15). Because we want our dark matter to be heavier than the neutrinos we therefore need M N M F , which can be seen from the following estimate In the above we used m 0 = O (1 TeV). Unlike the N which are much heavier the F and electrically charged components of D could be potentially be produced at future colliders and have a direct coupling to the SM like Higgs.

UV completion
In order to gauge the U(1) B-L symmetry and to explain the origin of the dimensionful coupling κ in the trilinear term (2.1) we introduce a second SM singlet scalar φ with the charge Q φ = −1 − Q 1 = 1 without any couplings to the fermion spectrum: We parameterize the new scalar as which allows us to identify κ = λ IV v φ ≡ λ IV v B-L . We do not depict an insertion of this vev in figure 1, because the neutrino and DM mass generation only requires a non-zero value of κ irrespective of its origin in the UV. The φ 0 I is the would-be-Goldstone-Boson that gets absorbed to become the longitudinal component of the massive U(1) B-L gauge boson that we call Z whose mass reads because φ is the only field with B-L charge that receives a vev. Direct searches at LEP place the following bound [31,32] on the mass of a new gauge boson v B-L = m Z g B-L > 6.9 TeV @ 95% C.L. (2.36) that couples to the conventional B-L charges of the SM fermions. Searches at the LHC exclude Z s below 0.2 − 3.5 TeV [33]. Since no scalar field that receives a vev is charged under both weak isospin/hypercharge or B-L there is no mass mixing between the Z and Z' bosons. However there can be gauge kinetic mixing [34], for instance generated at the loop level by self-energy graphs containing the (D L , D R ) or η fields, which are charged under both abelian symmetries and weak isospin. The additional scalar interactions contribute to the masses of the η 0 and σ 0 bosons by shifting the relations in (2.9) tõ Additionally the mixed quartic between H and φ leads to mass mixing between them: First we minimize the potential in each direction and find expressions to eliminate the parameters µ 2 H , µ 2 φ < 0. We find that the minimum in each direction can be obtained for and we arrive at with the eigenvalues In the limit v B-L v H we find at leading order The correction to the SM like Higgs mass can be understood as a tree level threshold correction to its quartic from integrating out the heavier field [35]. The mass eigenstates are determined from and at leading order in v H /v B-L this reduces to In the present study we will neglect this mixing completely. It is important to note that the discrete Z 5 symmetry we imposed will most likely be broken by quantum gravitational effects [36][37][38], which is why we assume it is e.g. a residual symmetry arising from the spontaneous symmetry breaking of a gauge symmetry [39]. This larger symmetry could also connect our choice of U(1) B-L with the rest of the SM gauge group, e.g. by unifying it with QCD into the Pati-Salam hypercolor SU(4) c [40]. Vector-like fermions such as our singlets (N L , N R ) and (S L , S R ), doublets (D L , D R ) as well as exotic vector-like downtype quarks arise in E 6 -GUTs [41,42]. This could provide an interesting route for further completing our model in the UV as the Pati-Salam model can be embedded in SO (10) which is a subgroup of E 6 .

Dark Matter
As previously mentioned our DM candidate does not mix with the active neutrinos. Hence the usually considered possibility of creating keV-scale neutrino DM via active-to-sterile oscillations [43], that can be enhanced in the presence of a chemical potential for neutrinos [44], are ruled out and we have to look into other avenues to produce DM. In the following we will briefly explain why we do not consider thermal production and focus on nonthermal scenarios. To study non-thermal production of DM we assume that the reheating temperature T RH of the universe is below both the masses of the particles in the loops of figure 1 and the mass of the B-L gauge boson Z This ensures that none of the new, potentially stable neutral particles, which are good thermal dark matter candidates, are present in the plasma. We can thus limit ourselves to the SM degrees of freedom augmented by the two ν R and the light DM.

Lyman bound for FIMPs
The Lyman-α forest consists of absorption lines in the spectra of quasars due to neutral hydrogen in the intergalactic medium. It provides a window into the matter power spectrum, which contains information on the Dark matter's free-streaming scale from the time of structure formation. One can use the existing data on the Lyman-α forest to set bounds on dark matter models affecting small scale structures such as the thermally produced warm DM (WDM). Numerically challenging simulations for WDM have been performed and lead to a lower limit of m Ly-α WDM = 5.3 keV at 95% confidence level (CL) [45,46]. Reference [47] argued that the aforementioned bound is too strong when taking into account systematics such as assumptions about the thermal history and instead they find m Ly-α WDM = 1.9 keV at 95% CL. In order to avoid such time consuming simulations for other DM production modes a bound mapping formalism has been devised in [48][49][50][51][52][53][54] and a recent reevaluation [55] found that the previously mentioned mass range m WDM (1.9 − 5.3) keV translates into a bound on the FIMP mass of m FIMP (4 − 16) keV.

Out of equilibrium Higgs decays
In the following we focus on the decay h → χ L χ R , χ R χ L , which is obtained from the second diagram in 1 by replacing one of the Higgs vevs with the radial excitation h, which leads Leading order diagrams for the decay h → χ L χ R in the mass basis. See the main text for more details.
to the two diagrams depicted in figure 2. By replacing both Higgs vevs one can compute the scattering process hh → χ L χ R , χ R χ L , however for our first estimate we will limit ourselves to the decays. We only consider the trilinear coupling from (2.1) and neglect all the decay modes to the same chiralities (LL or RR) which occur via the quartic couplings λ Hη1,2 , λ Hσ , from DM mass insertions on the external lines [17], or from mass mixing in the heavy scalar or fermion sector to focus on the parameters for the DM mass. This is also why we will work to the lowest order in the mixing angles sin(α) and sin(β) because in the mass basis there are 32 diagrams contributing and both neutrino masses were independent of the aforementioned angles in the radiative seesaw limit. We neglect the mixing sin(γ) between h and ϕ. In this approximation with cos(α) cos(β) 1 and S 1 η 0 R , S 2 σ 0 R , F 1 S, F 2 D 0 there are only four diagrams contributing. By dropping the final state DM mass we find The first set of graphs depicted on the left side of figure 2 is obtained by replacing the upper vev in 1 with h and only depends on sin(β). The amplitude is finite because it comes from a difference of terms due to the relative sign between the F 1,2 contributions. The corresponding effective Yukawa coupling is found to be at leading order in sin(β) from (2.12) and by making use of (2.29) Similarly the second diagram on the right side of figure 2 is obtained by replacing the lower vev in 1 with h. It is proportional to sin(α) from (2.27) and finite because here the difference arises due to the relative sign of the S 1,2 (A 1,2 ) contributions. The effective coupling is In both expressions we neglected the Higgs mass. The sum of both contributions can be re-expressed by comparison with (2.30) as which agrees with the EFT expectation that after EWSB the diagram on the right in figure  1 can be represented by an effective Weinberg-type operator [18] at energy scales below all the mediator masses The remainder of this section discusses how to produce DM from this effective operator and can be applied to other models that generate this operator as well. For the decay width we find after neglecting the phase space suppression and we emphasize that the only free parameter is the DM mass. The experimental limit on the branching ratio (BR) from searches for invisible Higgs decays beyond the SM is between 19% (CMS) and 26% (ATLAS) which translates to approximately Γ (h → Inv.) < 1. This bound is only one order of magnitude stronger than (3.8) due to quadratic dependence of the branching ratio on the DM mass. The invisible Higgs decays lead to the strongest terrestrial bound on the DM mass, however as we will see avoiding cosmological over-production of DM from Higgs mediated scatterings firmly requires the DM mass to be below the MeV-scale see (3.44).
In the following we will limit ourselves to the era of radiation domination and make extensive use of the Hubble rate and the entropy density where g * ρ and g * S are the effective number of degrees of freedom in energy and entropy respectively. Before we deal with non-thermal DM production let us take a look the thermal case first: The decay (3.7) will be in thermal equilibrium at T = m h provided that m DM 4.5 keV (we will show this later in (3.17)). Since during radiation domination we have Γ/H ∼ T −2 for decays at temperatures below the mass of the decaying particle, a decay never falls out of thermal equilibrium. Consequently we need to know when the inverse decay freezes-out in order to find the decoupling temperature of χ. The corresponding rate reads at T m h [62, 63] and the phase suppression is encoded in the Boltzmann factor. Numerically we find that this interaction freezes out at T FO 3 GeV for m DM 2 GeV. Of course there is also a scattering process hh → χχ, but since this requires two on shell Higgses the rate density will be double Boltzmann-suppressed below m h typically leading to an earlier freeze-out than the inverse decays. Since χ is relativistic at decoupling it would be a warm DM candidate, however it has long been known, that such a DM candidate would overclose the universe [64] if there is no release of entropy that dilutes the relic density to the observed value. Realizing the warm DM scenario requires additional degress of freedom in the plasma like long-lived particles that decoupled while relativistic whose decays generate the necessary entropy dilution [65]. For the sake of minimality we do not consider this idea further and focus on out-of-equilibrium-processes involving only SM states that are connected to the DM via the previously introduced BSM Yukawa and gauge interactions.
Next we investigate out of equilibrium Higgs decays. We use the notation of [66] to write down the Boltzmann equations for the DM production where Y DM ≡ n DM s , with s being the entropy density and z where we assumed that entropy is conserved. Note that away from thermal equilibrium the temperatures of the SM and DM baths are different so that γ h→χχ depends on T SM and γ χχ→h is a function of T DM . The freeze-in regime [67] is defined by the condition Y χ Y e.q. χ and the same for χ. If we use the fact that the SM Higgs is kept in thermal equilibrium Y h Y e.q.
where the thermally averaged decay width density reads [62] In this context g h = 1 is the spin degeneracy of the Higgs and K 1 (z) denotes a modified Bessel function of the first kind. To ensure that we are in the freeze-in regime the decay is not allowed to thermalize which leads to the condition that can be re-expressed as a bound on the DM mass m DM 4.5 keV · g * ρ (m h ) 100 1 4 , that is borderline compatible with the lower limit of the Lyman-α window. Under the assumption that there is no primordial abundance of DM we can integrate (3.15) to determine the DM abundance today at z 0 and the factor C h is a short hand for all microscopic and cosmological parameters We then use this to compute the energy density in dark matter by using the present day entropy density s 0 and the critical density ρ c [68] Here the factor of two arises because our DM candidate is a Dirac fermion. For a simple analytical estimate we can neglect the temperature dependence of the relativistic number of degrees of freedom in energy g * ρ (z) and entropy g * S (z) and replace them by their average values at the time of predominant dark matter production. This can be done because freeze-in production of DM is always sharply peaked around either T m h for the IR freeze-in [67,69] or at the reheating temperature T RH for UV freeze-in [70]. First let us suppose a standard big bang cosmology that corresponds to z RH → 0 which gives the maximally possible abundance that corresponds to where we used the maximum possible number of relativistic degrees of freedom above the EW phase transition in the SM. One can see that the correct relic density [71] is obtained for a DM mass that is in conflict with the more conservative Lyman-α bound that requires m DM > 4 keV. Since h 2 Ω DM ∼ m DM Y χ (z 0 ) we can allow for a larger DM mass by lowering the yield Y χ (z 0 ). This is most easily done by assuming z RH > 0 which lowers the relic abundance below (3.21). In doing so we introduce a second free parameter in the form of T RH . We find that we can decrease the abundance for z RH > 1, however our fix comes with two complications: On the one hand one needs to make sure that the SM Higgs is actually thermalized after reheating. Reference [72] found that particles charged under non-abelian gauge symmetries that are produced from inflaton decays during reheating thermalize before the end of reheating (which is not an instantaneous process) provided that the fine structure constant of the gauge interaction satisfies .

(3.23)
In this context m I is the inflaton mass and Γ I is its decay width, which we can trade for an expression involving T RH (see (5.2)). We find that which is definitely satisfied for the SM Higgs coupling to SU(2) gauge bosons where α 2 = g 2 4π with g 0.64. On the other hand the out of equilibrium condition (3.16) must be reevaluated at T RH < m h leading to The necessary z RH > 1 leads to DM masses which even violate the lower more conservative Lyman-α bound. In other words: If we tried to satisfy the Lyman-α window we would obtain a thermalized population of χ, which would actually be warm dark matter and this can only be made to work with additional processes that release enough entropy to dilute it. Since this channel leads to over-production of dark matter and the inclusion of 2 → 2 scattering processes will only increase the relic abundance further, we conclude that freeze-in from the SM Higgs via a Weinberg-type operator is not a viable production mode for keV-scale DM. Furthermore in order to avoid any contribution from Higgs decays we will only consider cosmologies with T RH < T FO

Super WIMP contribution
Another production channel for DM is the Super WIMP scenario [75] in which the DM is produced after the thermal freeze-out of the Higgs boson from its gauge and Yukawa interactions at T FO m h 25 . However the Higgs has decay modes to SM particles which are much faster than the decay to DM so the frozen out abundance of Higgses can not lead to a significant production of DM.

Gauge Scattering
We can also produce DM via the new gauge interaction [76,77]. In the limit s m Z the cross section for interconverting DM and SM fermions f i via Z exchange reads for massless fermions [78] σ where Q denotes the various B-L charges and N c is a color factor which equals three for quarks and one for leptons. The above was summed and not averaged over the initial state spins. Since m Z = g B-L v B-L the cross section is only sensitive to v B-L in the effective operator limit. Even though the DM mass in (2.31) formally depends on κ = λ IV v B-L we treat m DM and v B-L as independent parameters, because a larger v B-L can always be compensated by a smaller λ IV or by making the fermions running in the loop heavier.
where we applied the freeze-in approximation in the last step and for simplicity we compute the scattering rate densities via Maxwell Boltzmann-averaging [79,80] γ instead of averaging with Fermi-Dirac statistics. By neglecting the masses of the DM and SM fermions the simpler Maxwell-Boltzmann average allows us to find an analytical expression by employing the relation [80] ∞ Note that while the functional forms above are the same the densities depend on the different temperatures of the SM and DM baths. The fact that both densities are equal for equal temperatures reflects the principle of detailed balance, so that the right hand side of the Boltzmann equation vanishes in thermal equilibrium [81]. Owing to our previous simplifying assumptions we will only work with relativistic fermions in the SM plasma. Annihilations from non-relativistic fermions will be Boltzmann-suppressed at T < m f i and therefore less important than relativistic processes. From this we can deduce the more familiar interaction rate for relativistic SM fermions (g f i = 2) which agrees with the estimate based on dimensional analysis for an effective four fermion operator that leads to Γ ∼ T 5 /v 4 B-L . Our result is larger by only around 11% compared to the result [32] found by averaging over Fermi-Dirac statistics and also using massless fermions. This numerical difference agrees with the findings of [82] but we do not take percent level effects into account since what matters for freeze-in is the order of magnitude of the couplings and not their precise value. In the effective operator limit the scattering rate is UV dominated so its maximum value is found at the largest available bath temperature after completion of reheating given by T RH . As a consequence of our analysis for Higgs decays in 3.2 we will assume a reheating temperature 4 MeV T RH 5 GeV. Since the SM fermions also couple to non-abelian gauge interactions the estimate (3.24) is still approximately valid even though the SM fermions are not necessarily produced from inflaton decays. If we assume the inflaton decays to the SM like h, which definitely will be thermalized according to (3.24), and that h decays or scatters to produce the SM fermions it is plausible to expect a thermalized SM fermion bath immediately after reheating. Then in order to guarantee that we stay in the freeze-in regime the rate needs to satisfy and we can use this to constrain the B-L breaking vev to be v B-L 56.8 TeV · T RH 1 GeV which is a stronger constraint than the laboratory bound (2.35). For the above we summed over all the relativistic fermions because of the sum on the right hand side of (3.27). Moreover we used that for all SM leptons Q 2 l = 1, for quarks Q 2 q = 1 9 with N c = 3 colors and assumed all leptons and quarks except the top and bottom quark to be relativistic at T RH = 1 GeV. We compute the effective coupling of the relativistic SM fermions as Here we treat a fermions as relativistic as long as E 3T > m f . In the above definition the first 3 stands for the SM neutrinos and the contribution from the charged leptons and quarks is multiplied by a 2 because both chiralities produce DM. We only wrote out the contributions from the heavy quarks explicitly and the term 2θ (T − 200 MeV) is the contribution from u, d, s, whose mass is below the temperature of the QCD phase transition at T QCD 200 MeV. Below this transition all quarks hadronize and at least for a short period of temperature the light mesons are still relativistic and should be taken into account [78]. The inclusion of these particles requires the use of form-factors and we ignore them because they quickly become non-relativistic and hence the rate density becomes double-Boltzmann suppressed compared to the contributions from ν l and e − . Note that we can reuse this estimate to make sure that the same interaction does not equilibrate the ν R with charge Q 1 = −2: The cross section (3.26) also applies to ν R by replacing which is valid for both possible DM charge assignments (2.20) and (2.21). Therefore the ν R production rate is always smaller than the DM production rate so that (3.35) ensures that there is no thermal population of ν R . We proceed by analytically solving (3.28) (3.38) Here the reheating temperature T RH acts as a UV-regulator for the effective cross section and if we were to consider T RH → ∞ we would need to use the full kinematic dependence of the Z propagator to unitarize the rate. For the estimate we again replace the relativistic number of d.o.f with their values at T RH (see 3.2) so that (3.20) evaluates to .
Note that one can not set m DM to arbitrarily large values since we neglected the phase space suppression for the finite DM mass in (3.26). As a rule of thumb our results apply as long as m DM T RH . For the numerical evaluation we use the full temperature dependence of g * S and g * ρ by employing the fitting functions from [83], which agree up to less than one percent with the exact expressions except during the QCD phase transition and during e + e − annihilations, where the differences are about four percent. Figures 3 and 4 illustrate the behaviour of the DM abundance today for different values of the reheating temperature, DM mass and v B-L together with the observed relic abundance [71]. As previously alluded to one can see that the yield reaches its asymptotic value shortly after reheating and its final value strongly depends on T RH as expected for UV freeze-in. Before closing we would like to emphasize that the SM like Higgs can also mediate SM fermions annihilating to DM via the effective interaction in (3.6). The corresponding scattering rate density is found from (3.32) encodes the couplings of the relativistic SM fermions to the Higgs in analogy with (3.36). We neglect the coupling to the active neutrinos as it scales with their mass. The Higgs mediated interaction does not thermalize at reheating provided that 42) which is stronger than the bound from invisible Higgs decays (3.8). The estimate for the Higgs mediated relic abundance is straightforward and by comparing with (3.38) we arrive at . (3.43) If we demand that this additional contribution is smaller than the Z mediated one we obtain an upper limit on the DM mass of which was evaluated at T RH = 1 GeV, where all charged fermions except the top and bottom quark contribute. This represents the strongest upper bound on the DM mass and is the reason why we only consider DM with typical masses at the keV-scale. We depict contours in the T RH versus v B-L plane that reproduce the measured DM abundance today for multiple representative masses that agree with the Lyman-α bound and (3.44) in figure 7.

Dark matter phenomenology
Owing to our choice of Z 5 symmetry the DM is absolutely stable and does not mix with the SM neutrinos. Consequently the DM has no radiative decay mode to a ν L plus a photon. This decay constitutes the canonical signature of keV scale sterile neutrino DM that is being looked for via X-ray searches investigating the diffuse X-Ray background or dwarf galaxies [84][85][86][87] (see also [88] and references therein). The coupling to the Z and the Higgs induce velocity independent dark matter self interaction cross sections. However due to the small DM mass and v B-L m h m DM the resulting transfer cross sections [89][90][91] are far to small to help with the "cusp-core" and "too-big-to-fail"-problems [92][93][94][95] or even to come into conflict with bounds from the Milky way or the Bullet cluster [92,93,95]. The aforementioned self interactions could lead to efficient DM self scatterings which would lead to kinetic equilibrium of the DM population in the early universe. Because of the separation of scales between m h and v B-L we only focus on the individual contributions and ignore the interference term. For the Higgs mediated interaction we find in the limit and use the methods of section 3.4 to compute By comparing the interaction rate Γ h = γ h /n χ , where n χ is the DM number density with g χ = 4 degrees of freedom, to the Hubble rate evaluated at reheating we find that the DM  is not in kinetic equilibrium with itself at T RH as long as m DM 475 MeV · 5 GeV T RH 3 4 · g * ρ (T RH ) 85 1 8 .

(3.47)
Similarly to find Γ Z (χχ → χχ) we can reuse the result (3.33) by replacing the charges (see (3.36)) The Z mediated diagram does not equilibrate the DM with itself at reheating provided that v B-L 175 TeV · T RH 5 GeV 3 4 · 85 g * ρ (T RH ) 1 8 . (3.49) We conclude that scattering can lead to kinetic equilibrium of the DM at early times for certain choices of parameters. Since both rates arise from effective operators they decrease with temperature, which means that even if the DM was thermalized with itself initially it will fall out of kinetic equilibrium during the evolution of the universe. As a consequence of the constraint (3.8) we only investigate very light DM with typical masses below 2 GeV. Since nuclear recoil experiments basically have no sensitivity for sub-GeV DM due to kinematics, there has been a growing interest in studying atomic bound state electrons as targets for direct detection of light DM [96]. In order to estimate whether these targets can be used to find our DM candidate, we compute the Higgs and Z mediated cross sections for non-relativistic DM in the electron rest frame and expand to leading order in v DM 1: The Higgs mediated cross section comes with two more powers of both m e and m DM compared to the Z mediated one, because the couplings to the Higgs are proportional to the aforementioned masses. In the above we chose v B-L to reproduce the observed DM relic density for a given DM mass. The best current limit including form factors for bound state electrons is σ 10 −40 cm 2 [97,98]. One can see that direct detection via electrons is not a viable search strategy for our DM candidate owing to the small values of m DM and the large v B-L necessary for freeze-in.

Dark Radiation
The SM prediction for the number of relativistic neutrinos is [99][100][101][102][103][104] N eff. = 3.0440 ± 0.0002, (4.1) and the small deviation from the value of 3 expected for three generations of ν L comes from the fact that their decoupling from the SM bath is not instantaneous. Additional relativistic degrees of freedom are usually referred to as dark radiation (DR 0.28 @ 2σ C.L. .

(4.2)
Currently there is a lot experimental effort to improve this bound: The South Pole Telescope [105] and the Simons observatory [106] both aim to reach ∆N eff. 0.12 @ 2σ C.L. while the upcoming CMB Stage 4 (CMB-S4), experiment [107][108][109] and NASA's PICO proposal [110] have a sensitivity forecast of ∆N eff. = 0.06 @ 2σ C.L. There is also the planned CORE experiment by the ESA [111] with similar goals.

Dark Matter as dark radiation
Since the dark matter is out of equilibrium with the SM bath, its typical momentum after production can in principle be vastly different from the temperature of the SM. Even though the DM is non-relativistic today, it might have been relativistic at the time of BBN or CMB decoupling. One can find a condition for having a relativistic DM particle at the SM bath temperature T [55] Here T γ (t 0 ) is the photon temperature today and g * S (Tγ (t 0 )) g * S (T RH ) is the ratio in the number of relativistic degrees of freedom in entropy today versus the number at the time of DM production, which we approximate with T RH . a(T ) denotes the scale factor, whose value ranges from 10 −10 at the time of BBN (T 1 MeV) to 10 3 at the time of CMB decoupling (T 1 eV). Consequently our DM candidate can only be relativistic around BBN, but not at recombination. The contribution of the DM to ∆N eff. at BBN temperatures was found to be [55,112] ∆N eff. (T BBN ) 3.4 × 10 −4 · Ω DM h 2 0.12 · 10 keV m DM · 100 g * S (T RH ) 1 3

(4.4)
and is negligible compared to the expected sensitivities. Note that the above estimate relied on the FIMP being produced from a decay, however we do not expect production from scattering to significantly alter the order of magnitude of the result.

Right handed neutrinos as dark radiation
Due to their feeble effective Yukawa interaction with the left handed neutrinos (see (2.16)) the ν R never equilibrate with the SM [113] and the freeze-in of the aforementioned interaction contributes an even more negligible amount of [114] ∆N eff. 7.5 × 10 −12 · m ν 0.1 eV 2 (4.5) in standard Big Bang cosmology. Gauge annihilations of SM fermions via the Z can also create ν R . From section 3.4 we already know that if we want to produce the DM from freeze-in the ν R production will occur in the freeze-in regime as well. The corresponding cross section is given by (3.26) under the replacement (3.37) and α χi → α ν R i . We can write down the coupled Boltzmann equations for the evolution of the SM and DM energy densities [115] dρ SM dt + 3H (ρ SM + P SM ) = −C ρ , (4.6) where P denotes the pressure density. Adding both Boltzmann equations gives the result expected from the continuum equation Making use of the equation of state for radiation allows us to write The right hand side of the Boltzmann equations is known as the collision term and parameterizes the energy exchange between the SM and DM baths. It can be written as [82,116] where we neglect the back-reaction from the ν R bath in the freeze-in approximation in the second line. The quantities Eσ | v| are functions of the respective bath temperatures and are defined completely analogous to σ | v| in (3.29) as [79,82,116] δ (a + b → i + j + . . . ) = Eσ | v| n eq. a n eq. b (4.12) K 2 is the modified Bessel function of the second kind, which arises compared to the K 1 in σ | v| due to the presence of a factor of E in the thermal average. Again we use Maxwell-Boltzmann statistics instead of the correct Fermi-Dirac averaging to obtain simpler analytic results. By employing the relation we can compute the average for massless initial and final states (4.14) Note again that in the above one has to take into account that the rate densities depend on the different bath temperatures. The scaling of this energy exchange rate density is consistent with dimensional analysis as it scales like the rate density (3.32) for the DM abundance multiplied by another factor of T . Since we do not know the phasespace distribution function of the non-thermal ν R we do not know their temperature so we compute their energy density directly from solving the Boltzmann equation. If we neglect the energy loss of the SM bath, which is the basis of the freeze-in scenario and assume that the SM entropy is conserved we find [114] ρ ν R (T ) 2 · s SM (T ) in terms of the SM temperature T and use this to compute [114] where the first factor of two in (4.15) accounts for the fact that the ν R have g ν R = 2 spin polarizations and the second one in (4.16) for two generations of ν R . At temperatures below the electron mass e + e − annihilations heat the SM plasma compared to the decoupled species so that by using g * S (T < m e ) = 43 11 we recover the more familiar formula For the regime where the ν R were initially in thermal equilibrium with the SM until they decoupled at T FO before the ν L decoupling one would find [71,109] ∆N eq. eff. = 2 · g ν R 2 10.75 g * S (T FO ) 4 3 (4. 18) instead. Integrating the collision term in (4.15) is straightforward and we find (4.20) Our estimate for the additional number of relativistic species is in the limit T RH T .
As expected the abundance of non-thermal DR strongly depends on their production temperature T RH . Note that while it seems that the above expression can lead to arbitrarily large values of ∆N eff. one should keep in mind, that the present treatment relying on (4.15) breaks down as soon as one starts to violate (3.35) because the ν R thermalize. In that case one can use (4.18) to compute ∆N eff. from the freeze-out temperature and finds that it asymptotes to a value of two for two ν R . By plugging in the lower limit on v B-L from the DM production being out of thermal equilibrium in (3.35) we find that the freeze-in contribution of ν R via Z mediated scatterings is at least a factor of five below the sensitivities of the upcoming CMB experiments ∆N eff. < 1.2 × 10 −2 · 85 g * S (T RH = 5 GeV) 4 3 . (4.22) We conclude that the interplay of the tiny rates ∼ v −4 B-L together with the fact that we consider a cosmology with a low reheating temperature reduces the impact of ν R and χ on ∆N eff. below all current and future sensitivities. This opens up an interesting indirect way to test our model: Should observations ever point to ∆N eff. > 0.012 our scenario for DM production is excluded. For the numerical evaluation of (4.15) we proceed as in section 3.4. The temperature dependence of ∆N eff. was depicted in 5 and 6 together with the limit from Planck [71]. For better visibility of the final DR yield we chose values of v B-L below the bound (3.35). The curves in 5 demonstrate that the abundance strongly depends on the reheating temperature and 6 that it decreases with growing v B-L . Both plots show how the final yield is reached shortly after reheating as was the case for DM production. There is also a contribution to the annihilations of SM fermions to ν R via the exchange of an SM like Higgs. The corresponding rate density reads in terms of the coupling (3.41) . (4.24) The above was evaluated at T RH = 1 GeV, where all charged fermions except the top and bottom quark contribute. Figure 7 demonstrates the available parameter space for realizing the entire DM abundance from χs via freeze-in together with the predicted amount of dark radiation parameterized in terms of ∆N eff. . A few comments are in order: The gray region excluded by (3.35) has a more rugged contour because of the sequence of Heaviside functions in the expression for the sum of fermion charges (3.36). Additionally there is a noticeable kink in all DM and DR contours, which occurs around the temperature of the QCD phase transition at T QCD 200 MeV. The physical reason for this behaviour can be  found by inspecting the expressions for the DM and DR yields in (3.38) and (4.19): The integrands in both cases depend on inverse powers of g * S (T ) and g * ρ (T ) and the number of relativistic degrees of freedom in entropy and energy both decrease drastically when the quarks and gluons confine at T QCD . To keep the relic density or ∆N eff. fixed one needs to compensate this increase of the integrand by allowing for a larger value of v B-L , hence the contours appear to be shifted to the right below T QCD , which is why for illustration we chose to display a straight line at the corresponding temperature. One should not forget that the factor of i N i in both numerators also decreases sharply below T QCD , but is approximately cancelled by one of the factors in the denominator leaving one factor in the denominator leading to the previously explained behaviour. It is evident from 7 that the Planck constraint on ∆N eff. would only be relevant for DM masses far below 4 keV, which is already excluded by the Lyman-α constraints. Moreover it is clear that producing ∆N eff. 0.06 only occurs in regions where there is either too much DM or the freeze-in approximation for DM production is not applicable because the production rates from relativistic SM fermions thermalize. Moreover we see that for larger allowed DM masses there is actually less ∆N eff. . The reason is simply that larger m DM at constant T RH require larger v B-L to fix the relic density, which decreases ∆N eff. ∼ v −4 B-L . Consequently our scenario for FIMP DM predicts only a small value of ∆N eff. despite the fact that we introduce two ν R and a rather light DM candidate. This makes the present construction different from the cosmology of other (Dirac) neutrino mass models like e.g. the neutrino-philic Two-Higgs-Doublet model [117][118][119][120][121] or its gauged variations such as [122][123][124][125][126][127] which usually feature light mediators below the EW scale that unavoidably thermalize the ν R and themselves leading to ∆N eff. > O(0.1) [128,129]. Another interesting scheme is called "Common Origin of Warm and Relativistic Decay Products" (COWaRD) [130], where DM and DR are produced together from the decay of a parent particle and the amount of ∆N eff. is correlated with the warmness of DM. There a non-zero ∆N eff. can help to reduce the σ 8 -tension for large scale structures [131,132]. In a sense the COWaRD scheme is the opposite of our idea as it involves thermal DM and predicts a larger amount of DR. All of these models have in common that the more stringent limits on ∆N eff. will already constrain significant amount of their parameter space or even exclude them completely. The only ways to exclude our scenario would be CMB experiments in the far future with a sensitivity to even smaller values of ∆N eff. = O(10 −3 ) or the actual observation of a signal with 0.28 > ∆N eff. > 0.012, which by itself would be a smoking gun for different BSM physics.

Inflation and candidates for the inflaton
The assumed production mode for DM crucially relies on a low value of the reheating temperature 4 MeV T RH 5 GeV together with the assumption of no primordial DM abundance from e.g. inflaton decays during reheating. This puts non-trivial constraints on the explicit realization of inflation. Of course one can assume that the scalar field responsible for creating the inflationary phase of cosmic expansions is another scalar field with no couplings to the DM. However the present model already contains four different scalar multiplets so a minimal solution is to embed the inflaton into one of them. For concreteness we will assume that the candidate field for inflation is the real component of a complex scalar field ω. Recent Planck constraints [133] disfavour monomial inflation of the form Re(ω) p with p > 1 because their potential is too steep leading to a too large tensor to scalar ratio. This is why we only investigate scenarios with a non-minimal coupling of the inflaton to gravity: This scenario is known as Starobinsky-like inflation [134][135][136][137][138][139][140][141] and Figure 7. We depict the allowed combinations of the reheating temperature T RH and the scale of B-L breaking v B-L . The blue shaded area indicates where DM would overclose the universe and the blue contours reproduce the observed DM relic density for m DM ∈ 4, 16, 100, 10 3 keV. Furthermore we show the contours for generating ∆N eff. within the Planck bound [71], the estimated sensitivities of the South Pole Telescope [105], the Simons observatory [106] and for the CMB stage 4 experiment [107][108][109] as well as PICO [110] . The grey area is excluded because the interaction producing DM would equilibrate see (3.35) and searches fom LEP exclude v B-L < 6.9 TeV [31]. the action in the Jordan frame reads In this context we denote the determinant of the metric as g, the Ricci curvature scalar as R and ξ ω is a dimensionless coupling. On can single out a scalar ω field to play the role of the inflaton by imposing that the couplings of the other scalar fields satisfy λ ω /ξ 2 , where the λ denote the scalar self couplings. The remaining fields will be treated as spectator fields. We will use the constraints from reheating to find the appropriate inflaton candidate in our model. Due to the presence of additional scalars besides the inflaton there is the possibility of creating isocurvature perturbations as in multi-field inflation models [143,144], which could come into conflict with CMB bounds. Essentially the problem is that massless particles are sensitive to quantum fluctuations during inflation [145]. However large isocurvature fluctuations can be prevented if either the tree-level mass or the effective mass generated from inflaton oscillations during reheating is larger than the Hubble rate during inflation [146]. Since both η, σ have tree-level masses unconnected to any vev and potentially receive effective masses, we do not expect isocurvature perturbations in these directions. Similarly if we assume that the scale of B-L breaking v B-L is larger than the Hubble rate H I during inflation and U(1) B-L is never restored, then the would-be-Goldstone mode ϕ I corresponds to the longitudinal mode of the massive Z and not to a massless field. Reference [142] found that in the extension of the SM with an inert doublet η housing the inflaton there are only negligible isocurvature fluctuations. A detailed investigation of these fluctuations for the full model is beyond the scope of the present study and we will be content with just outlining how inflation could be realized. Note that we can also allow for a temperature at the end of inflation far above the MeV and GeV range if there is an additional long lived particle that dominates the energy budget of the universe. This leads to an intermediate matter dominated phase [65] which can end in a second radiation dominated epoch with a smaller temperature of the required order of magnitude.

The SM like Higgs
Using the SM like Higgs as the inflaton [139][140][141][147][148][149][150] is a very minimal scenario see [151] for a review. The main drawback of this approach is that the measured value of the Higgs self coupling λ H requires a rather large value of ξ H = O 10 4 , which might give rise to unitarity problems [152][153][154] at scales above M Pl /ξ H . The unitarity problem could for instance be cured by assuming a different coupling to gravity [155][156][157]. Another possibility is exploiting that the SM Higgs self coupling λ H becomes very small at large energy scales which flattens the potential and leads to ξ = O (10), which is known as critical Higgs inflation [158][159][160][161][162]. In terms of BSM physics there is also the attractive possibility to invoke additional scalars to modify the Higgs potential see e.g. [163][164][165]. While it would be interesting to see whether the additional scalars present in this model can solve the unitarity problem it would definitely require a dedicated analysis beyond this work. Consequently we do not consider Higgs inflation further and investigate the other scalar fields as inflaton candidates.

The B-L breaking singlet
The only singlet with a B-L breaking vev could be the inflaton too [166,167]. We neglect the mixing between h and ϕ because the EW gauge symmetry is restored at large temperatures [168,169] so the mixing term vanishes together with v H . For the same reason we compute the decays to the entire doublet H and not just h. For the purpose of finding estimates we work in the regime of perturbative reheating. We assume that all additional scalars, fermions and the Z are heavier than the inflaton so the only available decay modes are where the decay width to DM is obtained from (3.3) after converting it to ∼ m DM /v H and replacing v H → v B-L . Since the scalar will oscillate in its potential during reheating it develops an effective mass depending on its oscillation frequency and the same goes for all other scalar fields as well as the Z since they share quartic couplings with ϕ. Hence requiring that (5.2) are the only available decay modes and that e.g. ϕ → S 1 S 2 h is absent amounts to a bound on the effective field dependent masses and not on the tree-level masses that we have employed so far. For the sake of simplicity this first estimate will work exclusively with the tree level masses. If we want to avoid a primordial abundance of DM the first step is to make sure that decays to SM particles dominate the reheating process which sets bounds on the model parameters. Assuming BR 1 we can determine the reheating temperature from the decay to the SM Higgses, which themselves will decay to fermions creating a hot thermal bath. In this limit the reheating temperature is found to be The assumed range of reheating temperatures for DM production requires that either λ Hφ 1 or that m ϕ v B-L . However the second condition can not be realized because m ϕ is proportional to v B-L according to (2.41) and we can not make m ϕ arbitrarily heavy due to the perturbativity limit λ φ < √ 4π. Inflaton decays can produce DM as well and the corresponding Boltzmann equation during reheating reads [170] dn χ dt and we denote the energy density of the non-relativistic inflaton condensate as ρ ϕ . The DM yield today is found to be [170] Y χ (T 0 ) 3 4 T RH m ϕ BR (5.6) and the DM energy density today can be calculated with (3.20). For simplicity we assume g * ρ (T RH ) g * S (T RH ). We trade the inflaton mass via equation (5.4) for an expression involving T RH and v B-L , where the dependence on λ Hφ in Y χ (T 0 ) divides out. By using our limit on v B-L in (3.35) we derive an upper-limit on the relic abundance from inflaton decays It is evident that large m DM and low reheating temperatures could lead to an abundance that is larger than the FIMP one in (3.39). Demanding that the abundance from inflaton decays does not overclose the universe cuts away all the available parameter space in 7.
There is another reason why this channel is not suited for light DM production: Since the production mode is different from both freeze-in (which requires a thermal bath) and thermal production, the phase space distribution and hence the velocity distribution of the DM will be different assuming all of DM was produced via this single channel. This manifest itself in a modified Lymann-α bound [170,171] m DM 2 keV · m ϕ T RH · m WDM 3.5 keV , (5.8) which was recast from the bound for thermally produced DM with m WDM 3.5 keV [172] (which is the average of the two possible warm DM masses in section 3.1). If we assume that m ϕ = O (v B-L ) then we expect an inflaton with at least a TeV scale mass (see (2.36)), which is much larger than the assumed MeV-GeV reheating temperatures. Therefore the DM mass for inflaton production would be orders of magnitude larger than 2 keV and potentially violates the bound from invisible Higgs decays in (3.8) and could lead to overclosure. Thus we conclude that for our purposes ϕ can not be the inflaton, because it tends to produce too much DM. Therefore we assume that ϕ is too heavy to be produced during reheating.

The inert doublet or singlet scalars
As previously mentioned v H vanishes due to the restoration of the EW symmetry at large temperatures [168,169] . In this limit we can relate S 1 = η 0 R as well as S 2 = σ 0 R and consider each field as a candidate individually. Similar to Higgs inflation the inert doublet η can house the inflaton [142]. This scenario is free of the unitarity problem because the value of the η self coupling λ η is unconstrained by phenomenology. We can not just reuse the perturbative reheating estimate (5.4) from the previous section, because without a vev there is no tree level decay to Higgses like in (5.2) or to EW gauge bosons for η 0 R . In this model reheating occurs via quartic couplings to electroweak gauge bosons and SM Higgses [141,142] and we assume that the Z is too heavy to be produced. Reheating typically takes place through resonant gauge boson production which then annihilate to SM fermions. In this scenario the reheating temperature was found to be [142] T η RH 10 14 GeV λ Generating sub-GeV reheating temperatures is impossible in this regime, because it would require non-perturbative values of λ η . We conclude that another reheating channel is needed and hence consider an inflaton without SM gauge interactions: σ is an SM singlet and has no vev as well. If we assume that the effective field dependent mass of the Z is too large to be produced then creating SM Higgses via the quartic coupling λ Hσ in (2.6) is the only possibility left. Since this process depends on the new coupling λ Hσ instead of the known SM gauge couplings the reheating temperature will also depend on this unconstrained parameter. Subsequent decays and annihilations of the Higgs to SM states then seed the SM radiation bath. Reference [173] found that for resonant Higgs production T σ res.

RH
3 × 10 13 GeV λ σ λ 2 Hσ 1 4 . (5.10) The analysis [173] made the conservative assumption of having reheating occur during the quadratic phase of the potential before the quartic self-interaction of the inflaton becomes dominant, which can be expressed as λ σ > 0.25 λ Hσ [173]. If we drop this assumption, which [173] emphasizes is not ruled out, we can choose smaller values of λ σ λ Hσ and can at least in principle accommodate the range 4 MeV T RH 5 GeV. The authors of [173] also found that reheating can occur in another regime if inflaton excitations annihilate into pairs of Higgs bosons leading to the estimate T σ ann. RH 9 × 10 13 GeV · λ 1 4 σ . (5.11) The conservative assumption about reheating occurring in the quadratic regime of the potential would lead to λ σ > 0.019 [173], but again we need to drop this assumption and require λ σ 1 to obtain the phenomenologically favoured reheating temperatures. In the next section 6 we will introduce a decay of σ to exotic quarks, which might open up another possibility for realizing the required reheating temperature. Let us emphasize that there are bounds from vacuum stability and perturbativity on the quartic couplings [174], but since these bounds are usually obtained in models with a simpler scalar sector it requires a dedicated study to translate them to our construction, because of e.g. threshold effects from heavy scalars [35]. Note that at some point during reheating there will be the SSB of the EW symmetry generating a coupling of σ 0 R to the EW gauge bosons proportional to sin(α). But since the neutrino mass (2.15) does not directly depend on the mixing angle α in the radiative seesaw limit we can make this mixing small. If we assume that the F fermions are heavier than the σ there will be no inflaton decays to χ via the Yukawa interaction in (2.22). The only way to generate the unwanted primordial DM population would be annihilation processes of the form σ 0 R σ 0 R → χχ mediated by heavy F s. We do not expect this to lead to a significant DM abundance, because scattering is inefficient for non-relativistic excitations of the inflaton field and the production is suppressed by the heavy F mass. On top of that the DM production competes with the unsuppressed process for creating the SM radiation σ 0 R σ 0 R → H † H. Since the singlet scalar might not have decay modes, we need to ensure that the inflaton becomes a subdominant component of the universe's energy budget after reheating. The additional interactions like Higgs or Z mediated scatterings with the SM fermions could help thermalize the inflaton with the radiation bath, which is already in thermal equilibrium [72]. We conclude that the only possible inflaton candidate that is not in conflict with the cosmological DM and reheating requirements is σ 0 R .

Baryogenesis
The assumed low scale reheating is hard to reconcile with most known mechanisms [175][176][177] for Baryogenesis. Leptogenesis [178] for instance relies on producing a leptonic asym-metry that gets converted into a baryon asymmetry by electroweak sphaleron processes, which are in equilibrium only above the EW phase transition at T EW = O(100 GeV). On top of that since the SM neutrinos do not mix with any of the heavy new neutrinos N, F we can not realize leptogenesis via oscillations [179] as well. Thus we are left with mechanisms that do not rely on the sphaleron transition above the EW scale. One example is the spontaneous Baryogenesis [180,181] mechanism, which however needs reheating temperatures far above the assumed MeV-GeV scale window. Hence some other form of non-thermal baryogenesis during reheating seems to be the only possibility left if we insist that the temperature at the end inflation is indeed in the previously mentioned range. The Affleck-Dine mechanism [19] relies on baryon number charged scalars whose real and imaginary parts evolve non-trivially in time, which acts as a source term for baryon number. This scenario can in principle operate at low reheating temperatures if the initial field value of the Affleck-Dine field is very large compared to its mass. Since all of our scalars except H are charged under B-L this is an attractive possibility. For concreteness we will treat σ as the Affleck-Dine field; whether it can accommodate both baryogenesis and inflation at the same time like e.g. [146,[182][183][184][185][186][187][188][189][190][191][192] will be left for future investigation. An important ingredient is a small explicit Baryon number breaking interaction. Of course we can not break our gauged B-L explicitly but a term of the form λ AD (σ 4 + σ * 4 ) could arise after the spontaneous breaking of U(1) B-L . To do so we allow for the small Z 5 breaking term L ⊂ −λ σ 2 φ 2 + h.c. , (6.1) which after integrating out the heavy radial mode ϕ (we ignore the ϕ-Higgs mixing here) leads to an operator and we can identify λ AD = λ 2 v 2 B-L /m 2 ϕ λ 2 /(2 λ φ ) from (2.41). Quite interestingly this allows us to make λ AD small just by assuming √ λ λ φ . Since λ → 0 would restore the discrete symmetry the choice λ 1 is technically natural [193]. Of course assuming the existence of this operator begs the question why the other Z 5 breaking interactions are absent. The last missing ingredient is a way to transmit the σ-asymmetry to the quarks. To do so we introduce a pair of heavy vector-like quarks (Q L , Q R ) that are weak isospin singlets with the hypercharge Y = −2/3 (4/3) of the right chiral down (up) quarks. The quarks come with a B-L charge Q σ + 1/3 = −2/3 and transform as ω −4 under Z 5 , where ω = e 2iπ 5 , so that we can realize the operators Here d R can in principle also be replaced with u R ; we chose the hypercharge −2/3 to make the vector-like quarks resemble the down-type quarks which might help with unification [41,42]. The above interaction could also lead to inflaton mediated washout scatterings depleting the baryon asymmetry [189], which puts constraints on the coupling Y Qq . In order to prevent stable exotic quarks from forming relics [194] we have to demand that m Q > m σ so the Q can decay via the above operator to σ u R in the late universe. Alternatively one can also arrange for m σ > m Q > 2m h instead so that the decay of the vector-like quarks proceeds via off-shell σ as Q L → σ * + d R → 2h + d R ; the Higgses then further decay to SM states. In the early universe the field σ receives a potentially large effective mass from inflaton oscillations during reheating so for both aforementioned cases the CP -conserving decay σ → Q L u R would be possible and one can indeed transmit the asymmetry from the Affleck-Dine field to the quark sector. This decay could open up another interesting reheating scenario as well. In the following we will assume that T RH arises either from to the channels enumerated in the previous section 5.3 or via the aforementioned decay. An estimate for the baryon asymmetry leads to [195,196] n B s 10 −10 · λ AD 10 −2 · sin (4θ i ) 0.5 · r i /m σ 6 × 10 6 3 · r i 6 × 10 9 GeV · T RH 1 GeV . (6.4) Here we use the polar parameterization for σ, where r i is the initial value of the radial component and θ i denotes the initial angle needed for CP violation. This decomposition should not to be confused with the cartesian representation from (2.7). The initial angle can not be set to arbitrarily small values in order to avoid isocurvature perturbations [190], which is why we chose sin (4θ i ) 0.5. We see that very large initial field values are needed to compensate for the low reheating temperature. Such a high field value of r i /m σ 6×10 6 usually requires a very flat potential and could be an initial condition. Alternatively the non-minimal coupling to gravity might help to generate this field value dynamically [197]: It was found that this coupling together with the tree level mass squared creates an effective mass squared depending on the Hubble parameter. This effective mass is tachyonic at early times when H m σ and later turns real again, which can be understood as an inverted phase transition [198,199]. Afterwards the field, which can be visualized as an over-damped oscillator, is stuck in its previous non-trivial minimum corresponding to an initial value of [68] r i ξ σ λ σ m σ , (6.5) before it starts to relax to its true minimum σ = 0 as soon as the Hubble rate satisfies H ∼ m σ provided that λ AD λ σ . From this mechanism we can deduce that a scalar self coupling of λ σ 2.8 × 10 −14 · ξ σ · r i 6 × 10 9 GeV 2 · 1 TeV m σ 2 (6.6) would be required for the initial field value and a scalar mass in accord with our previous estimates (2.16) and (2.31). Note that this violates the previous assumption λ AD λ σ , but we can reconcile this by assuming that the heavy ϕ will only be integrated out at temperatures somewhat below the inverted phase transition so that the operator (6.2) is absent initially. On the level of estimates it seems that our scalar potential can reproduce the observed baryon asymmetry, but again we stress that it requires a separate study to work out the details especially in the inflationary context and considering the radiative stability. It is noteworthy that the operator (6.1) also sources a mass splitting ±λ v 2 B-L between the real and imaginary parts of σ, while for the neutrino mass generation we assumed that they are mass degenerate. Under the assumption that this additional mass splitting is small compared to the overall mass scale of the S 1,2 (A 1,2 ) and the mass splitting between the different generations of scalars our conclusions about the neutrino and DM masses are unchanged. If T RH is the temperature after an intermediate epoch of matter domination and the true temperature at the beginning of the first radiation dominated phase was far above the electroweak scale this allows for the other previously discussed mechanisms again. In that case the challenge is to generate enough entropy to dilute unwanted relics (such as thermally produced DM) while retaining enough baryon asymmetry [200].

Summary
We presented an extension of the Dirac scotogenic model [5,6] that creates the Dirac mass of a light fermionic DM candidate χ together with the active neutrino masses via one-loop diagrams. The model relied on a gauged U(1) B-L symmetry, whose anomaly-freedom determined the charges of the DM and two copies of ν R . We found that our symmetry based approach predicts that only two SM neutrinos are massive Dirac fermions, whereas the third one remains exactly massless, because there is no third ν R . In order to ensure the DM stability and to prevent unwanted operators that could affect the neutrino or DM mass generation we had to impose a separate Z 5 symmetry as well. Additionally one requires an inert scalar doublet η and an inert singlet σ together with the B-L breaking scalar singlet φ. Moreover we had to introduce a host of vector-like fermions to generate the necessary loop diagrams. It was found that the vector-like leptons F needed for the DM masses couple to the SM Higgs and are light enough to potentially be probed by next generation collider experiments. We then chose a minimal scenario where we assumed that only the SM degrees of freedom augmented by two ν R and χ are present after reheating. The constraint from invisible Higgs decays enforces m DM 2 GeV and the DM mass has to be larger than (4 − 16) keV due to the Lyman-α forest. After demonstrating that thermal production and out of equilibrium Higgs decays both lead to an over-production of DM, we were able to narrow the window of the allowed reheating temperatures down to the range between about 5 GeV and 4 MeV. Consequently we analyzed the joint production of DM χ and DR ν R from out of equilibrium annihilations of the SM fermions via the B-L gauge boson Z . The DM mass has to be smaller than O(MeV) in order to suppress DM production via diagrams with an intermediate SM like Higgs compared to Z mediated scatterings. We found a potentially viable parameter space with v B-L O (10 TeV) that leads to the correct observed DM abundance but predicts ∆N eff. 0.012. The amount of produced dark radiation decreases with the DM mass so in a sense m DM and ∆N eff. are anti-correlated. This is in striking contrast to other Dirac neutrino and DM mediator models which usually predict larger ∆N eff. . Thus while the aforementioned models can already be tested or ruled out by tightening the observational bounds on ∆N eff. , only the detection of ∆N eff. > 0.012 could falsify our DM production scenario in the near future.
Owing to the fact that we need a very low reheating temperature and want a negligible primordial DM abundance we were able to single out the real component of the σ field to play the role of the inflaton. In addition we found a way for how the σ field can also potentially realize Affleck-Dine baryogenesis if we introduce a small source of Z 5 -breaking in the scalar potential together with a pair of vector-like down quarks. We leave a detailed study of the inflationary predictions, reheating and non-thermal baryogenesis for future investigation. To summarize, we introduced a new abelian gauge theory that can simultaneously explain the active neutrino and fermionic dark matter masses via loop diagrams. Our construction produces the observed DM relic abundance together with minuscule amounts of dark radiation in the freeze-in regime and can potentially account for inflation, reheating and Baryogenesis.