When FIMPs Decay into Neutrinos: The $N_\mathrm{eff}$ Story

The existence of feebly interacting massive particles (FIMPs) could have significant implications on the effective number of relativistic species $N_\mathrm{eff}$ in the early Universe. In this work, we investigate in detail how short-lived FIMPs that can decay into neutrinos affect $N_\mathrm{eff}$ and highlight the relevant effects that govern its evolution. We show that even if unstable FIMPs inject most of their energy into neutrinos, they may still decrease $N_{\mathrm{eff}}$, and identify neutrino spectral distortions as the driving power behind this effect. As a case study, we consider Heavy Neutral Leptons (HNLs) and indicate which regions of their parameter space increase or decrease $N_{\mathrm{eff}}$. Moreover, we derive bounds on the HNL lifetime from the Cosmic Microwave Background and comment on the possible role that HNLs could play in alleviating the Hubble tension.


I. INTRODUCTION
The Standard Model (SM) of particle physics and cosmology has proven to be very successful in describing our observable Universe.Nevertheless, in its current form, it does not provide a physical origin for a number of observed phenomena.The baryon asymmetry of the Universe, the existence of neutrino masses and the evidence for the dark matter each establish that the SM picture is not complete and call for the addition of new physics, usually in the form of a new particle [1][2][3].While the landscape of viable candidates is almost limitless, extensions of the SM often involve feebly interacting massive particles (FIMPs) that possess small couplings to the SM sector, see e.g.[4][5][6][7] for reviews.Their inclusion in the primordial plasma comes with potential consequences on fundamental probes in the early Universe, such as primordial nucleosynthesis and the Cosmic Microwave Background (CMB) [8][9][10].Indeed, just their presence in the system would already contribute to the dynamics of the expanding Universe, let alone further implications that could follow from their decay.
A key parameter in this topic is the effective number of relativistic species N eff , given by: 11 4 4/3 where ρ rad and ρ γ are the total radiation and photon energy densities respectively.We define the change in this quantity as ∆N eff = N eff − N SM eff , where within the Stan- A boyarsky@lorentz.leidenuniv.nlM ovchynnikov@lorentz.leidenuniv.nlN nashwan.sabti@kcl.ac.ukV vsevolod.syvolap@nbi.ku.dk dard Model N SM eff 3.044 [11][12][13][14][15]. Any deviation from the SM value is regulated by weak interactions between neutrinos and electromagnetic (EM) particles, which are efficient enough at temperatures T 1 MeV to keep these species in equilibrium with each other.At lower temperatures, the interactions gradually go out of equilibrium and the energy exchange between the two sectors will stop.Decaying FIMPs can affect this delicate process in different ways, depending on whether they inject most of their energy into EM particles or neutrinos.
The impact of FIMPs predominantly decaying into EM particles has been extensively studied in the literature, see e.g.[16][17][18][19][20][21].Such particles heat up the EM plasma and consequently decrease N eff , independently of whether the decay happens during or after neutrino decoupling.On the other hand, for FIMPs that mostly decay into neutrinos, we naively expect that N eff would increase.This is indeed true for lifetimes τ FIMP t dec ν ∼ 0.1 − 1 s, where t dec ν is the time of neutrino decoupling, see e.g.[22].However, the situation for lifetimes τ FIMP ∼ t dec ν is different: the neutrinos are still in partial equilibrium and try to equilibrate with the injected neutrinos.This scenario has been considered before in [23][24][25] that arrived at different conclusions about the impact on N eff .The work [23] studied a reheating scenario in which all the SM particles are absent before FIMPs start decaying.In such a framework, all neutrinos have high energies, which means that they mainly thermalise via neutrino-EM interactions and N eff naturally decreases.References [24,25] considered Heavy Neutral Leptons (HNLs) with masses m N < m π and lifetimes τ N 1 s.Such HNLs are in thermal equilibrium in the early Universe, but decouple as the Universe expands and eventually decay mainly into high-energy neutrinos at MeV temperatures.These two works drew different conclusions about N eff : [24] reported ∆N eff > 0 for the whole studied mass range, whereas [25] presented in their Fig. 3 that ∆N eff < 0 for masses 60 MeV m N < m π and lifetimes arXiv:2103.09831v2[hep-ph] 20 Jul 2021 τ N 1 s.The sign of ∆N eff is not emphasised in these two papers; [25] did not comment on the contradiction with [24] on this issue and no physical discussion of this phenomenon was provided 1 .
In this paper, we aim to clarify the behaviour of N eff in the presence of FIMPs that decay mainly into neutrinos and have lifetimes τ FIMP ∼ t dec ν .Here we will assume that a thermal bath of SM particles is already present in the primordial Universe.We will first construct a simple model in Sec.II that provides us with a qualitative understanding of how such particles impact N eff , the findings of which we then confirm by using the Boltzmann code pyBBN [21].Our analysis shows that short-lived FIMPs that inject most of their energy into neutrinos may decrease N eff .This is because during the equilibration process, the injected high-energy neutrinos redistribute their energy among the neutrino and EM plasma.If the energy of the injected neutrinos is sufficiently large, the energy transfer to the EM sector occurs faster than the equilibration with the neutrino sector.This means that the EM plasma heats up more than the neutrino plasma, which eventually leads to ∆N eff < 0. We will find that this mechanism is especially relevant for FIMPs with masses larger than a few tens of MeV.We will then apply these general considerations to the well-motivated case of Heavy Neutral Leptons.This is done in Sec.III, where we will also briefly discuss their influence on the physics at the CMB epoch, derive bounds on their mass and lifetime, and comment on the implications their presence may have for the Hubble tension.Our main results are summarised in Figs. 2 and 5. Finally, we will present our conclusions in Sec.IV.Complementary details and simulation results are included in the appendices A−C.

II. HOW FIMPS IMPACT Neff
The goal of this section is to provide a physical understanding of the processes that govern the behaviour of N eff in the presence of decaying FIMPs.In particular, we focus on FIMPs with masses 1 MeV that decay when neutrinos are still in (partial) equilibrium.Such FIMPs can decay into high-energy neutrinos that then participate in interactions with thermal neutrinos and electrons/positrons.We will find that even if most of the FIMP energy is injected into neutrinos, these interactions may still cause a decrease in N eff .This feature ap-1 A more recent work [26] considered long-lived (i.e., decaying after e + e − annihilation) HNLs that could decay both into EM particles and neutrinos.In this case, N eff could both increase and decrease, as at such late times the injected energy densities from HNL decays can dominate over the SM densities of both the EM and neutrino sectors.Another recent paper [27], which appeared after our work was submitted, claims that ∆N eff ≥ 0 for all cases in which FIMPs decay mostly into neutrinos.We comment on it in Appendix D.
pears since the injected high-energy neutrinos get quickly converted into electrons/positrons and drag thermal neutrinos residing in the plasma along with them.During this process, neutrino-neutrino interactions lead to the presence of residual non-thermal distortions in the distribution functions of neutrinos (neutrino spectral distortions) that keep the balance of ν ↔ EM interactions shifted to the right till long after the injection (i.e., more energy is transferred from the neutrino plasma to the electromagnetic plasma than vice versa).The energy transfer from neutrinos to EM particles accumulated over time can then be sizeable enough, such that ∆N eff becomes negative.This effect diminishes with larger FIMP lifetime, as neutrino-EM interactions go out of equilibrium and neutrinos can no longer be converted into electrons/positrons.Therefore, FIMPs that decay into neutrinos after neutrino decoupling will lead to ∆N eff > 0.
In what follows, we will consider FIMPs that can decay into both neutrinos and EM particles, and construct a simple model that provides a semi-analytic description of the aforementioned effect.At the end of this section and in Appendix A, we will also highlight and further elaborate on the central role of neutrino spectral distortions in the dynamics of N eff .
A FIMP of mass m FIMP T can decay into neutrinos with energies that are larger than those found in the primordial plasma.Such non-equilibrium neutrinos manifest themselves as spectral distortions at high energies and will subsequently interact with the thermal neutrinos and electrons/positrons in the plasma.We assume that the amount of injected non-equilibrium neutrinos is only a small fraction of the thermal neutrinos in the plasma.The evolution of the injected neutrinos is then mainly governed by the following reactions: ν non-eq + ν therm → e + + e − where 'non-eq' and 'therm' refer to neutrinos with nonequilibrium and thermal energies respectively.Through these reactions, non-equilibrium neutrinos thermalise and quickly redistribute their energy among the neutrino and EM plasma.The energy loss rate of these nonequilibrium neutrinos is higher than the interaction rates of thermal particles [28]: where E inj ν is the average energy of the injected nonequilibrium neutrinos.Note that reactions between thermal particles also exchange energy between the neutrino and EM sectors, but this energy exchange is subdominant as far as Eq. ( 5) holds.
The amount of energy that ends up in the EM plasma has three contributions: 1) the direct decay of FIMPs into EM particles, 2) the energy transfer of non-equilibrium neutrinos to EM particles during thermalisation and 3) the energy transfer from thermal neutrinos to EM particles as a consequence of them being dragged by non-equilibrium neutrinos during thermalisation (reactions (2) and ( 3)).The first process injects a fraction ξ EM of the total FIMP energy into the EM plasma, while the latter two increase this fraction to: where ξ ν = 1 − ξ EM is the energy fraction that FIMPs directly inject into the neutrino sector and = non-eq + thermal is the effective fraction of ξ ν that went to the EM plasma during the thermalisation.The latter quantity can be split in a contribution from non-equilibrium neutrinos ( non-eq = E non-eq→EM ν /E inj ν ) and an effective contribution from thermal neutrinos Now, based on Eq. ( 6), if > 0.5, then ξ EM,eff > 0.5.This means that more than half of the FIMP energy eventually ends up in the EM plasma (i.e., EM plasma heats up more than the neutrino plasma), which results in ∆N eff < 0 independently of the value of ξ EM2 .This simplified energy redistribution picture only holds if the non-equilibrium neutrino energy is much larger than the average energy of thermal neutrinos.Once these two energies become similar in magnitude, backreactions cannot be neglected anymore and the evolution can only be accurately described with a system of Boltzmann equations.
We can make a simple estimate of as a function of the injected neutrino energy E inj ν and temperature T .We start with describing the thermalisation process of a single injected neutrino, which causes a cascade of nonequilibrium neutrinos.Such a cascade can result after the injected neutrino participates in the processes (2)−(4).We assume that in the processes (2) and (4) each nonequilibrium neutrino in the final state carries half of the energy of the non-equilibrium neutrino in the initial state.Thus, roughly speaking, the thermalisation occurs during N therm log 2 (E inj ν /3.15T ) interactions.In addition, the process (2) doubles the number of non-equilibrium neutrinos, while (3) makes neutrinos disappear and (4) leaves the number unchanged.Therefore, after the k-th step in the cascade, the average number of non-equilibrium neutrinos is given by: with N (0) ν = 1, and the total non-equilibrium energy is: where P νν→νν , P νν→ee , and P νe→νe are the average probabilities of the processes (2)−(4), respectively, and their sum equals unity.We define these probabilities as P i = Γ i /Γ tot ν , where Γ i is the interaction rate of each process and Γ tot ν is the total neutrino interaction rate.The relevant reactions and their corresponding matrix elements are summarised in appendix D of [21].Assuming a Fermi-Dirac distribution for neutrinos and averaging over neutrino flavours, we find: Finally, the value of non-eq that accounts for the energy transfer from non-equilibrium neutrinos to the EM plasma is given by: In addition to the transferred non-equilibrium energy, the non-equilibrium neutrinos catalyse the energy transfer from thermal neutrinos to the EM plasma via the processes ( 2) and ( 3).In other words, during the thermalisation process non-equilibrium neutrinos drag thermal neutrinos along with them, which leads to part of the energy stored in the thermal neutrino sector to end up in the EM sector.We assume that each reaction (2) transfers an energy amount of 3.15T from the thermal neutrino sector to non-equilibrium neutrinos, which then via (3) ends up in the EM plasma.Moreover, each reaction (3) contributes to another energy transfer of 3.15T from thermal neutrinos to the EM plasma.The effective contribution coming from this transfer is therefore: where the first term in the round brackets is the contribution from the process (3) and the terms in the square brackets are the contribution from the process (2).Note that the factor of 2 in the second sum accounts for the doubling of non-equilibrium neutrinos in the process (2).We find that thermal is at least 5 times smaller than non-eq , which makes this a sub-dominant effect.
As the Universe expands and the temperature decreases, weak reaction rates start to compete with the Hubble ν is injected at a temperature Tinj.At high temperatures of order of Tinj E inj ν , the injected neutrinos are thermal-like, and hence is small.Once the temperature decreases, we enter the regime E inj ν 3.15Tinj and neutrinos transfer a significant amount of their energy to the EM plasma while thermalising.With further decrease of Tinj, weak reactions go out of equilibrium and the energy transfer becomes less and less efficient, which results in a quick drop-off of .rate H.The energy transfer from neutrinos to the EM plasma therefore becomes less and less efficient, and tends to zero.In order to incorporate this effect, we multiply the probabilities in (9) with a factor min[Γ i /H, 1], where Γ i = Γ i (E inj ν /2 k ) is the interaction rate of any of the processes (2)−(4).The resulting energy fraction of neutrinos that is transferred to the EM plasma = non-eq + thermal is shown in Fig. 1 for a number of injected neutrino energies E inj ν .We find that can exceed 0.5 for E inj ν 60 MeV.This means that when FIMPs decay into neutrinos with such energies at temperatures of a few MeV, the majority of the injected neutrino energy will end up in the EM plasma during the thermalisation.This then leads to a decrease of N eff , independently of how much energy the FIMPs inject into the EM sector.Now that we are able to estimate , we can compute the correction to N eff for some benchmark FIMP scenario.It is worth noting here again that only depends on the energy of the injected neutrino and the temperature at which the injection happens.This means that is an independent quantity of the FIMP model considered, in contrast to ξ EM and ξ ν , which do depend on the choice of the model.As an illustrative example, we assume that ξ EM = 0, i.e., the FIMP injects all of its energy into neutrinos (ξ ν = 1).Given that in our simple model neu-trinos thermalise very quickly, we assume that they have a thermal-like distribution with a temperature T ν and follow the approach in [12,29] to obtain the time evolution of T ν and T EM in the presence of decaying FIMPs (see Appendix B, where we provide the relevant equations).In this benchmark example, we consider a generic FIMP of mass 500 MeV that can decay only into three neutrinos and show ∆N eff as a function of its lifetime in Fig. 2. In order to compare the accuracy of our simple model, we also include in this figure the evolution of ∆N eff as obtained from the publicly available Boltzmann code pyBBN3 [21].The grey band in this figure indicates the current sensitivity of N eff by Planck, which at 2σ reads4 N CMB eff = 2.89 ± 0.62 [30,31].We see that N eff can significantly decrease as a result of the thermalisation of the injected neutrinos.This decrease of N eff would only be further amplified if the FIMPs were also to inject some of their energy into the EM plasma.Importantly, we find that when using the Boltzmann equation N eff decreases more than predicted by our semi-analytic model.This is because our model assumes that the remaining fraction 1 − of the injected neutrinos is perfectly thermal.In fact, this is usually not the case and the remaining non-equilibrium neutrinos manifest themselves as residual spectral distortions in the distribution function of neutrinos that further lead to a transfer of energy from the neutrino sector to the EM sector.We see that in some cases the inclusion of this effect can make the difference between being excluded by current data or not.We elaborate more on the effect of spectral distortions in Appendix A. In short, the semi-analytic model is useful in providing a qualitative understanding of the behaviour of N eff in the presence of decaying FIMPs.On the other hand, if the aim is to obtain accurate predictions for N eff (relevant for setting bounds and forecasting), it is crucial to use the Boltzmann equation to track the evolution of the neutrino distribution functions.As such, we will use pyBBN in the remainder of this paper to simulate the impact of FIMPs on N eff .
As a final point, we can make a rough model-independent estimate for which neutrino energies the decrease of N eff happens.In the particular FIMP scenario considered here, we find that this effect occurs for masses higher than ∼70 MeV (see Fig. 2).Given that in this case the neutrinos are created via 3-body decays, this would correspond to an average injected neutrino energy of roughly E inj ν ∼ m FIMP /3 ∼ 25 MeV.Therefore, as long as a FIMP injects most of its energy into neutrinos around neutrino decoupling, N eff could decrease if neutrinos with energies of at least E inj ν ∼ 25 MeV are produced.The initial FIMP abundance is assumed to be nχ/s = 0.01 at T = 1 GeV, where s is the total entropy density of a universe consisting of photons, neutrinos and electrons/positrons.The solid lines are the result of our semi-analytic model, while the dotted lines are obtained with the Boltzmann code pyBBN.The grey band is the current sensitivity by Planck (see text for details).The golden curves roughly indicate the lowest FIMP mass for which N eff can decrease due to the thermalization of the injected neutrinos.The stronger decrease of the blue, dotted curve as compared to the solid curve highlights the significance of residual neutrino spectral distortions in the evolution of N eff (see Appendix A for more details).

III. CASE STUDY: HEAVY NEUTRAL LEPTONS
A specific class of FIMPs that has seen an increasing interest over the last few years involves Heavy Neutral Leptons (HNLs).HNLs can be regarded as the right-handed companions to the SM neutrinos, which in a natural way give rise to neutrino masses [32,33].Besides this, they could play the role of the dark matter [34,35] and provide a mechanism for the generation of the baryon asymmetry in the Universe [36][37][38].We consider HNLs that couple to the SM through the neutrino portal: where α = {e, µ, τ }, β = 1, 2, 3, . .., F αβ are dimensionless Yukawa couplings, L α is the SM lepton doublet, H ≡ iσ 2 H * is the Higgs doublet in the conjugated representation and N β are the HNLs.This portal induces an effective interaction with SM particles that is similar to that of active neutrinos, but with an additional suppression in the form of small mixing angles [39].The mixing angle also parametrises with which neutrino flavour generation HNLs mix.An extension of the SM with three HNLs is also known as the neutrino Minimal Standard Model (νMSM) [40][41][42].This framework includes one light HNL, that plays the role of the dark matter, and two heavier ones, which are able to account for neutrino masses and the baryon asymmetry.
In what follows, we will consider two quasi-degenerate HNLs [43,44], with masses up to m N ∼ 1 GeV and lifetimes down to 0.01 s.Higher masses are not considered, as we do not include decays into multi-meson final states that become relevant in this mass range [39] and for which currently there is no adequate description of the corresponding decay widths 5 .Moreover, HNLs with lifetimes shorter than 0.01 s change N eff well below our accuracy (which is at the 1% level).We make use of pyBBN [21] to simulate their impact on the cosmological history, in particular on N eff and the primordial helium abundance Y P .We examine the region of parameter space in which HNLs inject most of their energy into neutrinos, but where ∆N eff is negative, illustrating the effect described in the previous section.Finally, we derive bounds from the CMB and comment on the possible role of HNLs in alleviating the Hubble tension.In Appendix C, we include fitting functions for N eff as obtained from the simulations.

A. Impact on Cosmological History
HNLs alter the cosmological history through their contribution to the total energy density of the Universe and their decay into SM particles.HNLs that decay well before the decoupling of active neutrinos will leave no traceable impact.On the other hand, if HNLs live long enough, they could alter several physical quantities, such as N eff and the primordial abundances of light elements [22,24,25,[45][46][47].Indeed, strong limits have been set on their mass and lifetime by considering their impact on Big Bang Nucleosynthesis and the Cosmic Microwave Background, see e.g.[21,26,48] for recent works on this subject.HNLs inject (eventually) all of their energy either into the neutrino or electromagnetic plasma.The fraction of the HNL energy that is injected into each of these two sectors is mass-dependent and shows a significant shift to the EM plasma once HNLs can decay into neutral pions (∼135 MeV), see Fig. 3.This plot shows that HNLs below the pion mass decay mainly into neutrinos and, therefore, one would naively expect that in this mass range N eff increases.However, we find that HNLs are able to decrease N eff for masses already above ∼70 MeV, while for smaller masses an increase of N eff The fraction of HNL mass that is injected into the neutrino plasma.Contributions to this fraction from unstable HNL decay products (mesons and muons) are included and we assume that the kinetic energy of all created charged particles goes into the EM plasma.For mN 135 MeV, HNLs can decay into neutral pions, which in their turn decay into two photons.This causes the sudden decrease of ξν around that mass.At higher masses, ξν keeps increasing in the case of τ -mixing, which is due to the absence of HNL decays into charged mesons (such decays are possible in the other two mixing cases).
is observed.The origin of this sign change in ∆N eff at m N 70 MeV (rather than m N > m π as one would guess from Fig. 3) lies in the energy transfer from the neutrino plasma to the electromagnetic plasma that is induced by the injected non-equilibrium neutrinos, as discussed earlier in Sec.II.We run pyBBN simulations to examine in which region of parameter space this sign change happens 6 .We show ∆N eff as a function of the HNL lifetime in Fig. 4 for a number of benchmark masses.The grey band in this figure indicates the current sensitivity by Planck.Included in this figure is an HNL of mass 110 MeV, which decreases N eff for lifetimes below τ N 0.6 s and increases N eff for longer lifetimes.Such a lifetime (τ N ∼ 0.6 s) roughly corresponds to the time of neutrino decoupling, beyond which thermalisation between the neutrino and EM plasma is not efficient anymore and the injected neutrinos remain in the neutrino sector.This exemplifies the ability of HNLs below the pion mass to diminish N eff , even when neutrinos are on the verge of being completely decoupled.With the cur- rent sensitivity of Planck, however, this initial decrease of ∆N eff for this mass falls within the error range and is thus not observable.Nevertheless, a number of upcoming and proposed CMB missions, such as the Simons Observatory [49] and CMB-S4 [50], could provide a determination of N eff around the percent-level and probe this effect.
We depict the region of HNL parameter space where ∆N eff changes sign in the top panel of Fig. 5.This is shown for the case of pure mixing with tau neutrinos only, as the parameter space where HNLs mix purely with electron and muon neutrinos is excluded in the lower mass range (where ∆N eff can be positive) by BBN, the CMB and experimental searches [21,48,52].In these latter two cases, ∆N eff can only be negative in the unconstrained parameter space.This top panel shows that there is a large region of HNL parameter space, where these particles inject most of their energy into neutrinos and still decrease N eff .The behaviour of negative ∆N eff continues for short-lived HNLs with masses m N > 1 GeV, since the neutrino energy increases with the HNL mass.On the other hand, for HNLs with lifetimes τ N 1 s, it depends on how much energy they inject into the neutrino plasma.Indeed, such HNLs decay long after neutrino decoupling, when non-equilibrium effects are not important anymore and the injected neutrinos remain in the neutrino sector.Mixing with tau neutrinos only is considered here.Top panel : Regions of the HNL parameter space that predict an increase (blue) or decrease (red) of N eff with respect to the SM value.The horizontal lines at the bottom of the plot indicate the mass ranges where HNLs inject most of their energy into neutrinos (ξν > 0.5) or the EM plasma (ξν < 0.5).In the former case, HNLs can still decrease N eff as a result of the efficient transfer of energy from neutrinos to EM particles.Bottom panel : Regions of the HNL parameter space that are excluded by BBN abundance measurements (green) and CMB observations (yellow).The ∆N eff = {0, ±0.4} contours give an indication of by how much HNLs can change N eff at the most in the unconstrained region.The BBN bound is from [21] and uses a central value for the primordial helium abundance of YP = 0.245 [51] with an error of 4.35% (see [48] for a discussion on how this error is obtained).For masses higher than the eta-meson mass (∼550 MeV), the meson effect from [48] is included in the analysis.The CMB constraint is obtained using the approach as detailed in [21] (see Sec. III B for more details on the CMB bound).This panel also shows that there is only a relatively small unconstrained region of parameter space left that can increase N eff and where HNLs could play a role in alleviating the Hubble tension.
This means that the sign of ∆N eff is simply determined by the value of ξ ν .As a result, for masses where ξ ν > 0.5 (see Fig. 3) this would mean that eventually ∆N eff > 0 and vice versa (see Fig. 6 for an illustration).

B. Bounds from the CMB
The CMB anisotropies are mainly sensitive to N eff through its impact on the damping tail [51,[53][54][55].For example, a larger number of relativistic degrees of freedom causes a stronger suppression of the power spectrum at high multipoles, as temperature anisotropies below the scale of the photon diffusion length are more damped by the increased expansion rate.This effect is, however, degenerate if the primordial helium abundance Y P is also considered as a free parameter [55].Y P is related to the number density of free electrons 7 , n e ∝ (1 − Y P ), which in its turn enters in the CMB damping scale.A larger Y P leads to a lower electron density, a larger electronphoton interaction rate, a larger photon diffusion length and thus a stronger damping.
We extend the CMB constraint on HNLs for masses up to 1 GeV using the same approach as detailed in [21,56] and show the result in the bottom panel of Fig. 5. Also included in this panel are the contours where ∆N eff = ±0.4,which give an indication of by how much HNLs can change N eff at the most, given the current constraints imposed by BBN and the CMB.The CMB bound is only stronger than the BBN bound in the lower mass range, as this is where N eff strongly increases.HNLs with short lifetimes and masses around O(10) MeV decouple while being non-relativistic and thus have a suppressed number density.They can therefore survive beyond the decoupling of SM weak reactions, without significantly affecting the primordial abundances.However, since the HNL energy density here falls off as (scale factor) −3 , the HNLs could eventually dominate the total energy density of the Universe.As can be seen in Fig. 3, HNLs in this lower mass range inject most of their energy into neutrinos, which remains in the neutrino sector after neutrino decoupling.The result is then a significant increase in N eff , which can be constrained with the CMB.On the other hand, for masses higher than ∼70 MeV, N eff starts decreasing.This decrease is relatively small in magnitude, especially in the region that is not constrained by BBN, where N eff − N CMB eff 0.4.In addition, the error in the determination of Y P by the CMB is larger than the one by BBN [30,56].These two properties make the CMB a weaker probe of HNLs in the higher mass range.As mentioned before, future CMB experiments could improve upon this result.Semi-analytic estimate of ∆N eff as a function of HNL mass and lifetime in the case of pure tau mixing.This plot is obtained using the method described in Sec.II (and is therefore only accurate up to a factor 3 − 4 for short lifetimes, when neutrinos are still in partial equilibrium).Nevertheless, it allows for a qualitative understanding of the behaviour of ∆N eff at lifetimes larger than considered in the main analysis (Fig. 5).Importantly, for lifetimes well beyond the time of neutrino decoupling (O(1) s), non-equilibrium effects are absent and the sign of ∆N eff is thus completely determined by the fraction of HNL energy ξν that is injected into the neutrino plasma, see Fig. 3.We see that HNLs with low masses and long lifetimes can still considerably affect N eff , while in the higher mass range ∆N eff tends to 0. This is because lowmass HNLs are more abundantly produced in this region of parameter space [26], where their mixing angles are relatively large.

C. Implications for the Hubble Parameter
An increase or decrease of N eff subsequently also changes the Hubble parameter.As such, HNLs could play a role in alleviating the longstanding tension between local determinations of the current day Hubble rate H 0 and the one as inferred from the CMB 8 [58,59].The usual approach involves increasing N eff , while keeping the angular scale of the sound horizon θ s = r s /D A fixed, see e.g.[60][61][62].Here, r s is the comoving sound horizon and D A is the comoving angular diameter distance to the surface of last scattering.Both of these quantities depend on the Hubble parameter: 8 This question has been considered before in [57].Importantly, this study used the results of [24], where the assumption was made that any change in the primordial helium abundance is due to ∆N eff .In contrast, here we find that neutrino spectral distortions are the driving power behind ∆N eff for short-lived HNLs.As a consequence, the results presented in our work and in [57] are rather different.
where z * is the redshift of the last-scattering surface and c s (z) is the speed of sound of the baryon-photon fluid in the early Universe.The Hubble rate in Eq. ( 13) depends mainly on the radiation (photons and neutrinos) and matter energy densities, while the one in Eq. ( 14) is the late-time Hubble rate and depends mostly on the dark energy and matter energy densities.This means that increasing N eff only results in a larger early-time Hubble rate and a r s .In order to keep θ s fixed, the comoving angular diameter distance must satisfy D A = r s /θ s , which then also decreases if r s decreases.Looking at Eq. ( 14), such a decrease can be accomplished by increasing the dark energy density ω Λ , or equivalently, H 0 (as Ω Λ = 1 − Ω m ).Since local measurements find a larger value of H 0 than the one inferred from the CMB within the Standard Model, this approach provides a way to reduce the Hubble tension.
This method, however, does not take into account the increased Silk damping induced by a larger N eff [51,[53][54][55].Therefore, a price must be paid when alleviating the Hubble tension in this way: An increase of N eff leads to a larger disagreement with the CMB itself.Given our CMB constraint in Fig. 5, we see that HNLs can increase N eff by at most ∆N eff ≈ 0.4.This gives us an indication of the extent to which unconstrained HNLs could increase H 0 and ameliorate the Hubble tension.We estimate the corresponding H 0 by running Monte Python [63,64] with the Planck 2018 baseline TTTEEE+lowE analysis.

IV. CONCLUSIONS
In this work, we have studied how heavy, unstable FIMPs that can decay into neutrinos impact the number of relativistic species N eff in the early Universe.A particularly interesting effect that could occur with these parti-cles, is when they inject most of their energy into neutrinos but still decrease N eff .This could happen if FIMPs decay when neutrinos are still in (partial) equilibrium (τ FIMP ∼ O(0.1) s) and is a direct consequence of the thermalisation process of the injected high-energy neutrinos (see Sec. II for a semi-analytical treatment of this effect).Here we identify neutrino spectral distortions as the driving power behind this effect, since they lead to an efficient transfer of energy from the neutrino plasma to the electromagnetic plasma (see Figs. 2 and 7).Some of the injected neutrino energy gets quickly transferred to the EM plasma, while the remaining will stay as residual spectral distortions in the neutrino distribution functions.These spectral distortions keep the energy transfer balance of ν ↔ EM reactions shifted to the right till long after FIMP decay.In order to accurately account for this effect, it is therefore important to solve the Boltzmann equation and track the evolution of the neutrino distribution functions.Using a thermal-like distribution for neutrinos as an approximation can lead to incorrect results, e.g., that N eff can never decrease when FIMPs inject most of their energy into neutrinos.
From our simulations, done with the publicly available Boltzmann code pyBBN [21], we find that this mechanism is especially relevant for FIMPs that can decay into neutrinos with average energies E inj ν 25 MeV.In case such neutrinos are created via 2-or 3-body decays, this roughly corresponds to FIMP masses m 2-body FIMP 50 MeV and m 3-body FIMP 70 MeV respectively.This is in agreement with the results presented in [23].As such, this effect may be relevant for many classes of FIMPs11 , including Higgs-like dark scalars [66], dark photons [67], neutralinos in supersymmetric models with broken R-parity [68], vector portals coupled to anomaly-free currents [69] and short-lived neutrinophilic scalars [70].
As a case study, we have considered Heavy Neutral Leptons and illustrated for which masses and lifetimes the aforementioned effect occurs.We show this in the top panel of Fig. 5 for the case of pure mixing with tau neutrinos.Such particles can decrease N eff for masses already above ∼70 MeV, even if they inject most of their energy into neutrinos.HNLs that mix purely with electron or muon neutrinos can only decrease N eff , as the parameter space in which they increase N eff is already excluded by BBN, CMB and experimental constraints [21,48].Therefore, in the pure mixing cases, only HNLs that mix with tau neutrinos have an unconstrained region of parameter space left where ∆N eff can be positive (bottom panel of Fig. 5).Such HNLs can increase N eff by at most ∆N eff ≈ +0.4, as higher values are excluded by BBN and the CMB.Given this maximum allowed value of ∆N eff , we then estimated what current-day Hubble rate it would correspond to.Using the Planck 2018 baseline TTTEEE+lowE analysis, we find that HNLs can increase the Hubble rate to at most H 0 = 70.5±0.7 km s −1 Mpc −1 and therefore moderately alleviate the Hubble tension.The simple model described in Sec.II relies on the assumption that the remaining fraction 1− of the injected neutrino energy is perfectly thermal.In reality, this may not be the case and the full thermalisation would occur during a much larger number of interactions than N therm log 2 (E inj ν /3.15T ).Therefore, this simple model underestimates the energy fraction that goes into the EM plasma 12 .The remaining non-equilibrium neutrinos will manifest themselves as residual non-thermal spectral distortions in the distribution function of neutrinos.These spectral distortions keep the energy exchange balance of ν ↔ EM reactions shifted to the right till long after FIMP decay.As a result, more neutrino energy will be transferred to the EM plasma and N eff can further decrease.There is a subtlety here that the remaining 1 − nonequilibrium neutrinos are only slightly hotter than the thermal neutrinos, and we cannot describe their thermalisation as an instant process: The corresponding rate is comparable to the thermal energy exchange rate.As such, the energy transfer process is extended in time, and a proper study of this effect requires solving the Boltzmann equation for the neutrino distribution function.
To study the impact of neutrino spectral distortions on the ν → EM energy balance shift, we consider a simple scenario where high-energy neutrinos are instantly injected into the primordial plasma.We make use of the publicly available Boltzmann code pyBBN13 [21] to simulate this process and to track the evolution of the neutrino distribution functions.Within this setup, neutrinos with energy E inj ν = 70 MeV are instantly injected at T = 3 MeV.They amount for a fixed percentage of the total neutrino energy density and are equally distributed over the three neutrino flavours.All Standard Model interactions as specified in [21] are included, but with neutrino oscillations turned off (without any loss of generality).In order to highlight the importance of neutrino spectral distortions, we perform this procedure a second time, but with neutrino spectral distortions turned off.In that case, the neutrino distribution function is given by a Fermi-Dirac distribution with temperature , where ρ να and g να = 2 are the energy density (of both neutrinos and anti-neutrinos) and number of degrees of freedom of neutrino flavour α respectively.
The evolution of the ratio ρ νe /ρ EM (relative to the one in the SM) is shown in the left panel of Fig. 7 for different amounts of injected neutrino energy.In agreement with the story in Sec.II, we observe a fast drop-off in the ratio right after the injection, which signifies the quick transfer of energy from the neutrino plasma to the EM plasma.After reaching the SM value (which naively corresponds to an equilibrium state), the ratio continues decreasing.This is the effect of the extended thermalisation due to neutrino spectral distortions, as caused by the remaining fraction 1 − of non-equilibrium neutrinos.Eventually, the ratio will be smaller than the SM value and ∆N eff becomes negative.In this plot, the dashed lines correspond to the same simulations but with a thermal-like distribution for the neutrinos.It is clear that without spectral distortions, the energy transfer from the neutrino sector to the EM sector is much less efficient.
Another way to look at this shift in the energy transfer balance from the neutrino plasma to the EM plasma is to ask the question: Which temperature T EM,eff is the EM plasma trying to reach after the injection?As we will see, depending on whether neutrinos have a non-equilibrium or a thermal-like distribution, this temperature can be either larger than or equal to the neutrino temperature 14 .In the former case, it means that the EM plasma temperature can exceed the neutrino temperature (and thus ∆N eff can be negative), while in the latter case ∆N eff cannot be negative.
In more technical terms, the exchange of energy between neutrinos and EM particles is regulated by the Boltzmann collision integral I coll , which encodes all interactions between the species.For neutrinos that participate in reactions of the form ν+2 ↔ 3+4, the collision integral is given by [73]: where f i and P i are the distribution function and four momentum of species i respectively, and |M| 2 is the unaveraged squared matrix element summed over degrees of freedom of initial and final states.The energy transfer rate between the neutrino and EM plasma can be written as: where we consider I coll to be a function of the EM plasma temperature T EM .There exists a temperature T EM,eff 14 In all cases, with 'neutrino temperature' we refer to the quantity , where gν = 2 and ρν is the energy density of both neutrinos and anti-neutrinos.
T ν (thermal) for which this rate is equal to 0. This corresponds to the temperature the EM plasma tends to during thermalisation, since then the system would be in equilibrium.In the case where neutrinos would have a thermallike spectrum with temperature T ν , the rate vanishes when T EM,eff = T ν .On the other hand, when a nonequilibrium neutrino spectrum is considered, we find that T EM,eff > T ν when Γ = 0.In the former case N eff cannot decrease, while in the latter case the EM plasma temperature can exceed T ν and N eff can thus decrease.We show the evolution of T EM,eff and T ν as obtained from the instant neutrino injection simulations in the right panel of Fig. 7.
The conclusion here is that neutrino spectral distortions play a central role in transferring energy from the neutrino sector to the EM sector.When considering shortlived FIMPs that can decay into neutrinos, the impact of these distortions on the evolution of N eff should not be neglected.After our work was submitted, the paper [27] appeared that studies the impact of HNLs with masses m N < m π on N eff .The authors of this work used numerical simulations in order to obtain N eff and disagree with our conclusion that N eff can decrease even if most of the HNL energy is injected into neutrinos.They have presented an analytic argument in their Appendix C which aims to demonstrate that our conclusion on N eff is wrong.They start with a toy model in Eq. (C.1) that describes the evolution of the distribution function of neutrinos f ν : D1) where x = ma (with a the scale factor and m = 1 MeV), H is the Hubble rate, ς is a constant and S(x) > 0 is the source term from decays of HNLs.The second term in the brackets describes the interactions between neutrinos and EM particles, where f eq is the equilibrium distri-bution function resulting from the interaction dynamics of neutrinos and EM particles in the presence of HNLs.Their argument as to why N eff cannot decrease goes as follows: as far as the source injecting rate S(x, E ν ) and the collision rate G 2 F T 4 E ν are much higher than the Hubble rate, the solution of Eq. (D1) may be given in terms of the quasi-static solution: In the limiting case S G 2 F T 4 E ν , the solution is just f ν = f eq , while in the opposite case f ν f eq .The authors conclude that in any case f ν ≥ f eq and thus ∆N eff ≥ 0. However, while this argument may be applicable at very early times when neutrinos are in perfect equilibrium, it is no longer valid at temperatures T = O(1 MeV), when they start to decouple.During decoupling process, the dynamics of the equilibration between neutrinos and EM particles, i.e., the energy transfer between the two sectors, becomes very important and is not captured by Eq. (D1).
We reiterate our argument as to why N eff can decrease when FIMPs inject most of their energy into neutrinos, but now from the point of view of the neutrino distribution function (see also the right panel of Fig. 7 and the surrounding text for a similar discussion).Before the decay of the FIMP, the neutrino distribution function is the same as the equilibrium distribution, f ν = f eq .Right after the decay of the FIMP, the neutrino distribution at high energies becomes f ν > f eq , while at low energies it is still f ν = f eq .During the thermalisation, highenergy neutrinos interact with both low-energy neutrinos and EM particles.In this process, the temperature of the equilibrium distribution function f eq increases.Now, neutrinos in the high-energy tail of f ν interact efficiently, see Eq. ( 5), and f ν −→ f eq for such neutrinos.But at low energies, neutrinos do not interact efficiently anymore to catch up with the increase of f eq , which eventually leads to f ν < f eq in this energy range.Given that these lowenergy neutrinos contribute the most to N eff , it means that ∆N eff can become negative.

Figure 1 :
Figure1: Estimate of the fraction of injected neutrino energy (both thermal and non-equilibrium) that gets transferred to the EM plasma during thermalisation (see text for details).The three curves indicate the value of when a neutrino of energy E inj ν is injected at a temperature Tinj.At high temperatures of order of Tinj E inj ν , the injected neutrinos are thermal-like, and hence is small.Once the temperature decreases, we enter the regime E inj

Figure 2 :
Figure 2: ∆N eff as a function of the lifetime of a FIMP χ that can only decay into neutrinos through χ → νe + νµ + νµ.The initial FIMP abundance is assumed to be nχ/s = 0.01 at T = 1 GeV, where s is the total entropy density of a universe consisting of photons, neutrinos and electrons/positrons.The solid lines are the result of our semi-analytic model, while the dotted lines are obtained with the Boltzmann code pyBBN.The grey band is the current sensitivity by Planck (see text for details).The golden curves roughly indicate the lowest FIMP mass for which N eff can decrease due to the thermalization of the injected neutrinos.The stronger decrease of the blue, dotted curve as compared to the solid curve highlights the significance of residual neutrino spectral distortions in the evolution of N eff (see Appendix A for more details).

Figure 3 :
Figure3: The fraction of HNL mass that is injected into the neutrino plasma.Contributions to this fraction from unstable HNL decay products (mesons and muons) are included and we assume that the kinetic energy of all created charged particles goes into the EM plasma.For mN 135 MeV, HNLs can decay into neutral pions, which in their turn decay into two photons.This causes the sudden decrease of ξν around that mass.At higher masses, ξν keeps increasing in the case of τ -mixing, which is due to the absence of HNL decays into charged mesons (such decays are possible in the other two mixing cases).

Figure 4 :
Figure4: ∆N eff as a function of HNL lifetime for a number of benchmark masses.Mixing with electron neutrinos only is considered here.The curves illustrate three cases of how HNLs can affect N eff : 1) they can decay mostly into neutrinos and simply increase N eff (30 MeV curve), 2) they can decay mostly into neutrinos and either decrease or increase N eff depending on their lifetime (110 MeV curve), and 3) they can decay mostly into EM particles and simply decrease N eff (200 MeV curve).HNLs with masses mN 70 MeV that decay mainly into neutrinos around neutrino decoupling, show an initial decrease of ∆N eff as a result of the thermalisation of the injected high-energy neutrinos.The grey band is the current sensitivity by Planck.

4 − 0 . 4 Figure 5 :
Figure5: How HNLs change ∆N eff as a function of their mass and lifetime.Mixing with tau neutrinos only is considered here.Top panel : Regions of the HNL parameter space that predict an increase (blue) or decrease (red) of N eff with respect to the SM value.The horizontal lines at the bottom of the plot indicate the mass ranges where HNLs inject most of their energy into neutrinos (ξν > 0.5) or the EM plasma (ξν < 0.5).In the former case, HNLs can still decrease N eff as a result of the efficient transfer of energy from neutrinos to EM particles.Bottom panel : Regions of the HNL parameter space that are excluded by BBN abundance measurements (green) and CMB observations (yellow).The ∆N eff = {0, ±0.4} contours give an indication of by how much HNLs can change N eff at the most in the unconstrained region.The BBN bound is from[21] and uses a central value for the primordial helium abundance of YP = 0.245[51] with an error of 4.35% (see[48] for a discussion on how this error is obtained).For masses higher than the eta-meson mass (∼550 MeV), the meson effect from[48] is included in the analysis.The CMB constraint is obtained using the approach as detailed in[21] (see Sec. III B for more details on the CMB bound).This panel also shows that there is only a relatively small unconstrained region of parameter space left that can increase N eff and where HNLs could play a role in alleviating the Hubble tension.

Figure 6 :
Figure6: Semi-analytic estimate of ∆N eff as a function of HNL mass and lifetime in the case of pure tau mixing.This plot is obtained using the method described in Sec.II (and is therefore only accurate up to a factor 3 − 4 for short lifetimes, when neutrinos are still in partial equilibrium).Nevertheless, it allows for a qualitative understanding of the behaviour of ∆N eff at lifetimes larger than considered in the main analysis (Fig.5).Importantly, for lifetimes well beyond the time of neutrino decoupling (O(1) s), non-equilibrium effects are absent and the sign of ∆N eff is thus completely determined by the fraction of HNL energy ξν that is injected into the neutrino plasma, see Fig.3.We see that HNLs with low masses and long lifetimes can still considerably affect N eff , while in the higher mass range ∆N eff tends to 0. This is because lowmass HNLs are more abundantly produced in this region of parameter space[26], where their mixing angles are relatively large.
Appendix A: Effect of Residual Non-equilibrium Neutrino Distortions

Figure 7 :
Figure 7: Evolution of the neutrino and EM plasma after the instant injection of neutrinos with energy E inj ν = 70 MeV at T = 3 MeV.Left panel : The ratio of electron neutrino energy density to electromagnetic energy density, relative to the SM prediction.Three fractions of the injected energy density are considered: ρ inj νe /ρ tot ν = {0.2%,1%, 5%}.The solid lines are obtained by taking into account the full non-equilibrium spectrum of neutrinos, whereas the dashed lines correspond to the evolution assuming that neutrinos always have a thermal-like spectrum with temperature Tν ∝ ρ 1/4ν .Right panel : Evolution of the neutrino temperature (dashed) and effective EM plasma temperature (solid) for which the energy transfer rate in Eq. (A2) vanishes.An injected fraction of ρ inj νe /ρ tot ν = 5% is considered here.The solid and dashed lines indicate when non-equilibrium and thermal-like neutrino distributions are used respectively.
m N is the HNL mass in MeV and τ N is the HNL lifetime in seconds.The change in N eff is with respect to the SM value of N SM eff = 3.026.The fitting functions are tested for masses 100 MeV ≤ m N ≤ 1 GeV and lifetimes 0.02 s ≤ τ N ≤ 0.05 s, and have a maximum deviation from the simulated data of roughly ∼3%.
Appendix D: Comment on "Massive sterile neutrinos in the early universe: From thermal decoupling to cosmological constraints" by Mastrototaro et al.