Probing the Scotogenic Dirac Model with FIMP Dark Matter and ΔN_eff

The phenomenon of neutrino oscillations [1−6] provides strong evidence that neutrinos must have non-zero masses. For neutral fermions such as neutrinos, two types of mass are theoretically possible: Dirac and Majorana. On the experimental front, extensive efforts have been devoted to searching for lepton number violation—a featured signature of Majorana neutrinos. Experiments at colliders [7, 8] and neutrinoless-double-beta-decay ($ 0\nu\beta\beta $) [9−14] searches, in particular, have so far yielded no positive signals. As a result, the alternative scenario has received increasing attention in recent years [15−54]. To generate the tiny Dirac neutrino masses, an attractive approach is to attribute them to a radiative mechanism realized through loop diagrams. The loop suppression naturally explains the smallness of the masses, while the particles running in the loop can carry quantum numbers under a discrete symmetry. This same symmetry can stabilize the lightest particle in the loop, making it a viable dark matter (DM) candidate.

The weakly interacting massive particle (WIMP) paradigm, featuring electroweak-scale particles with weak cross sections ($ \sigma \sim10^{-26} $ cm$ ^3 $/s), naturally yields the observed DM density ($ \Omega_{\rm DM}h^2\approx0.12 $) via thermal freeze-out (see [55] for a review). However, this paradigm faces significant challenges from null results in DM direct detection [56, 57]. In response, the feebly interacting massive particle (FIMP) scenario has emerged as a compelling alternative [58−62]. Unlike thermal WIMPs, FIMPs are produced via extremely weak couplings, avoiding current experimental bounds while naturally explaining the observed dark matter abundance.

The Scotogenic Dirac Model, originally proposed in Ref. [20] and later developed with an alternative symmetry realization in Ref. [63], extends the SM with right-handed neutrinos $ \nu_R $, vector-like fermions N, a scalar doublet Φ, and a scalar singlet χ. Previous work in Ref. [63] explored the scenario where the lightest fermion $ N_1 $ acts as a WIMP dark matter. This scenario faces significant challenges due to stringent constraints from lepton flavor violation (LFV) experiments, which severely limit the Yukawa couplings required for efficient thermal relic production. Consequently, achieving the observed relic density requires substantial tuning of parameters and relies on annihilation channels mediated by the $ y_\chi $ coupling. In this work, we propose a novel realization of the Scotogenic Dirac Model in which the lightest neutral fremion $ N_1 $ serves as a FIMP DM candidate. This scenario is defined by the extreme suppression of the Yukawa couplings $ (y_\Phi)_{\alpha 1} $ and $ (y_\chi)_{\alpha 1} $, which ensures that $ N_1 $ never enters thermal equilibrium with the Standard Model (SM) plasma throughout cosmic history. This FIMP framework not only provides a natural mechanism for DM production but also has a profound impact on neutrino physics, leading to a rank-2 structure for the neutrino mass matrix. This structure predicts one nearly massless neutrino state, which is allowed by the oscillation data. We systematically investigate two complementary production mechanisms for $ N_1 $. The first is the conventional freeze-in mechanism, where $ N_1 $ is slowly produced via the extremely small Yukawa couplings through decays of heavier $ Z_2 $-odd scalars and fermions. The second is the super-WIMP mechanism, where $ N_1 $ is produced by the late decays of a thermally frozen-out next-to-lightest odd particle (NLOP). The identity of the NLOP—whether it is $ N_2 $, $ \phi_{R,I} $, $ \phi^\pm $, or χ—plays a crucial role in shaping the model's phenomenology, directly influencing the DM relic abundance and cosmological observables. A central focus of our analysis is the effective number of relativistic species, $ \Delta N_{{\rm{eff}}} $. We show that the deviation $ \Delta N_{{\rm{eff}}} $ receives contributions from two distinct sources: a thermal component from the primordial bath of right-handed neutrinos, and a non-thermal component arising from energy injection due to late NLOP decays. Our comprehensive parameter scan demonstrates that viable regions of parameter space exist for all possible NLOP candidates, satisfying all current observational constraints from DM relic density, lepton flavor violation, BBN, CMB, and $ \Delta N_{{\rm{eff}}} $.

The remainder of this paper is organized as follows: Section II introduces the FIMP Scotogenic Dirac Model and the mechanism for Dirac neutrino mass generation. Section III details the dark matter relic density calculations. Section IV analyzes the cosmological implications for the effective number of relativistic species. Finally, we make our conclusion in Section V.

The Scotogenic Dirac Model was originally proposed in Ref. [20], and a different symmetry realization was put forward in [63] by one of the authors. In Table 1 we list out the relevant particles in this model. Apart from the Standard Model lepton doublet $ F_L = \left(\nu_L, \ell_L\right)^T $ and Higgs H, additional species are the right-handed neutrinos $ \nu_R $, three vector-like fermions N, an extra scalar doublet $ \Phi\equiv \left(\phi^+, (\phi_R + i\phi_I)/\sqrt{2}\right)^T $ and scalar singlet χ. Tree-level Dirac neutrino masses and possible Majorana masses of $ \nu_R $ and N are forbidden due to the existence of $ Z_3 $.

	$ F_L $	H	$ \nu_R $	N	Φ	χ
$ SU(2)_L $	2	2	1	1	2	1
$ U(1)_Y $	$ -1/2 $	$ 1/2 $	0	0	$ 1/2 $	0
$ Z_3 $	0	0	ω	ω	ω	0
$ Z_2 $	+	+	+	-	-	-

Table 1.  Particles and symmetries in the scotogenic Dirac model.

The introduction of new particles in Table 1 would bring in the following interactions:

where $ \tilde \Phi \equiv i \sigma_2 \Phi^\star $. Relevant scalar potential terms are given by

The scalar singlet χ is set as real, for simplicity. The μ-terms in Eq. (2) is the source which softly break the $ Z_3 $ symmetry, hence it's natural to keep it a small value. The one-loop neutrino mass (as illustrated in Fig. 1) is then generated as

Figure 1.  The one-loop neutrino mass in Scotogenic Dirac Model.

Lower indices $ \alpha(\beta) $ and k stand for different generations of leptons and Ns. The angle θ in Eq. 3 originates from mixing between the real component of scalar doublet Φ and singlet χ. When the SM Higgs H acquires its vacuum expectation value v after electroweak spontaneous symmetry breaking, the angle is quantitatively fixed as

The $ m^2_{1(2)} $ in Eq. 3 indicate eigenvalues of the $ \chi-\phi_R $ mixing mass-square matrix

Treating μ as a tiny parameter results in small mixing, justifying the approximation $ m_{1(2)} \approx m_{\chi(\phi_R)} $. Hence we will use $ m_{\chi} $ and $ m_{\phi_R} $ hereafter. The masses of the other scalars are $ m_h^2 = 2 \lambda_1 v^2, \; m_{\phi^\pm}^2 = \mu_\Phi^2 + \lambda_4 v^2, \; m_{\phi_I}^2 = \mu_\Phi^2 + (\lambda_4 + \lambda_7) v^2 $. Here h is the observed 125 GeV boson.

The assignment of $ Z_2 $ in Table 1 would make particles in loop of Fig. 1 to be in dark-sector, and the lightest one is naturally served as the dark matter particle candidate. In Ref. [63], we had made a detailed study on the situation that $ N_1 $ is a WIMP-type dark matter particle. In this work, however, we explore an alternative scenario where the dark matter candidate $ N_1 $ is a feebly interacting massive particle. The FIMP nature of $ N_1 $ implies that the Yukawa couplings $ (y_\Phi)_{\alpha 1} $ and $ (y_\chi)_{\alpha 1} $ are extremely small. This setup significantly distinguishes the model from the conventional WIMP Dirac scotogenic framework. From Eq. 3, we observe that the neutrino mass matrix $ M_\nu $ is approximately of rank-2, a direct consequence of the suppressed couplings $ (y_\Phi)_{\alpha 1} $ and $ (y_\chi)_{\alpha 1} $. This structure predicts that one of the three neutrino mass eigenstates is nearly massless, which remains consistent with the oscillation data requiring at least two massive neutrinos.

The Yukawa interactions $ y_\Phi $ could mediate processes leading to LFV, which is extremely suppressed in SM. On the experimental side, the most stringent limits arise from the rare decays $ \ell_\alpha \to \ell_\beta \gamma $. The decay branching ratio is calculated as

here $ j(r) = \left(1-6r + 3r^2 + 2r^3 -6r^2\ln r\right)/\left(12(1-r)^4\right) $. The latest result from MEG II reported the limit as $ \text{Br}(\mu \to e \gamma) \lesssim 1.5 \times 10^{-13} $ [64]. This could be translated into limit on $ y_\Phi $ as $ (y_\Phi)_{\alpha 2,3} \lesssim 0.01 $, if we assume the new masses are not far away from the electroweak scale. When combined with the neutrino mass in Eq. 3, and noting that the mixing angle θ is small (as it reflects the $ Z_3 $-breaking parameter μ in Eq. 4), we deduce that $ (y_\chi)_{\alpha 2,3} $ must be comparable to or larger than $ (y_\Phi)_{\alpha 2,3} $. For example, with $ \sin2\theta \sim 10^{-5} $ and $ (y_\Phi)_{\alpha 2,3} \lesssim 0.01 $, reproducing the observed neutrino masses requires $ (y_\chi)_{\alpha 2,3} \gtrsim 0.01 $.

The residual $ Z_2 $ symmetry confines loop particles to the dark sector, with the lightest $ Z_2 $-odd particle serving as dark matter candidate. In [63], we analyzed the scenario where the fermion $ N_1 $ acts as WIMP dark matter. We found that the dark matter annihilation via $ y_\Phi $ is difficult to achieve at the correct relic density. This limitation arises because the Yukawa coupling $ y_\Phi $ faces stringent constraints from LFV. The $ y_\Phi $ compatible with these LFV constraints substantially suppresses $ N_1 $ annihilation cross-sections, resulting in dark matter overproduction that is excluded by cosmological observations. Consequently, the dominant annihilation must proceed through $ y_\chi $.

In this work, however, we would try a FIMP realization of dark matter $ N_1 $. The feeble interaction strength would retain $ N_1 $ always out of equilibrium. Thus the key ways to produce $ N_1 $ are decays of heavier $ Z_2 $-odd particles. The dominant contributions are from decays of scalars:

Additional production channels involve decays from $ N_{2,3} $:

Fermionic decay contributions remain subdominant due to both three-body phase space suppression and additional Yukawa coupling factors. Moreover, decays into $ \bar{N}_1 $ and decays of $ \bar{N}_2 \to N_1(\bar{N}_1) $, whose widths are similar as in Eqns. 7 to 11, should also be taken into account to get the total abundance of dark matter. To ensure $ N_1 $ remains out of thermal equilibrium, the decay rate should be smaller than the Hubble expansion rate. A rough estimate requires the Yukawa coupling $ (y_{\Phi,\chi})_{\alpha 1} \lesssim 10^{-7} $, for scalar masses at $ {\cal{O}}(\rm TeV) $.

The evolution of dark matter and relevant particles are governed by the following Boltzmann equations:

Here $ X = \phi_R, \phi_I, \phi^\pm , \chi, N_{2,3}(\bar N_{2,3}) $ stands for all the $ Z_2 $-odd particles that could decay into dark matter. For a specific X, its evolution is affected by three kinds of processes: decays of heavier A into X, decays of X into lighter B, and SM pairs annihilating into X pairs. The dimensionless parameter z is defined as the ratio of the dark matter mass to the evolution temperature, i.e. $ z \equiv m_{N_1} / T $. The parameter $ k^{*} $ is given by $ k^{*}=\sqrt{45/(4 \pi^3 g^*)} M_{\rm Pl} /m_{N_1}^{2} $, $ g^* $ is the effective number of degrees of freedom of the relativistic particles, $ M_{{\rm Pl}} = 1.2 \times 10^{19}\; \rm{GeV} $ is the Planck mass. Thermal decay width $ \tilde{\Gamma}_{i\to j} $ is defined as $ \tilde{\Gamma } _{i\to j} =\Gamma_{i\to j} K_1/K_2 $, where $ K_1 $ and $ K_2 $ denote the first and second order of modified Bessel functions of the second kind, respectively. $ Y_i^{\text{eq}} $ represents the abundance in equilibrium for species i.

In practice, we write down the model using $ \mathrm{FeynRules} $ [65] and the output model files are used to calculate relic density by $ \mathrm{micrOmegas} $ [66]. We present the results of this calculation in Fig. 2, which illustrates the comoving yield of $ N_1 $ as a function of the cosmological variable $ z = m_{N_1}/T $ for several benchmark values of the feeble Yukawa couplings $ (y_\Phi)_{\alpha 1} $ and $ (y_\chi)_{\alpha 1} $. For simplicity in this analysis, we set these couplings to be equal, i.e., $ (y_\Phi)_{\alpha 1}=(y_\chi)_{\alpha 1} $, while the remaining model parameters are fixed to the values specified in Table 2. The resulting evolution curves exhibit the characteristic hallmark of the freeze-in mechanism. At high temperatures ($ z \ll 1 $), the abundance of $ N_1 $ is negligible. As the universe cools, $ N_1 $ is slowly but steadily produced from decays of the heavier $ Z_2 $-odd particles. The yield grows progressively until the source particles become thermally suppressed and depleted, at which point the comoving number density of $ N_1 $ "freezes in" and remains constant. For reference, the observed dark matter relic density, $ \Omega_{\rm DM} h^2 \approx 0.12 $, is indicated by a horizontal dashed line. We demonstrate a viable parameter point, corresponding to the specific coupling strength $ (y_\Phi)_{\alpha 1} = (y_\chi)_{\alpha 1} = 5 \times 10^{-12} $, successfully saturates the observed dark matter density.

Figure 2.  (color online) The evolution of dark matter through freeze-in. The parameters used are listed in Table 2.

$ m_{N_1} $	$ m_{N_2} $	$ m_{N_3} $	$ m_{\phi_R} $	$ m_{\phi_I, \phi^\pm} $	$ m_\chi $	$ (y_{\Phi,\chi})_{\alpha 1} $	$ (y_{\Phi})_{\alpha 2,3} $	$ (y_{\chi})_{\alpha 2,3} $	$ \lambda_5 $	Figure
$ 10 $	$ 800 $	$ 800 $	$ 800 $	$ 800 $	$ 800 $	$ 5\times10^{-10,-11,-12} $	$ 10^{-3} $	$ 10^{-2} $	$ 10^{-3} $	Fig. 2
$ 10 $	$ 600/770 $	$ 800 $	$ 800 $	$ 800 $	$ 800 $	$ 5\times10^{-13} $	$ 10^{-2} $	$ 10^{-2} $	$ 10^{-3} $	Fig. 3a
$ 10 $	$ 300 $	$ 800 $	$ 800 $	$ 800 $	$ 800 $	$ 5\times10^{-13} $	$ 10^{-2} $	$ 0.1/0.3 $	$ 10^{-3} $	Fig. 3b
$ 500 $	$ 800 $	$ 800 $	$ 750 $	$ 800 $	$ 800 $	$ 5\times10^{-13} $	$ 10^{-3} $	$ 10^{-2} $	$ 10^{-3} $	Fig. 3c
$ 10 $	$ 800 $	$ 800 $	$ 800 $	$ 800 $	$ 550 $	$ 5\times10^{-13} $	$ 10^{-3} $	$ 0.5/10^{-3} $	$ 10^{-3}/0.02 $	Fig. 3d

Table 2.  Parameter sets used in Fig. 2 and Fig. 3, masses are given in GeV.

In addition to the conventional freeze-in mechanism, dark matter can also be produced via the late decay of a thermally frozen-out next-to-lightest odd particle (NLOP), a scenario known as the "super-WIMP" mechanism [67]. In this framework, heavier $ Z_2 $-odd particles decay predominantly into the NLOP in the early universe, as these transitions are not suppressed by small couplings. The NLOP is initially in thermal equilibrium with both SM and other $ Z_2 $-odd particles, but as the universe cools, its interaction rate drops below the Hubble expansion rate, leading to decoupling and the freeze-out of its comoving number density. Subsequently, the NLOP undergoes a highly suppressed decay into the stable dark matter, due to an extremely weak coupling. This results in a long-lived NLOP and a delayed, non-thermal production of dark matter that persists until the NLOP population fully decays. The resulting dark matter relic abundance is totally fixed from the NLOP's freeze-out density:

This production mechanism is illustrated in Fig. 3, which showcases several distinct realizations of the NLOP within our model.

Figure 3.  (color online) The evolution of dark matter and NLOP under the "super-WIMP" mechanism. $ N_2 $ is chosen as NLOP in the upper two subfigures, while in the lower two subfigures $ \phi_R $ and χ are selected as NLOP, respectively. The parameters used are listed in Table 2.

Figs. 3a and 3b depict scenarios where the NLOP is the fermionic state $ N_2 $. Its freeze-out is driven by annihilations into SM leptons, mediated by the Yukawa couplings $ (y_\Phi)_{\alpha 2} $ or $ (y_\chi)_{\alpha 2} $. However, as shown in [62], LFV constraints severely limit $ (y_\Phi)_{\alpha 2} $, suppressing the annihilation cross section and leading to an overproduction of $ N_2 $ if this coupling dominates. To achieve the correct relic abundance, additional mechanisms are required. One such mechanism is coannihilation: when $ N_2 $ is nearly degenerate in mass with other $ Z_2 $-odd states, their combined interactions enhance the effective annihilation rate. Fig. 3a compares two benchmarks: $ m_{N_2} = 600\; \rm{GeV} $ and $ 770\; \rm{GeV} $, with the other odd particles fixed at $ 800\; \rm{GeV} $. The heavier $ N_2 $ case benefits from stronger coannihilation effects, yielding a reduced freeze-out abundance consistent with observations. Alternatively, in our model, the coupling $ (y_\chi)_{\alpha 2} $ opens an additional t-channel annihilation via exchange of the scalar singlet χ, producing neutrino final states. As demonstrated in Fig. 3b, increasing $ (y_\chi)_{\alpha 2} $ significantly enhances the annihilation rate, lowering $ \Omega_{N_2}^{\text{freeze-out}} $ and bringing the final dark matter density into agreement with cosmological data.

In scenarios where the scalar doublet Φ serves as the NLOP, we focus on its neutral real component $ \phi_R $ as the representative case (with $ \phi_I $ and $ \phi^\pm $ exhibiting similar dynamics). $ \phi_R $ annihilates efficiently through electroweak gauge interactions (e.g., $ \phi_R \phi_R \to W^+W^-, ZZ $), resulting in strong thermal coupling. The correct relic abundance is only achieved in mass region of $ m_{\phi_R} > 500\; \rm{GeV} $ (or $ m_{\phi_R} < m_W $) [68−71]. In Fig. 3c, we consider a benchmark mass $ m_{\phi_R} = 750\; \rm{GeV} $. A moderate mass hierarchy between $ N_1 $ and $ \phi_R $ is required (hence we set $ m_{N_1} = 500\; \rm{GeV} $), as the freeze-out abundance of $ \phi_R $ is typically not significantly larger than the observed dark matter density [62].

For completeness, we also consider the phenomenology of the charged component $ \phi^\pm $ when it acts as the NLOP. Being electrically charged and long-lived (due to suppressed decays into $ N_1 $), $ \phi^\pm $ behaves as a long-lived charged particle. The ATLAS Collaboration has recently searched for such states in $ 140\; {\rm{fb}}^{-1} $ of $ pp $ collision data at $ \sqrt{s} = 13\; {\rm{TeV}} $, using signatures based on high specific ionization energy loss and time-of-flight measurements [72]. Interpreting their results in the context of stau-to-gravitino decays, ATLAS excludes masses up to $ 560\; {\rm{GeV}} $ for such long-lived charged particles. Furthermore, in Fig. 4 we give the assessment of the projected sensitivity of future long-lived particle detector MATHUSLA [73]. Fig. 4 shows the projected reach in the mass-decay length plane of a long-lived $ \phi^\pm $ [73], assumed to be pair-produced via Drell-Yan processes at the 14 TeV HL-LHC and to decay into dark matter. The reconstruction efficiency is assumed in 0.5–1. Contours indicate regions where at least four events are expected with 3 ab$ ^{-1} $. In Fig. 4 we also show the Yukawa coupling strength as a function of mass and decay length. Nevertheless, for the parameter space that yields the correct dark matter relic density, MATHUSLA's projected coverage remains limited to only a marginal segment.

Figure 4.  (color online) The projected detecting capabilities of MATHUSLA to a long-lived $ \phi^\pm $, in the mass-decay length plane. The band shows the efficiencies vary from 0.5 to 1. Colored lines correspond to different strengths of Yukawa coupling between $ \phi^\pm $ and the dark matter.

Finally, we consider the scenario in which the scalar singlet χ serves as the NLOP. The annihilation of χ in the early universe is governed by three key parameters: the Yukawa coupling $ y_\chi $, the Higgs portal coupling $ \lambda_5 $, and the trilinear scalar coupling μ. Annihilation through $ y_\chi $ proceeds via t-channel exchange of the heavy neutral fermions $ N_{2,3} $, yielding neutrino final states. Alternatively, χ can annihilate into SM particle pairs—such as $ b\bar{b} $, $ WW $, or $ ZZ $—through s-channel Higgs exchange, mediated by the $ \lambda_5 (H^\dagger H) \chi^2 $ interaction. Motivated by the requirement of naturally small neutrino masses and the preservation of an approximate global symmetry in the scalar sector, we treat the coupling μ as hierarchically suppressed throughout this analysis. This naturalness argument renders μ negligible for thermal processes, leaving $ y_\chi $ and $ \lambda_5 $ as the dominant drivers of χ's annihilation dynamics. To illustrate the viable parameter space, we examine two representative benchmark scenarios. The first features a large Yukawa coupling $ (y_\chi)_{2,3} = 0.5 $ with a small Higgs portal $ \lambda_5 = 10^{-3} $, favoring annihilation into neutrinos. The second adopts a larger portal coupling $ \lambda_5 = 0.02 $ while keeping $ (y_\chi)_{2,3} = 10^{-3} $ small, enhancing annihilation into SM states via Higgs mediation. As shown in Fig. 3d, both configurations yield a χ freeze-out abundance that, after subsequent decay into $ N_1 $, reproduces the observed dark matter relic density. This demonstrates the flexibility of the singlet χ NLOP scenario in achieving the correct dark matter yield through distinct but equally viable annihilation pathways.

The Dirac nature of neutrinos necessitates both left-handed ($ \nu_L $) and right-handed ($ \nu_R $) chiral components. The existence of $ \nu_R $ contributes additional relativistic degrees of freedom, thereby increasing the radiation energy density in the early universe. This contribution is commonly parameterized by the effective number of relativistic species, defined as

Within the SM, the effective number of relativistic species is precisely calculated to be $ N_{\rm eff}^{\rm SM} = 3.045 $ [74−76], incorporating effects from neutrino oscillations, non-thermal spectral distortions, and finite-temperature corrections. Any deviation from this value due to physics beyond the Standard Model (BSM) is typically expressed as $ \Delta N_{\rm eff} \equiv N_{\rm eff} - N_{\rm eff}^{\rm SM} $. The Planck2018 results provide the most accurate and stringent constraint to date, yielding $ \Delta N_{\rm eff} \leq 0.285 $ at the $ 2\sigma $ confidence level [77]. While the upcoming CMB-S4 experiment is projected to significantly improve sensitivity, with a forecasted reach of $ \Delta N_{\rm eff} \leq 0.06 $ at the same confidence level [78].

In the early hot universe, all particles except the FIMP $ N_1 $ were in thermal equilibrium with the SM plasma. As the universe expanded and cooled, the dark plasma---comprising the scalar singlet χ, fermion singlets $ N_{2,3} $, and right-handed neutrinos $ \nu_R $---gradually decoupled from the SM plasma. However, partial thermal equilibrium was maintained within the dark plasma due to the large Yukawa coupling $ y_\chi $. Interactions between the dark plasma and the SM plasma can occur through either $ N_{2,3} $ or χ. For $ N_{2,3} $, interactions with the SM are mediated by the coupling $ y_\Phi $, with the scattering amplitude scaling as $ |y_\Phi|^2 $. In contrast, χ couples to the SM either via t-channel Higgs exchange or through a contact interaction with Higgs, with the scattering amplitude proportional to $ \lambda_5 $. Due to stringent constraints from LFV, $ y_\Phi $ is typically too small to support significant interactions. As a result, the dominant portal between the dark plasma and the SM is mediated by χ. When χ eventually decouples from the SM plasma, the temperature of the right-handed neutrinos $ \nu_R $ begins to deviate from that of the SM thermal bath. The decoupling temperature $ T_{\text{dec}} $ is determined by the condition: $ \Gamma_{\text{el}}(T_{\text{dec}}) = H(T_{\text{dec}}) $, where $ \Gamma_{\text{el}} \equiv \sum n_{\text{SM}} \langle \sigma v \rangle_{\chi\,\text{SM} \to \chi\,\text{SM}} / (m_\chi / T) $ denotes the effective elastic scattering rate between χ and SM particles, accounting for the number of scatterings required to transfer energy of order T. In Fig. 5, we show the ratio $ \Gamma_{\rm el}/H $ as a function of temperature for different values of $ \lambda_5 $. The dashed line, corresponding to $ \Gamma_{\rm el}/H = 1 $, indicates the condition for decoupling. The temperature ratio $ T_{\nu_R}/T_{\rm SM} $ can be determined using entropy conservation:

Figure 5.  (color online) (Left panel) The ratio of $ \Gamma_{\rm el}/H $ as a function of temperature for different choice of $ \lambda_5 $. Here we have chosen $ m_\chi = 100\; \rm{GeV} $ for demonstration. (Right panel) The evolution of ratio $ T_{\nu_R}/T_\gamma $ for different dark plasma decoupling temperature.

where $ g_X^{*s} $ denotes the effective number of relativistic degrees of freedom, with respect to entropy density, in the X plasma (with "DP" referring to the dark plasma). The evaluation at $ T_{\rm dec} $ indicates that the ratio is set at the dark plasma decoupling temperature. The evolution of $ T_{\nu_R}/T_{\rm SM} $ is shown in Fig. 5 for different choices of the dark plasma decoupling temperature. The later the dark plasma decouples (which corresponds to a larger $ \lambda_5 $, and $ \nu_R $ is heated more efficiently by the SM plasma over a longer period), the higher temperature ratio we will arrive at. This finally results in a larger contribution to the effective number of relativistic species.

In addition to thermal contributions from $ \nu_R $, the delayed decay of the NLOP—particularly during the BBN or CMB epochs—can leave observable imprints on cosmological data. Several studies in the literature have examined the late decays of heavy relics into electron-positron pairs [79] or neutrinos [80], and constraints have been derived based on their effects on BBN, CMB anisotropies and spectral distortions. In our model, the NLOP candidates include $ N_2 $, $ \phi_{R,I} $, $ \phi^\pm $, and χ (The roles of $ \phi_R $ and $ \phi_I $ as the NLOP are entirely analogous, therefore the following discussion for $ \phi_R $ as the NLOP applies equally to $ \phi_I $). Among these, $ N_2 $ and $ \phi^\pm $ can decay into SM charged leptons. Exotic electromagnetic (EM) energy injection is tightly constrained by cosmology. CMB observations (e.g., Planck) limit such injection as it alters recombination and distorts temperature and polarization spectra, while COBE-FIRAS rules out spectral distortions from early injection. During BBN, excess energy disrupts light-element abundances, conflicting with observed D and $ ^4 $He. Long-lived $ N_2 $ and $ \phi^\pm $ with lifetimes between $ 3\times 10^4\,{\rm{s}} $ and recombination are strongly disfavored unless nearly degenerate with the dark matter (hence minimizing energy release) [79, 81−83]. In this work, we impose a conservative upper bound, requiring the lifetimes of $ N_2 $ and $ \phi^\pm $ to be less than $ 3\times 10^4\,{\rm{s}} $. For scenarios in which $ \phi_{R} $ or χ is the NLOP, their delayed decays can inject energy into the left- or right-handed neutrinos, thereby altering the neutrino temperatures. The evolution of $ T_{\nu_L} $, $ T_{\nu_R} $, and $ T_\gamma $ is governed by the following system of differential equations:

Here, ρ and p denote energy and pressure densities, respectively. $ \rho_{\rm NLOP} $ stands for the energy density of NLOP, it evolves as

$ {\delta \rho}/{\delta t} $ in Eq. 17 represents the energy transfer rate between the SM neutrinos and EM plasma due to weak interactions, as detailed in Refs. [84, 85]. The parameter $ \varepsilon^{\nu_L(\nu_R)}_{\rm NLOP} = (m_{\rm NLOP}^2 - m_{N_1}^2)/(2m_{\rm NLOP}^2) $ quantifies the fraction of the decaying NLOP's rest energy that is effectively deposited into $ \nu_L(\nu_R) $ ¹. The initial conditions for this system are set at a time when $ \nu_L $ and the photon bath were in thermal equilibrium ($ T_{\nu_L} = T_\gamma $), and the NLOP has not yet decayed. In practice, we begin the numerical calculation at $ T_\gamma = 10\; {\rm{MeV}} $, where $ \nu_L $ is still coupled to the EM plasma, and any decays occurring earlier would contribute negligibly to the entropy of $ \nu_L $ or $ \nu_R $. The starting temperature of $ \nu_R $ is fixed by Eq. 16 with setting $ T_{\rm SM} = 10\; \rm MeV $. Following Ref. [80], we define a dimensionless parameter $ f_{\rm NLOP} \equiv \Omega_{\rm NLOP} / \Omega_{\rm DM} $, where $ \Omega_{\rm NLOP} $ represents the hypothetical present-day relic abundance of the NLOP if it were stable. Consequently, the initial energy density of the NLOP in Eq. 18 can be expressed as a function of $ f_{\rm NLOP} $. In principle, a similar differential equation system should be solved for $ N_2 $ or $ \phi^\pm $ as the NLOP. However, we have verified that the BBN and CMB constraints (i.e. $ \tau < 3 \times 10^4\,{\rm{s}} $) ensure that any late-decay effects on $ \Delta N_{\rm eff} $ are negligible.

We show the evolution of $ T_{\nu_L}/T_\gamma $ (solid lines) and $ T_{\nu_R}/T_\gamma $ (dashed lines) in Fig. 6, where the first row corresponds to $ \phi_R $ as the NLOP and the second row to χ as the NLOP. The input parameters for Eqs. (17) and (18) are the NLOP energy fraction $ f_{{\rm{NLOP}}} $, its lifetime $ \tau_{{\rm{NLOP}}} $, and the temperature of the right-handed neutrino bath at $ T_\gamma = 10\; {\rm{MeV}} $, denoted $ T_{\nu_R,10} $. We illustrate how different choices of these parameters influence the resulting temperature ratios. In the first column, we fix $ f_{{\rm{NLOP}}} = 100 $ and $ T_{\nu_R,10} = 5\; {\rm{MeV}} $. The NLOP lifetimes are set to $ 10^5\; {\rm{s}} $, $ 10^6\; {\rm{s}} $, and $ 10^7\; {\rm{s}} $ for illustration. As the lifetime of $ \phi_R(\chi) $ increases, the delayed decay injects energy into $ \nu_L(\nu_R) $ at a later cosmic time, leading to a more pronounced heating effect on $ \nu_L(\nu_R) $. In the second column, we fix $ \tau_{{\rm{NLOP}}} = 10^7\; {\rm{s}} $ and $ T_{\nu_R,10} = 5\; {\rm{MeV}} $, and vary $ f_{{\rm{NLOP}}} = 10 $, $ 100 $, and $ 1000 $. A larger $ f_{\phi_R}(\chi) $ implies a higher relic density after freeze-out, resulting in greater energy injection into $ \nu_L(\nu_R) $ once the decays are complete. In the third column, we examine the impact of varying $ T_{\nu_R,10} $ on the temperature ratio evolution. When $ \phi_R $ is the NLOP, its late decay does not affect the thermal history of $ \nu_R $, and hence there is no heating of $ \nu_R $. When χ is the NLOP, variations in $ T_{\nu_R,10} $ do not influence $ T_{\nu_L}/T_\gamma $, since χ decays exclusively into $ \nu_R $. A higher $ T_{\nu_R,10} $ results in a larger $ T_{\nu_R}/T_\gamma $ at the completion of χ decay.

Figure 6.  (color online) The evolution of the temperature ratios $ T_{\nu_L}/T_\gamma $ (solid lines) and $ T_{\nu_R}/T_\gamma $ (dashed lines). The upper panel shows the case where $ \phi_R $ is the NLOP, and the lower panel shows the case where χ is the NLOP.

With the time-dependent temperatures $ T_{\nu_L} $, $ T_{\nu_R} $, and $ T_\gamma $ determined, $ \Delta N_{{\rm{eff}}} $ can be computed using Eq. (15). From Eq. (17), it follows that the energy injection from NLOP decay into $ \nu_L $ or $ \nu_R $ proceeds in an identical manner in both scenarios. As a result, the evolution of $ \Delta N_{{\rm{eff}}} $ is the same whether the NLOP is $ \phi_R $ or χ. We present the resulting $ \Delta N_{{\rm{eff}}} $ in Fig. 7, where the three subfigures show the effects of varying $ \tau_{{\rm{NLOP}}} $, $ f_{{\rm{NLOP}}} $, and $ T_{\nu_R,10} $, respectively. The two dashed lines stand for the current experimental limits from Planck [77] and future detecting capability from CMB-S4 [78]. A longer NLOP lifetime or a larger $ f_{{\rm{NLOP}}} $ leads to a more pronounced heating effect, resulting in a higher $ \Delta N_{{\rm{eff}}} $. A hotter $ \nu_R $ bath, originating from stronger and more prolonged interactions between the dark plasma and the SM plasma, also yields a larger $ \Delta N_{{\rm{eff}}} $. For completeness, we also show the evolution of $ \Delta N_{{\rm{eff}}} $ in the scenario where the contribution from the late decay of the NLOP is negligible, as depicted in the final subfigure of Fig. 7. This occurs when either the NLOP abundance is small or the NLOP lifetime is sufficiently short. In this limit, $ \Delta N_{{\rm{eff}}} $ is dominated solely by the thermal contribution from the decoupled $ \nu_R $ bath. As expected, a higher $ T_{\nu_R,10} $ results in a larger $ \Delta N_{{\rm{eff}}} $, since the $ \nu_R $s contribute more to the radiation density.

Figure 7.  (color online) Evolution of $ \Delta N_{{\rm{eff}}} $ for varying $ \tau_{{\rm{NLOP}}} $, $ f_{{\rm{NLOP}}} $, and $ T_{\nu_R,10} $. The final subfigure illustrates the scenario in which $ \Delta N_{{\rm{eff}}} $ is dominated by thermal contributions.

Finally, we perform a comprehensive scan over the full parameter space, imposing constraints from dark matter relic density, neutrino mass, LFV, and cosmological observations such as BBN, CMB, and $ \Delta N_{{\rm{eff}}} $. The scanning ranges for the model parameters are as follows:

Here $ m_{\text{NLOP}} $ is the mass of the NLOP, and $ m_{\text{others}} $ refers to the masses of heavier odd-sector states. The lower bound on $ m_{\phi^{\pm}} $ is set to 560 GeV, based on collider searches for long-lived charged particles [72]. The dark matter relic density is required to lie within the $ 3\sigma $ range of the Planck2018 measurement [77], i.e., $ \Omega_{N_1} h^2 \in [0.117, 0.123] $. The odd states masses are scanned up to energy scale within reach of current and near-future experiments. The couplings $ (y_{\Phi,\chi})_{\alpha 1} $ are bounded from above to ensure that dark matter remains out of thermal equilibrium. The couplings $ (y_{\Phi})_{\alpha 2,3} $ are constrained by LFV bounds to be less than $ 10^{-2} $, while $ (y_{\chi})_{\alpha 2,3} $ are chosen to reproduce observed neutrino masses. The mixing parameter $ \sin 2\theta $ is restricted to below $ 10^{-5} $ on naturalness grounds—since $ \theta = 0 $ (or equivalently $ \mu = 0 $) enhances the model's symmetry. The quartic coupling $ \lambda_5 $ is scanned down to $ 10^{-3} $, as the thermal contribution to $ \Delta N_{\rm eff} $ is no longer significant for smaller values.

As discussed in the previous subsection, we impose a conservative upper limit of $ \tau < 3 \times 10^4 \,\text{s} $ on the lifetime of the NLOP when it is either the charged scalar $ \phi^\pm $ or the neutral fermion $ N_2 $. This bound arises because EM energy injection from late decays can alter light element abundances during BBN and distort CMB anisotropies. However, when the NLOP is the neutral scalar $ \phi_R $, the constraints are significantly relaxed, as its dominant decays inject energy primarily into neutrinos rather than the EM plasma. In this case, we adopt the bounds from Ref. [80], which depend on the effective energy transfer parameter $ f_{\phi_R} \varepsilon_{\phi_R} $, and display them as the blue curve in Fig. 8. In the scan, parameter points satisfying the dark matter relic density condition are colored pink. Among these, points that additionally satisfy constraints from neutrino mass and LFV processes are highlighted in green. Points that further pass the BBN and CMB bounds are shown in red. One can observe that the BBN and CMB limits exclude long-lived $ \phi_R $ with lifetimes $ \gtrsim 10^4 $–$ 10^7\,\text{s} $, although the constraints weaken significantly when the energy transferred to the neutrino sector is small. When the NLOP is the scalar χ, which decays directly into right-handed neutrino, there are no direct BBN or CMB constraints, as the decay products are sterile and do not interact with the thermal bath.

Figure 8.  (color online) Scan results in the $ f_{\phi_R} \varepsilon_{\phi_R} $–$ \tau_{\phi_R} $ plane for the case where $ \phi_R $ is the NLOP. The blue curve shows the combined constraints from BBN and CMB as derived in Ref. [80].

In Fig. 9, we present the scanning results in the $ m_{{\rm{NLOP}}} $-$ m_{N_1} $ plane, considering different particles as the NLOP. Constraints from dark matter relic density, neutrino mass and LFV, BBN and CMB, and $ \Delta N_{{\rm{eff}}} $ are applied successively. The color coding follows the same progression as in Fig. 8. Finally, points that also satisfy the $ \Delta N_{{\rm{eff}}} $ constraint are marked in blue, indicating full consistency with all observational requirements. The dark matter relic density receives contributions from both the freeze-in and super-WIMP mechanisms, as discussed in detail in Sect. III. Constraints from LFV, particularly the stringent experimental limit on the branching ratio of $ \mu \to e\gamma $, exclude a significant portion of the parameter space. The impact of LFV constraints is also visible in the subfigure where χ is the NLOP, even though χ itself does not directly participate in LFV processes. This occurs because the same Yukawa couplings that govern LFV amplitudes also influence the dark matter production rate and thus must satisfy both relic density and flavor constraints simultaneously. The bounds from BBN and CMB have been discussed previously and are applied depending on the identity of the NLOP. In all cases—whether $ \phi_R $, χ, $ \phi^\pm $, or $ N_2 $ serves as the NLOP—we identify regions of parameter space that satisfy all observational constraints, including the latest limits on $ \Delta N_{{\rm{eff}}} $. This demonstrates the viability of the FIMP Dirac scotogenic Model as a unified framework for explaining both neutrino mass generation and the observed dark matter relic abundance. Future experiments will play a crucial role in testing this scenario. Long-lived particle searches at proposed facilities such as MATHUSLA [73] could probe the metastable NLOP states, while next-generation CMB observations from CMB-S4 [78] will improve sensitivity to $ \Delta N_{{\rm{eff}}} $ and energy injection during cosmic evolution. These experiments will either confirm the predictions of the model or further constrain its surviving parameter space.

Figure 9.  (color online) Results of parameter scan in the $ m_{\rm NLOP}-m_{N_1} $ plane under combined constraints from DM relic density, LFV, neutrino mass, CMB, BBN, and $ \Delta N_{\rm eff} $.

In this work, we have studied the FIMP Dirac Scotogenic Model, a well-motivated extension of the SM that simultaneously addresses the origin of neutrino masses and the nature of dark matter. The lightest neutral fermion $ N_1 $ is a feebly interacting massive particle DM candidate, produced out of equilibrium. Neutrino masses are generated radiatively at one loop, resulting in a rank-2 mass matrix due to highly suppressed couplings to $ N_1 $, predicting one nearly massless neutrino. Stringent LFV constraints, require the Yukawa couplings $ (y_\Phi)_{\alpha 2,3} $ to be small, which in turn necessitates larger $ (y_\chi)_{\alpha 2,3} $ to fit the neutrino mass scale. This coupling structure crucially impacts the phenomenology, governing the freeze-out dynamics of the NLOP and the thermal history of the right-handed neutrino bath.

We have systematically analyzed two distinct production mechanisms for $ N_1 $: direct freeze-in via decays of heavier $ Z_2 $-odd scalars and fermions, and the "super-WIMP" mechanism, where $ N_1 $ is produced from the late decay of a thermally frozen-out NLOP. The latter scenario leads to rich phenomenology, depending on whether $ N_2 $, $ \phi_{R,I} $, $ \phi^\pm $, or χ serves as the NLOP. We demonstrate that the correct relic abundance can be achieved across various benchmarks, with coannihilation effects and enhanced annihilation channels playing crucial roles in regulating the NLOP's freeze-out density, offering flexible pathways to match the observed DM density.

Another central focus of this work is the cosmological impact on the effective number of relativistic species. Our analysis reveals two distinct contributions to this quantity. The first is a thermal contribution, arising from the primordial right-handed neutrino bath. This component is determined by the decoupling temperature of the $ \nu_R $ bath, which is set by the strength of interactions between the dark sector and the SM plasma. The second contribution is non-thermal, originating from the late decays of the NLOP. The magnitude of this non-thermal contribution is highly sensitive to the NLOP's properties: both a longer lifetime $ \tau_{{\rm{NLOP}}} $ and a larger energy fraction $ f_{{\rm{NLOP}}} $ lead to a more pronounced heating effect and thus a higher $ \Delta N_{{\rm{eff}}} $.

We have performed a comprehensive parameter scan, sequentially applying constraints from DM relic density, neutrino mass and LFV, BBN/CMB, and $ \Delta N_{{\rm{eff}}} $. Viable parameter space exists for all NLOP candidates that satisfy all current observational bounds. Future experiments setup for long-lived particle searches and precision $ \Delta N_{{\rm{eff}}} $ measurements will provide powerful probes to either discover this framework or further constrain its remaining parameter space.

Finally, It is instructive to contrast the two possible dark matter production mechanisms in our work: thermal freeze-out (WIMP) [63] versus freeze-in (FIMP). The distinction is most pronounced in the size of the Yukawa couplings $ (y_{\Phi,\chi})_{\alpha 1} $ that govern the interaction of the lightest odd particle $ N_1 $ with the Standard Model sector. In the WIMP scenario, these couplings must be relatively large to ensure efficient annihilation and avoid overproduction of dark matter. In the FIMP regime, by contrast, the couplings are extremely suppressed, $ (y_{\Phi,\chi})_{\alpha 1} \lesssim 10^{-7} $, keeping $ N_1 $ out of thermal equilibrium throughout cosmic history. The relic density is then generated via the slow decay of heavier states, and this suppression naturally leads to a rank-2 neutrino mass matrix, predicting one nearly massless neutrino—a feature not generic in WIMP realizations. Phenomenologically, the two scenarios are also sharply distinguished by experimental constraints. WIMP parameter space is heavily constrained: collider searches exclude light dark matter masses, while direct and indirect detection experiments probe the TeV-scale region. In contrast, FIMPs evade these conventional probes due to their feeble couplings. Their collider signatures are limited to displaced decays of the NLOP (e.g., displaced vertices), and current long-lived particle searches provide only weak limits. Direct and indirect DM detection signals are negligible. Instead, the primary observational window for FIMPs lies in cosmology—particularly through their contribution to $ \Delta N_{{\rm{eff}}} $ from late decays of long-lived particles—which offers a powerful and complementary test of the freeze-in paradigm.

We thank Zhi-Long Han, Ang Liu and Xiao-Dong Ma for helpful discussions.