Massive sterile neutrinos in the early universe: From thermal decoupling to cosmological constraints

We consider relatively heavy neutrinos $\nu_H$, mostly contributing to a sterile state $\nu_s$, with mass in the range 10 MeV $\lesssim m_s \lesssim m_{\pi} \sim 135$ MeV, which are thermally produced in the early universe in collisional processes involving active neutrinos, and freezing out after the QCD phase transition. If these neutrinos decay after the active neutrino decoupling, they generate extra neutrino radiation, but also contribute to entropy production. Thus, they alter the value of the effective number of neutrino species $N_{\rm eff}$ as for instance measured by the cosmic microwave background (CMB), as well as affect primordial nucleosynthesis (BBN), notably ${}^4$He production. We provide a detailed account of the solution of the relevant Boltzmann equations. We also identify the parameter space allowed by current Planck satellite data and forecast the parameter space probed by future Stage-4 ground-based CMB observations, expected to match or surpass BBN sensitivity.


I. INTRODUCTION
Feebly interacting particles characterised by extremely suppressed interactions with the Standard Model particles have received growing interest in the last decade (see [1] for a recent review). In this context, a fourth neutrino mass state ν H with mass ∼ O(100) MeV, mostly contributing to an electroweak singlet neutrino state ν s due to Z-width constraints [2], emerges rather naturally in extensions of the Standard Model, like dynamical electroweak symmetry breaking [3] or the Neutrino Minimal Standard Model (νMSM) [4,5]. In the latter case, such particles can be related to fundamental problems of particle physics like the origin of neutrino mass, the baryon asymmetry in the early universe and the nature of dark matter.
The parameter space of a fourth neutrino in this mass range is strongly constrained by collider and beam-dump experiments for a dominant mixing with either ν e and ν µ [6,7], but it is significantly less constrained if mixed with ν τ , with bounds at high masses coming from searches of decays of D mesons and τ leptons [8] and SuperKamiokande data [9]. Furthermore, ν H can be emitted by a core-collapse supernova. In this context, limits have been placed from the SN 1987A observation, requiring the SN core may not emit too much energy in the ν s channel, since this additional energy-loss would shorten the observed neutrino burst [10][11][12][13]. Additionally, in [12,14] it has been discussed a possible role of heavy sterile neutrinos in enhancing supernova explosions.
Further and complementary constraints on heavy sterile neutrinos can also be placed from cosmological arguments. Indeed, ν s can be produced in the early universe via collisional processes involving active neutrinos, and then decay into lighter species. In particular, for masses in the range 10 MeV m s m π ∼ 135 MeV the main decay channels are ν s → ν ανβ ν β (with branching ratio of O(60 − 90%), depending on the mixing) and ν s → νe + e − (with branching ratio of O(10 − 40%), depending on the mixing). The decay products of the sterile neutrinos are injected into the primordial plasma, with the timescale of the event determining its phenomenological impact. If the decay is over when active neutrinos are still tightly coupled to the electromagnetic (e.m.) plasma dominated by photons and e ± pairs, full equilibrium conditions are quickly established and no effect remains, but for an unobservable renormalisation of the baryon to photon ratio η 1 . If significant decay happens after the active neutrinos have decoupled, part of the energy injected will end up in extra neutrino radiation and part heats the e.m. plasma up. The latter adds to the eventual heating due to e + e − annihilation and increases the photon to neutrino temperature beyond its standard value of T γ /T ν (11/4) 1/3 . Together with the former effect, this alters the effective number of neutrino species N eff , with the two processes going in opposite directions. Also, the non-thermal ν e andν e spectra enter weak interactions, alteringtogether with N eff -the neutron-to-proton ratio which rules the abundance of the primordial yields [10,11,[15][16][17][18], affecting in particular the 4 He abundance encoded in the primordial Helium mass fraction parameter Y p .
The aim of our paper is to perform a detailed calculation of heavy sterile neutrino decoupling in the early universe, in particular computing the effects on N eff and Y p and assessing the impact of the approximations presented in the seminal works [10,11]. Compared instead to the most recent calculations such as [18,19], we put more emphasis on the low-mass sterile neutrino range, where the number of effects is limited, hopefully offering a pedagogically complete and contained treatment, besides clarifying some points of disagreement. At high masses, new channels involving pions open up and are responsible for extra effects on primordial nucleosynthesis as well as N eff (some of these interesting effects have been described in [17][18][19]). In the same spirit, we will limit ourselves to neutrinos that decouple after the QCD phase transition, which translates on the chosen parameter space. One of our primary goals is to compare the effect on BBN with the latest constraints by the Planck satellite experiment as well as forecasts of future Stage-4 (S4) ground-based CMB observations [20]. We shall illustrate the shifting cosmological constraining power, from a BBN-dominated one to a CMB-dominated one, expected to be basically completed by the S4 era.
The plan of our work is as follows. In Sec. II we present the heavy sterile neutrino model we will use as a benchmark. In Sec. III we discuss and solve the kinetic equations describing the sterile neutrino evolution in the early universe. In Sec. IV we characterise the impact of heavy sterile neutrino decays on active neutrinos and on N eff . In Sec. V we present the current constraints and forecasts on sterile neutrino parameter space from BBN and CMB data. Finally, in Sec. VI we summarise our results and conclude. In Appendix A we report useful analytical approximations of the decay and scattering rates involving sterile neutrinos. Appendix B details the steps involved in the dimensional reduction of the collision integrals used in the numerical integrations. Appendix C is devoted to a comparison of our results with others previously reported in the literature.

II. HEAVY STERILE NEUTRINO MODEL
We consider heavy sterile neutrinos with masses 10 MeV m s 135 MeV, mixed dominantly with one active neutrino ν α (α = e, µ, τ ) as ν α = cos θ αs ν + sin θ αs ν H , ν s = − sin θ αs ν + cos θ αs ν H , where ν and ν H are a light and a heavy mass eigenstate, respectively, and θ αs 1, i.e. ν is mostly active and ν H is mostly sterile. We can relate the mixing angle to the unitary mixing matrix U , where Through neutral-current interactions, ν H can decay into a ν and a pair of other light leptons. It can also scatter with other species in the plasma. With a little abuse of notation, we shall refer to ν H as ν s and to ν as ν α . Unless stated otherwise, we shall consider neutrinos to be Dirac particles. In vacuum, in the mass range of our interest, the decay rate of sterile neutrinos is dominated by three lepton final states. In what follows, we shall neglect terms of the order m ν /m W,Z in the matrix elements. Their relevant decay processes and matrix elements are presented in Table I. Unless stated otherwise, we adopt the case of mixing with ν τ as our benchmark. In some cases, we will consider mixing with ν e in the mass range m s < m e + m µ in order to compare with previous literature. 2 The decay rate of sterile neutrinos into three neutrinos (summed over all flavours) is given by (see e.g. [21]) I: Squared matrix elements for sterile neutrino decay processes (assuming mixing with the species α, and β = α), summed over initial and final states and divided by the spin dof of the sterile neutrino. The particles involved in each decay are enumerated as 1 → 2 + 3 + 4. In the last lines,gL is replaced by gL in case of mixing with νe. The symmetry factor S = 1/2! is included, when two identical particles are present in the final state (first row). Process TABLE II: Squared matrix elements for sterile neutrino scattering processes (assuming mixing with the species α, and β = α), summed over initial and final states and divided by the two spin dof of the sterile neutrino. The particles involved in each reaction are enumerated as 1 + 2 → 3 + 4. In the last line,gL is replaced by gL in case of mixing with νe. The symmetry factor S = 1/2! is included, when two identical particles are present in the final state (second row). Process while the decay into neutrino plus e + e − pair, in the limit where m e /m s 1 is neglected, is where and, in case of mixing with ν e ,g L is replaced by g L . As a result, the total decay width writes Note that for Majorana neutrinos, the widths for the exclusive decay processes would be the same as for the Dirac case, but the inclusive decay is a factor of two larger, since for each final state accessible to a Dirac particle, two channels are present (i.e. a channel and its L-number conjugate).
Two-body reactions affect the sterile neutrino chemical and kinetic equilibrium. Those of interest for us are reported, together with the relevant squared matrix elements, in Table II. Our results agree with those reported in table 5 and  table 6 of Ref. [18], bearing in mind their different definition of |M | 2 , just summed over helicities of initial and final states, instead of averaged over the initial state as in our case. 3 In relation instead to the Ref. [10,22], we find systematically lower matrix elements compared to their Tables 1, 2, by a factor 2. We believe that this is likely a typo in the quantities reported rather than in their results since we do find agreement with the rates they compute.

A. Equations of motion
Following [23], in order to describe the time evolution of the sterile neutrino ensemble in the early universe, it proves useful to define the following dimensionless variables which replace time, momentum and photon temperature, respectively where m is an arbitrary mass scale which we set equal to 1 MeV. Note that the function a can be normalized, without loss of generality, so that z = 1 at the largest temperature of interest here, when all particles in the plasma are in equilibrium with each other, and neutrinos also share the same temperature T . In terms of these variables, we can write the equations of motion (EoMs) for the heavy sterile neutrinos distribution function f νs as [10,11] A similar equation holds for the evolution of active neutrino species f να : In the previous expressions, H denotes the cosmic expansion Hubble rate given by the Friedmann equation as H 2 = 8πρ/(3m 2 Pl ) where ρ is the total energy density and m Pl = G −1/2 N is the Planck mass in terms of the Newton constant G N . We will consider the plasma to be initially thermally populated by pions and all lighter particles, neglecting the nuclei contribution to ρ and assuming equal distributions of particles and antiparticles.
The right-hand-side (r.h.s.) term in Eq. (8)-(9) contains the sum of the collisional and decay terms for sterile and active neutrinos, respectively and reads with |M f i | 2 the sum of the squared-matrix elements over initial and final states, divided by the spin multiplicity of the state of interest. In the case of I[f νs ], it contains decay and scattering processes, as shown in Table I and II, respectively. In the case of I[f νa ], it contains scattering processes analogous to those of Table II, apart for the replacement m s → 0 and |U τ s | 2 → 1 − |U τ s | 2 or 1, depending if one is dealing with the mixed flavour or not; it also contains the sterile neutrino decay source term, which is calculated referring to the processes in Table I. We use squared-matrix elements for the collisional processes for active neutrinos consistent e.g. with the results of Table 3 of [18], modulo our different definition of |M f i | 2 , while we find once again that the results reported e.g. in [10,22] are a factor two larger than the correct ones. The statistical factor writes in general as where f i,f are the distributions of the particles in the initial (i) or final (f ) states, the (-) sign ("blocking") refers to fermions and (+) sign applies to bosons ("stimulated" effect). Only fermions are however present in the processes of interest for us. Due to the different timescales of the active neutrino oscillation processes (see e.g. [24]), we take into account oscillations among active neutrinos by post-processing the flavour distributions via where P αβ are the time-averaged transition probabilities (see e.g. Eq. (26) in [17]). Medium modification of the mixing parameters are also included, similarly to ref. [18]. To get the time evolution of sterile and active neutrino distributions, we have to complete the set of equations [Eqs. (8)-(9)- (23)] with the continuity equation, stating the conservation of the total energy density whereρ andP are the comoving energy density and pressure of the primordial plasma, respectivelȳ Note that only massive components in the plasma contribute to the r.h.s. of Eq. (13), while (neglecting plasma corrections) relativistic species at equilibrium cancel out. From Eq. (13) by using the expressions in [23], and remembering that we consider equal particle and antiparticle distributions, where i runs over the three pions, π + , π − , π 0 , and is either the muon or the electron. Note that the Hubble function can be now expressed as x 2 (ρ γ +ρ e +ρ µ +ρ π +ρ νa +ρ νs ) 1/2 .
Equation (13) gets contributions from all species, and can be recast into the equation for the z(x) relation. Let us specify the different contributions. For photons one has: For the electrons, if setting κ e ≡ m e /m, one finds: where the functions F ± 1 and F ± 2 are defined as with the signs ± that take into account the boson and fermion nature of the particle respectively. A similar expression holds for muons, with κ e → κ µ . For each pion species i, one has instead: For sterile neutrinos (κ s ≡ m s /m): For active quasi-massless neutrinos, if α = e, µ, τ , one has Down to a few MeV, the active neutrinos are coupled to the rest of the plasma, which means that at sufficiently high temperatures (low x) we can write f να = 1/(exp(y/z) + 1) and As a result, in computing the z = z(x) relation we can save considerable computer time by considering two different regimes: Eq. (30) for x < x d , while numerically computing f να from the Boltzmann Eq. (9) for x > x d , where x d represents any epoch before neutrino decoupling, but otherwise arbitrary. In terms of the step function Θ, collecting all terms for photons, electrons, pions, active and sterile neutrinos and isolating dz/dx, Eq. (13) can be written as: where we defined: and Together with z = 1 as initial condition, Eq. (31) gives the "time-temperature" evolution. Provided that x d is sufficiently small, roughly x d 0.2 (i.e. T 5 MeV), the computed behaviour is insensitive to the choice of x d , as we illustrate in Fig. 1, where the extra comoving neutrino energy density evolution is computed using x d = 0.1 (red solid curve) and x d = 0.2 (black dotted curve), for parameters m s = 100 MeV and τ s = 0.045 s. The results are almost equal except for small numerical differences when 0.1 < x < 0.2. In the following, we fix x d = 0.1. Note that the parenthesis at the r.h.s. of Eq. (31) describes the "heating" of the e.m. plasma due to the sterile neutrino entropy release: At early times, all of the decay products end up in the coupled plasma of photons and active neutrinos, raising z compared to standard expectations. At late times, a sizable fraction decays into (decoupled) active neutrinos, hence the term in parenthesis largely cancels out. The only non-trivial evolution of z is then due to the finite mass term of e ± , affecting their annihilation at late time.

B. Evolution of heavy sterile neutrinos
In [11], an analytical solution of Eq. (8) was provided under the following assumptions: (i) The equilibrium distribution functions "inside" the collisional integral [Eq. (10)] are taken in the Boltzmann (classical) approximation, with Pauli blocking factors correspondingly neglected.
(ii) Electrons are considered ultrarelativistic, with terms in m 2 e neglected throughout. We repeated the derivation of [11] under these approximation, finding: and in agreement with their results quoted as where f eq νs is the Fermi-Dirac equilibrium distribution of the sterile neutrinos, τ s is the sterile neutrino lifetime, and E s = m 2 s + (y/x) 2 is the sterile neutrino energy. Details of the reduction of integrals in Eq. (10) under approximations (i) and (ii) are given in the Appendix A.
We also solve the sterile neutrino kinetic equations numerically, relaxing the approximations (i) and (ii) in Eq. (10), for both sterile neutrinos and active neutrinos. Following the well-known technique developed in [25], it is possible to analytically reduce the nine-dimensional collisional integral into a two-dimensional one, which is then integrated numerically. We developed an equivalent technique for the decay processes. We report the details in Appendix B.
Our solutions are obtained by assuming initial thermal equilibrium for all species, starting from a temperature T = min[2 m s , 150 MeV]. To compare the difference in using numerical results vs. the analytical approximation, in Table III we report the sterile neutrino freeze-out temperature T D , for a few representative points in parameter space, according to the condition I(T D ) = H(T D ), with I given by Eq. (10). We find typical differences at a few percent level, and in all cases below 10%. Although we are using the numerical results in the following, when requiring only moderate precision on the sterile decoupling, the analytical approximation seems largely sufficient, and allows one to significantly gain in computing time.

IV. IMPACT ON COSMOLOGICAL OBSERVABLES
After the distribution functions and temperature evolution are found, we relate them to the observables N eff (notably at the CMB epoch) and Y p (notably at the BBN epoch) to derive some constraints.
A. Impact on effective number of active neutrinos N eff Heavy ν s affect the total energy density in non-electromagnetic species. This is usually quantified in terms of the effective number of neutrinos, N eff , which is defined from the density in all species but electromagnetically interacting ones, as (see for instance [24]) where the r.h.s. is specific for our 4 neutrino model, where ∆ρ να are the changes in the neutrino energy densities with respect to ρ ν0 , the energy density in the instantaneous decoupling limit, due to the non-equilibrium effects. Note that at early times around T 100 MeV when all species are relativistic and share the same temperature, N eff → 4. Asymptotically, when all sterile neutrinos have disappeared and the e + e − annihilation is complete: i) z 0 → (11/4) 1/3 1.4, the asymptotic Standard Model photon-neutrino temperature ratio in the instantaneous  decoupling limit; ii) z → z fin , the actual final photon/neutrino temperature; iii) ρ νs /ρ ν0 → 0, since all sterile neutrinos have decayed away. At large x, we expect N eff 3 since, due to the branching fraction, the contribution due to extra radiation in the neutrino sector largely compensates the entropy transfer, which would tend to lower N eff below 3 via the z-dependent pre-factor at the r.h.s. of Eq. (37). Indeed this behavior can be seen in Fig. 2 where we have plotted the N eff evolution for m s = 30 MeV, τ τ s = 0.15 s (in solid red) and m s = 100 MeV, τ τ s = 0.055 s (dashed black), assuming mixing with ν τ . Note how the more massive neutrino decays deeper in its non-relativistic regime and well after its decoupling: The initial decline is due to the still-coupled massive neutrino, experiencing 'Boltzmann suppression'. Soon later, it decouples and its contribution to the plasma rises, with N eff that follows this trend. The growth is somewhat less steep than naively expected in absence of decays since this process partially counteracts the sterile neutrino energy density relative growth. This growth turns into a sharp decline around x 0.25, when decay takes over. If all its entropy were transferred to the active neutrinos only, N eff would stay constant. The partial redistribution to the e.m. plasma causes a minor decline in N eff , basically complete by x ∼ 1. Modulo quantitative differences, the lighter and longer lived neutrino depicted with the red line follows the same stages but shifted to the right since it stays coupled longer and decays later. The contribution of the sterile neutrino to this dynamics is more clearly visible in Fig. 3, while Fig. 4 shows the entropy dilution effect entering z, which partially counteracts the extra energy density in neutrinos and also affects N eff . All other parameters being the same, the effect on z is more pronounced, as expected due to the larger branching ratio in e.m. species. Note that a further, "standard" enhancement in z, due to e ± annihilation, happens at x 1 and is not visible in the figure. In Fig. 5 we show the evolution of the active neutrino energy densities. Note how they already depart from equilibrium somewhat by x ∼ 0.2 (for m s = 100 MeV), when the major enhancement happens due to the bulk of the decays.

B. Impact on Yp
Another important parameter affected by a massive sterile neutrino scenario is the Y p value, i.e a proxy for the primordial 4 He mass fraction. The effect arises due to the modified expansion history, i.e. via the role that N eff and z(x) play in the Hubble function H(x), but above all via the distortions to the electron (anti)neutrino distribution entering the isospin changing reactions between neutrons and protons. In Fig. 6, we illustrate the typical ν τ and ν e spectra (ν µ being intermediate between the two) associated to heavy sterile neutrino decays, in the case of mixing with ν τ : Despite mixing among the active species, the largest distortion remains in the ν τ species; also, heavier neutrinos lead to more energetic residual distortions. Alterations due to heavy sterile neutrino decays also affect other nuclei such as deuterium, but these are subleading compared to the effect on Y p (see e.g. Fig. 3 in [18]) and for simplicity we will limit ourselves to model the modifications on Y p . Both CMB and BBN are sensitive to Y p , but astrophysical determinations of Y p and thus a comparison with primordial nucleosynthesis predictions is currently more constraining.
A precise standard model calculation, Y prec p,SM , for the best-fit cosmological parameter Ω b h 2 = 0.02225 [26] with a number of subtle effects included (see [27,28] for details), is available from Parthenope [29,30]. This result is rescaled via the ratio of the Born estimate of the Y p for the ν s model, Y Born p,νs , over the Born standard model calculation Y Born p,SM , as Each term of the fraction at the r.h.s. can be estimated as (see e.g. [31]): where t on 180 s corresponds to the onset of the BBN (i.e. deuterium bottleneck opening around T 0.08 MeV) 4 , τ n is the neutron lifetime and ω B are the rates in the Born approximation of the processes in Table IV (with ∆ ≡ m n − m p 1.29 MeV) that can be written as: The distributions entering in F , as well as H(x), are taken from the numerical solutions of our system of equations. In the last column of Tab. III, we show the results for some relevant sterile neutrino parameters.

V. CONSTRAINTS AND FORECASTS
In order to obtain constraints on heavy sterile neutrinos, we compare our results on N eff and on the modification on Y p with both the latest CMB and BBN measurements. For BBN, we use the current bound at 2σ [32] Y p = 0.245 ± 0.006 .
Concerning CMB, if limiting oneself to N eff , the latest measurements of the Planck collaboration provide a value N eff = 2.99 ± 0.17 [26]. Therefore we could exclude at 2σ extra-radiation leading to ∆N eff > 0.33. In practice, the massive sterile neutrino model under consideration here leads to changes in both N eff and Y p , and the CMB is sensitive to both (albeit much less to Y p than BBN, at the moment). Hence we infer the CMB constraints using a reduced Gaussian likelihood matrix involving N eff and Y p , of the form [18]: (σ 1 , σ 2 ) = (0.2650, 0.0177) , To obtain results at 2σ, we have to consider a value of χ 2 = 6.18, value obtained by requiring that the integral of the χ 2 -distribution with 2 dof is equal to 0.9545. Our results from CMB measurements and from BBN based on Y p are shown in Fig. 7 for the constraints on the decay time, and in Fig. 8 for the constraints on the mixing parameter, for the most interesting case of mixing with ν τ . As a side result and in order to allow for cross-checks with past literature, in Fig. 9 and Fig. 10 we also report the corresponding results for a mixing with ν e . Besides the current constraints,  we also show the sensitivity forecast of the future CMB-S4 observations, with uncertainties (σ 1 , σ 2 ) = (0.062, 0.0053) according to [33], considering the same Θ obs as in Eq. (45).
We conclude that the CMB provides already the best constraints for m s 50 MeV, while BBN takes over at larger masses. However, we expect that CMB-S4 will attain leading constraining power in the whole range of parameter space considered here, if performing close to expectations. Qualitatively similar consideration apply for mixing with ν e , although the transition mass is around m s ∼ 20 MeV, and the future improvement of CMB-S4 over BBN is less significant. Note that, in particular for the CMB, at the same mass and lifetime the bounds are more stringent for a mixing with ν τ than one with ν e : This is due to the fact that the bound is dominated by N eff and, due to the larger b.r. in neutrinos for the case of mixing with ν τ , the growth of neutrino density via non-thermal injection is only mildly compensated by the entropy effect. For the case of mixing with ν e , there is instead a substantial compensation via the growth of z.
In the case of the BBN bound, however, the leading effect is due to ν e distortions, which are larger when the mixing is with ν e ; the effect of N eff > 3 altering H is however more relevant when the mixing is with ν τ , so that the two constraints are closer to each other in this case.
Our BBN constraints are largely consistent with recent calculations presented in [18], while our CMB constraints are consistent with theirs at low masses, while more stringent than theirs at high masses. A more detailed comparison and discussion are reported in Appendix C, where we identify the origin of the difference in the estimate of ∆N eff . In particular, contrarily to the results of [18], we always obtain ∆N eff > 0 for the parameter space of interest. We have supplemented our numerical calculations with a qualitative study of the Boltzmann equation in the analytical approximation of [11] to further support our conclusions.

VI. CONCLUSIONS
Heavy sterile neutrinos with masses O(MeV-GeV) are predicted in extensions of the Standard Model such as the Neutrino Minimal Standard Model (νMSM). Besides affecting collider and supernovae observables, their mass and mixing angle parameters can be also constrained with cosmological observables, notably CMB and BBN.
We have numerically studied the evolution of sterile neutrinos with 10 MeV m s 135 MeV in the early universe and set constraints on the mixing angles or lifetimes using the N eff and Y p observables. In order to achieve these results, we have solved the exact Boltzmann equation for sterile and active neutrino evolution while taking into account the temperature evolution of electrons and photons. Also, we checked the correctness of analytical approximations in the literature and verified that they are adequate to describe sterile neutrino decoupling at better than 10% level.
For the least constrained (and thus phenomenologically most interesting) sector of mixing with ν τ , at m s 50 MeV these cosmological bounds surpass the traditional benchmark of 0.1 s lifetime often considered in the literature, up to about 0.03 s for the highest masses considered. While currently CMB is more constraining at low masses and BBN dominates at high masses, we expect the future CMB-S4 experiments to yield the dominant constraining power by the end of the decade, unless the systematic error affecting the astrophysical determinations of Y p can be significantly reduced. The work of L.M. is supported by the Italian Istituto Nazionale di Fisica Nucleare (INFN) through the "QGSKY" project and by Ministero dell'Istruzione, Università e Ricerca (MIUR). The computational work has been executed on the IT resources of the ReCaS-Bari data center, which have been made available by two projects financed by the MIUR (Italian Ministry for Education, University and Research) in the "PON Ricerca e Competitività 2007-2013" Program: ReCaS (Azione I -Interventi di rafforzamento strutturale, PONa3 00052, Avviso 254/Ric) and PRISMA (Asse II -Sostegno all'innovazione, PON04a2A).

Appendix A: Reduction of the decay and collisional integrals
In this Appendix, we show how, under the approximations adopted in [11], we can analytically solve the integrals appearing at the r.h.s. of Eq. (8), namely with |M f i | 2 the sum of the squared-matrix elements for the decay and scattering processes and Please note that our results are obtained with a numerical integration of the collisional integrals in Eq. (A1), after the dimensional reduction recapped in Appendix B, not with the analytical approximations. Also, None of the approximations reported below is new and can be skipped unless one is interested in reproducing the analytical approximations. However, for pedagogical purposes, we thought it is useful to report them here in great detail, not to force readers to go back to the decades-old original literature.

Decay Integral
Following the approximations of [11] we assume that in Eq. (A1) the energy distributions are represented by Maxwell-Boltzmann distributions, instead of Fermi-Dirac ones. Moreover, we neglect the Pauli blocking factor, assuming (1 − f i ) 1. Then one gets where the label 1 indicated the sterile neutrino and |M f i | 2 is the sum over the dominant decay processes with α = e, µ, τ . Performing the integral over d 3 p 4 in eq. (A3) using the delta function enforcing E 1 = E 2 + E 3 + E 4 , one obtains that e −E2/T e −E3/T e −E4/T = e −E1/T = f eq 1 . Moreover, using the property of the delta function in eq. (B1), we can write where the sterile neutrino decay rate Γ D is given by [10] Thus, we end up with where τ s = 1/Γ D is the sterile neutrino lifetime. We could also write the corresponding source term for an active species at the same level of approximation (i.e. neglecting the blocking factor and the inverse decay) following the treatment detailed in Ref. [13], which we address the reader to for details. This yields where B i is the branching ratio of the i-th exclusive reaction, γ = 1 + (y s m) 2 /(xm s ) 2 , β = y s / (m s x/m) 2 + y 2 s , and F a,i is the double-differential distribution (with respect to y and to the angular variable θ) of the daughter particle a in the reaction i in the sterile neutrino rest frame, normalized to 1. The integral kernel in the integral above acts as a constraint picking the "right" momentum for the daughter neutrino, weighting it for the occupation factor of the parent sterile species. In practice, we always use the full numerical integration rather than this approximate expression.

Collisional Integral
In order to evaluate the collisional integral in Eq. (A1), we assume below that they only mix with ν τ . The relevant collisional processes are shown in Table II and have a squared interaction matrix given by where C = 2 4 G 2 F |U τ 4 | 2 (1 +g 2 L + g 2 R ) (remember that in our definition of |M | 2 we sum over all the degrees of freedom and include the average over the relevant state). We indicate with 1 the sterile neutrino state, and evaluate the quantities in the center of momentum frame; since particles 2, 3, 4 are relativistic, and Thus Hence: Similarly, Hence: These equalities among Lorentz invariants hold in any frame. As a result, let us write: We thus have where θ and µ are the azimuth and polar angle respectively, cos θ 32 = − cos θ 31 due to the fact we evaluated the integral in the center of momentum frame, and we used Eq. (A29), with a = E 1 + E 2 and b = E 2 − |p 1 |, as well as Let us evaluate the integral over particle 2 (expressed in terms of Lorentz-invariants) by the explicit replacement p 1 · p 2 = E 1 E 2 − p 1 p 2 cos θ 12 : Similarly Finally, taking into account that the integrals for the C−term are symmetric under the relabelling 3 ↔ 4, we have with a similar procedure of Eq. (A16) (A20) Hence Summing all contributions, we have: Assuming in this last step Fermi-Dirac distribution for particle 2, following the procedure in Ref [11], which implies The above result agrees with what reported in ref. [11] in the same limit: with τ s given in Eq. (6).
Finally, we report below some relations used in the numerous integrations: so that, when p 2 1 = m 2 s and p 2 2 = p 4 3 = p 2 4 = 0, Also, we used some notable integrals: implying that implying that implying that sin β i is found as ±(1 − cos 2 β i ) 1/2 , where and Q ≡ m 2 1 + m 2 2 + m 2 3 − m 2 4 . The equation for cos β has two solutions, but we can account for them by multiplying by two and using as integration's interval [0, π]. The limits of integration in d cos α come from demanding that cos 2 β ≤ 1, meaning that (2p 2 p 3 sin α sin θ sin β) 2 ≥ 0 .
(B12) This is the same requirement that (df (β)/dβ) 2 ≥ 0. Therefore we can write Introducing the following definitions: the derivative can be written as: All possible matrix elements only include products of the four-momenta. All the products are analytically integrable over d cos α and can be carried out by using these relations: The step function comes from demanding a real integration interval. This also ensures that the roots of ax 2 + bx + c are not outside the fundamental integration interval of [−1, 1]. Integration over dµ is trivial because there is no dependence on this parameter. All the possible products of these momenta are calculated below: To integrate over d cos θ, the solutions of b 2 − 4ac are important for the the integration interval: If there is to be a real integration interval, both of these solutions must be real and we will refer to them as cos θ min and cos θ max . The real integration limits are α = sup[−1, cos θ min ] and β = inf[+1, cos θ max ] with α ≤ β. Finally, with these conditions, it is possible to calculate numerically the collision integral left: where A is the parameter space allowed, Λ comes from the following analytical integral: and F (f i , f f ) is the expression defined in Eq. (11) and evaluated for p 4 given in Eq. (B2). Finally, for the decay processes, it is possible to follow a similar procedure modulo the fact that we do not integrate analytically over cos α and we define: (B20)

Appendix C: Comparison with previous results
To validate our code, we have first compared our results with pioneering results obtained in ref. [11], finding excellent agreement with respect to several outputs, see e.g. Fig 11 showing the comoving density of the sterile species. The most recent constraints on heavy decaying sterile neutrinos have been reported in Ref. [18]. In Fig. 12 and 13 we compare BBN bounds from Ref. [18] with our BBN bounds, while in Fig.s 14 and 15 we compare CMB bounds from Fig. 11 in Ref. [18] with CMB results from our code, obtained using the same value of Θ Obs as a benchmark.   [11,18] and the one from our code.
While there is always a qualitative agreement, the quantitative agreement between the results is rather good only for the BBN case, while it shows some discrepancy in the CMB case for high masses. From inspection of ref. [18] (e.g., Section 4.1.2, Appendix B) we infer that the authors find a systematically lower value of N eff than us when masses are significantly large than ∼ 10 MeV, and that their N eff can also attain values below 3, while we always find ∆N eff 0 (see e.g. Fig. 2). This has been confirmed by private communication with the corresponding author of ref. [18]. This behaviour is also found in Ref. [17] (right panel of figure 3) and is indeed traced back by the authors of [18] to the solution scheme provided by the pyBBN code. As a result, in [18] CMB bounds at large masses are dominated by Y p (to which CMB is less sensitive) rather than by N eff , and are weaker than ours. Since BBN bounds from Y p depend mostly on spectral distortions of the electron-type neutrinos, it is not surprising that the agreement is much better  in this observable channel. For the case of mixing with ν τ , where N eff plays a slightly bigger role, the agreement is not as excellent while remaining good. Note that the inclusion of deuterium constraints is subleading to the effect on Y p (see Fig. 3 in [18]), so neglecting it has no major impact on our BBN bounds.
In the recent paper [34], a physical interpretation of the effect on N eff found in [18] is discussed: It is claimed that ∆N eff < 0 is a quasi-generic outcome of injection of energy after neutrino decoupling, even for decay modes mostly in neutrinos, as a result of a dominant entropy transfer to e ± via non-thermal neutrino interactions with particles of the thermal bath, see Eqs. 3, 4, 5 in [34]. Although the authors of [34] draw some analogy of their ∆N eff < 0 effect with the results of the pioneering work [35] on a low-reheating scenario, we believe that this is not very instructive, since the thermal neutrino bath is obviously suppressed if the universe starts at temperatures comparable or lower than the neutrino decoupling, a situation very different from the one under study.
Since we could not confirm numerically these surprising results, we thought useful to discuss their plausibility with a qualitative study of the Boltzmann equation, in the analytical approximation derived above, matching e.g. the one in [11]. If denoting with f 1 the non-thermal neutrino distribution under exam, it obeys an equation of the form: The first term S ∝ (f s − f eq s )/τ s f s /τ s at the r.h.s is always positive, since describing the injection due to decays of the sterile neutrinos. More precisely, S is an integral where (f s − f eq s ) enters as a kernel, see e.g. Eq. (A7). The second (collisional) term, where ς 2 is a positive numerical constant, is initially zero since neutrinos are at equilibrium, but as a result of the source term S, f 1 grows above the equilibrium value and leads to a negative value of the collisional term, linear in (f eq − f 1 ). A depletion of the neutrino distribution to "sub-thermal" values, as apparently found in [34], would imply reversing the sign of the collisional term, i.e. obtaining a positive value of (f eq − f 1 ). However, this second term being controlled by the source term S, it can at most grow negative to the point of compensating the first term, thus reaching an equilibrium between injection and collisional redistribution of the energy. At that moment, the derivative of f 1 is driven to zero and f 1 becomes constant. In practice, unless none of the rates is fast compared to the Hubble expansion H, the evolution follows one of the following paths: i) The first term at the r.h.s of Eq. (C1) dominates over the collisional term, which means that f 1 grows significantly larger than f eq , with the scattering incapable of fully compensating it.
ii) The second term at the r.h.s of Eq. (C1) is dominant: As a result, f 1 tends to f eq , annihilating the collisional term, or more precisely settling to a slightly larger-than-thermal value to compensate for the injection.
Hence, we conclude that either f 1 freezes out at a value larger than the equilibrium one (implying ∆N eff > 0), or in the 'worst' case it tends to the equilibrium distribution ( ∆N eff → 0 + ), in contradiction with the conclusions of ref. [34]. In the above discussion, we neglected the effect of the growing temperature as a result of the transfer of entropy from the sterile neutrino decays: At this level of approximation, however, its effect is to make the collisional term increase via the T 4 factor, as well as to increase f eq , so that the equilibrium distribution the particles are driven to is not the same as the initial one. For the sake of the argument, let us assume that as a result of collisions, the second term at r.h.s of Eq. (C1) starts becoming positive at some instant, i.e. f eq > f 1 : Then, f 1 would start growing again (because the r.h.s would be positive) and hence the r.h.s. is brought closer to zero. The feedback is such that any collisionally induced depletion of neutrinos would be immediately compensated, preventing a depletion of f 1 in the circumstances under exam. For that to happen, one needs a situation in which f eq grows due to the growth of temperature, with a minimal impact on f 1 . This can be easily obtained if the injection of energy from decay is dominantly in the e.m. sector and the collisional terms is negligible, i.e. a situation of type i) but in the e.m. sector, as naively expected. We checked that this is indeed the case, namely one obtains ∆N eff < 0 when artificially pushing the b.r. into neutrinos to sub-dominant values, and the decay happens sufficiently late. These regimes are illustrated in Fig. 16: If the sterile state decays early, the collisional term is capable of restoring equilibrium and ∆N eff → 0, as argued in the situation ii) described above, independently of the channel in which energy is injected. For a decay happening later and later, collisions are less and less efficient in restoring equilibrium and ∆N eff > 0 for a dominant b.r. into neutrino states, while ∆N eff < 0 for a dominant b.r. into the e.m. channel.
All our numerical results qualitative agree with these conclusions that can also drawn from Eq. (C1). Since Eq. (C1) is an approximation, one may wonder if the results of ref. [34] are due to some features not captured by Eq. (C1) (and, for some unknown reason, also missed by our numerical calculation). We believe that this is not the case as we argue in the following: • Eq. (C1) neglects quantum statistics effects, which appear however irrelevant to the arguments of ref. [34], and are anyway present for both electrons and neutrinos, without causing a qualitative change in one sector compared to the other.
• Eq. (C1) assumes T ν = T . However, this is not crucial to the conclusions above. The presence of a "two temperatures background" is essentially equivalent to split the term (f eq − f 1 ) into a linear combination of (f eq ν − f 1 ) and (f eq e − f 1 ), each weighted by a positive factor. Both these functions are negative under the effect of neutrino energy injection: One has a fortiori (f eq ν − f 1 ) < 0 since f eq ν < f eq e when T > T ν (and electrons are relativistic), while (f eq e − f 1 ) < 0 must be satisfied if we ask for T to grow larger than T ν as a result of injection of energy in the neutrino sector. The T evolution equation is indeed controlled by the opposite of the neutrino collisional term, i.e. goes as (f 1 − f eq e ), as can be checked via Eq. (31). • Eq. (C1) does not account for the quadratic terms involving the non-thermal parts of the neutrino distributions, i.e. depending from (f ν − f eq ν ) 2 , since thermal distributions have been used in the kernels of the collisional integrals, see Appendix A. Note however that these are not the processes claimed in [34] to be responsible for the effect on N eff , since their Eq.s 3, 4, 5 explicitly indicate reactions of the non-thermal neutrinos on thermal background species, which are included in the above treatment. Yet, let us entertain the possibility that the results outlined by the authors of [34] are physical and due to these non-linear effects, and that it is only their interpretation/attribution to be incorrect. One should then expect that ∆N eff < 0 shows a "threshold" behaviour with respect to the amount of energy injected in neutrinos: The less non-thermal energy is injected, the less likely these non-linear interactions among non-thermal particles should become when compared to interactions with the thermal background. Hence, the authors of [34] should have found that the effect kicks in only above some fraction of the background energy injected in the medium, growing quadratically above this value. Instead, not only they claim an effect also in the idealized case of 'single neutrino injection' (see their Fig. 1), but clearly show an effect that is roughly linear in the injected non-thermal energy (see Fig. 7, left: The change is from slightly below -0.2% to about -0.8% when moving from an injection of 1% to 5%), inconsistent with this hypothetical explanation.
In conclusion, to the best of our knowledge, we attribute our departure from the results in [18] on the CMB constraints at large masses to some unidentified systematics, probably the same effects responsible for the "universal" ∆N eff < 0 outcome described in [34] which appears unphysical, for the reasons detailed above. The toy model advocated in [34] to support their findings is also untrustworthy, since it does not account for reverse reactions (only reactions transferring energy from neutrinos into e.m. particles are included) and does not implement the physical requirement that only excess energy (above the thermal value) can be effectively transferred in collisions. As a consequence, it is for instance incapable of predicting equilibration and ∆N eff → 0 when energy is injected at early times.
Note added: After our paper was accepted for publication, a new version of ref. [34] appeared including an appendix D where the authors criticise the generality of the conclusions reported above. We believe that their arguments are incomplete and incorrect, as explained below.
First, they claim that our discussion based on Eq. (C1) and leading to the conclusion that ∆N eff ≥ 0 only applies to a situation close to equilibrium, i.e. at high-T . The authors seem to believe that we only consider a limiting condition i) where the solution to Eq. (C1) writes f (E, x) = f eq + S(E, x) which holds when collisional processes are fast compared to the Hubble rate, thus annihilating the r.h.s. of Eq. (C1).