Improved BBN constraints on Heavy Neutral Leptons

We constrain the lifetime of thermally produced Heavy Neutral Leptons (HNLs) from primordial nucleosynthesis. We show that even a small fraction of mesons present in the primordial plasma leads to the over-production of the primordial helium. This puts an upper bound on the lifetime of HNLs $\tau_{N}<0.02$ s for masses $m_{N}>m_{\pi}$ (as compared to 0.1 s reported previously). In combination with accelerator searches, this allows us to put a new lower bound on the HNLs masses and defining the"bottom line"for HNL searches at the future Intensity Frontier experiments.

We constrain the lifetime of thermally produced Heavy Neutral Leptons (HNLs) from Big Bang Nucleosynthesis. We show that even a small fraction of mesons present in the primeval plasma leads to the over-production of the primordial helium-4. This constrains the lifetime of HNLs to be τN < 0.02 sec for masses above the mass of pion (as compared to 0.1 sec reported previously). In combination with accelerator searches, this allows us to put a new lower bound on the HNLs masses and to define the "bottom line" for HNL searches at the future Intensity Frontier experiments.
Introduction. Heavy neutral leptons (HNLs or righthanded neutrinos) are hypothetical particles capable of explaining neutrino masses and oscillations [1] and resolving other beyond-the-Standard-Model phenomena: the origin of the baryon asymmetry of the Universe [see e.g. 2] and the nature of dark matter [3]. Extending the Standard Model with exactly three HNLs with masses below the electroweak scale and not other heavy physics allows explaining all three beyond-the-Standard-Model phenomena [4][5][6]. An attractive feature of the models with MeV-GeV scale HNLs is that they can be tested by a combination of accelerator searches at the LHC ([see e.g. 1, 7] and refs. therein) and future intensity frontier experiments [see e.g. 8].
The single requirement that HNLs contribute sizeably to the masses of neutrinos combined with existing accelerator exclusions does not allow putting a lower bound on the HNL mass. Such a lower bound can be obtained if we add to the picture constraints from Big Bang nucleosynthesis (BBN). Primordial abundances of 4 He and D are measured with high accuracy [9][10][11][12][13][14]. Their agreement with the Standard Model predictions [15][16][17][18][19] serves as one of the "pillars" of modern cosmology.
The success of the Standard Model-based BBN (SBBN) predictions permits using it to constrain hypothetical particles with lifetimes as small as O(10 −2 ) s ([see e.g. 20] for review). In particular, decays of HNLs in MeV-temperature plasma affect two observable quantities: (i) the abundances of light elements (in this paper we consider only 4 He); (ii) the effective number of relativistic species, N eff . This was used to put an upper bound on HNL lifetime [21][22][23][24][25][26][27][28][29][30][31][32]. These bounds were obtained using decays of HNLs into electromagnetic particles and neutrinos.
The goal of this paper is to derive BBN bounds for HNLs that can decay into mesons. The effect of mesons on the MeV plasma is qualitatively different as they interact with protons and neutrons via strong interactions. Although lifetimes of mesons, τ meson ∼ 10 −8 s, are orders of magnitude smaller than any relevant BBN time scales, they can be present in the plasma as long as HNLs are still abundant and decay. p → n and n → p reactions driven by mesons have a cross-section O(10 16 ) larger than that of weak processes that are normally responsible for the p ↔ n conversion in the absence of mesons (e.g. in Standard Model BBN). Moreover, these reactions, in both directions, have no threshold and roughly equal cross-sections (due to isotopic symmetry).
Therefore, the presence of even a small fraction of mesons in the plasma quickly equilibrates the number densities of protons and neutrons. Once HNLs decay, mesons also disappear (instantaneously as compared to the time scales relevant for BBN), and the neutron-to-proton ratio n n /n p relaxes solely due to the SM (weak) processes. If the HNL's lifetime is short enough, n n /n p relaxes to its SM value before weak reactions freeze out, leaving no observable effect. However, if HNLs (and mesons) survive until T 1.5 MeV and below, there is not enough time to completely relax down to the SBBN value. This residual effect leads to a strong upper bound on the HNL lifetime.
In this work, we demonstrate for the first time that the meson-driven effect strengthens the upper bound on the HNL lifetime by a factor ∼ 5 (down to 0.02 s) as compared to the previous work [22]. In the context of dark scalars, the meson driven p ↔ n conversion and its influence on BBN were studied in [33][34][35][36][37][38].
For temperatures corresponding to such short lifetimes, all Standard Model particles are in thermal equilibrium. This makes all other effects of HNLs on BBN irrelevant and allows deriving the bounds purely analytically, avoiding any computations of complicated Boltzmann equations.
Meson-driven p ↔ n conversion. Sufficiently heavy HNLs can decay into mesons h = π, K, etc (see [39,40] or Appendix C). Charged pions drive the p ↔ n conversion via [33] π − + p → n + π 0 /γ, π + + n → p + π 0 . (1) The cross-section of these reactions is very large: present in the plasma in the amounts at least comparable with that of baryons, they drive the number densities of protons and neutrons to equal values, n n /n p σ π p→n v / σ π n→p v 1. 1 The effect of kaons is qualitatively similar, but leads to a slightly different neutronto-proton ratio (Appendix D 2).
The impact of this effect on primordial 4 He abundance depends on how long mesons remain present in plasma in significant amounts. Once mesons are created, they can (i) scatter and lose energy; (ii) decay; (iii) participate in p ↔ n conversion. The corresponding rates are very different: at MeV temperatures and below, Γ h scat Γ h decay Γ h p↔n (see [37]). The instantaneous number density of mesons is an interplay between their production (via decays of HNLs) and their decays: Here, Br N →h is the branching of HNLs into mesons (Appendix C). n N (T ) is the number density of HNLs. We consider here HNLs that were produced thermally and decouple at some temperature T dec (Appendix B). 2 1 For each of the processes (1), there are no inverse reactions. Indeed, π 0 decays very fast, whereas γs quickly lose their energy. Therefore, the conversion (1) is highly non-equilibrium, and the corresponding value of nn/np is not given by the usual Boltzmann exponent. 2 Besides thermal production, out-of-equilibrium production Therefore, where n dec N is the HNL number density at decoupling, and a(T ) (a dec ) is the scale factor at temperature T (correspondingly, at HNL decoupling).
The number of p ↔ n reactions per nucleon occurring after time t τ N (or below some corresponding temperature T (t)) is thus where n B is the baryon number density, the sum goes over meson species and P conv is the probability for a single meson to interact with nucleons before decaying: At O(1 MeV) temperatures, P conv ∼ 10 −2 − 10 −1 , see Appendix D 1.
mechanisms of HNLs exist [see e.g. 41,42], that we leave for future works (see, however, [31] where a part of this parameter space was explored).  The gray region is excluded as a result of this work (for masses below pion threshold we use the results of [43]). The magenta shaded region corresponds to the domain excluded in [31]. We considered HNLs that are produced thermally and sufficiently short-lived, such that they do not survive until the onset of nuclear reactions (Appendix E). Right panel : comparison of the results of this work (thick blue line) with the results of the previous works [22,25,32] (purple lines) assuming mixing with electron flavor only. Notice that other works have adopted different values for the maximally admissible 4 He abundance when deriving their bounds: Yp,max = 0.2696 in [22,25] and Yp,max = 0.253 in [32] as compared to Yp,max = 0.2573 in this work (see text for details).
The meson driven conversion keeps the value n n /n p 1 roughly until a temperature T 0 when the number of reactions drops below one, and weak SBBN reactions start to relax the n/p ratio down to its SBBN value, see Fig. 1 (left panel). However, if T 0 is close enough to the freeze-out of weak p ↔ n processes, occurring at T n 0.8 MeV, the relaxation is not complete (Fig. 1, right panel). This leads to a positive correction ∆(n n /n p ) as compared to the SBBN case, which translates to an increase of the 4 He abundance ∆Y p .
In this way, the upper bound on the 4 He abundance Y p,max is translated to the lower bound T 0 ≥ T min 0 . Together with the relations (5)-(7), this allow us to find an upper limit on the HNL lifetime τ N : Here, n γ is the number density of photons, η B is the baryon-to-photon ratio, and t(T ) is time-temperature relation. t(T ) is given by the Standard Model relation: Let us rewrite the logarithmic factor in (8) as 3 This is indeed the case for short-lived HNLs with τ N 0.1 s.
HNLs with m N m π and lifetimes τ N 0.1 s decouple while being ultrarelativistic, T dec m N (see Appendix B) implying n N,dec /n γ (T dec ) ≈ 3/2. In SBBN at temperatures T 1 MeV, all particles are at local equilibrium, which define the dynamics of the scale factor: Decays of heavy HNLs violate thermal equilibrium at O(1 MeV) and the scaling (10) is not valid. This leads to an additional decrease of this ratio by a factor of 0.1−0.6 for HNL masses m π m N 3 GeV (we will use 1 3 for normalization below), see Appendix B.
This results in (11) Using values of Br N →h , P conv and the scale factors ratio (Appendices C, D 1, B correspondingly), we conclude that the logarithm term in (11) is O(1) for HNLs in the mass range m N = O(1 GeV) and affects the overall bound very weakly. Therefore, the bound depends only on T min 0 . The presence of mesons increases the 4 He abundance. Therefore, to fix T min 0 (m N ), we need to adopt an upper bound on the primordial 4 He abundance, Y p,max , that is consistent with measurements [44]. The smallest error bars come from measuring Y p in low-metallicity interstellar regions and extrapolating its value to zero metallicity (pioneered in [45]). Several groups [10-14] have determined Y p using this method, albeit with different data and assumptions. The resulting scatter between results is larger than the reported error bars. We treat this difference as an additional systematic uncertainty and adopt the maximal value Y p,max = 0.2573 (see Appendix A). The maximally allowed relative deviation is therefore To relate ∆Y p and T min 0 , we study how the n n /n p ratio is relaxed below T 0 . The relaxation occurs solely via the SBBN reaction, HNLs on Y p . At temperatures T T 0 , for HNLs with lifetimes τ N 0.02 s, all other quantities that are relevant for BBN dynamics -η B , time-temperature relation, the nuclear reactions chain -remain the same as in SBBN, which is because most of HNLs are no longer left in the plasma at these temperatures (see also Appendix D 2). As a result, a value of where T bbn ≈ 84 keV is the temperature of the onset of nuclear reactions in SBBN [15]. The maximal admissible correction (12) is reached for T min 0 = 1.50 MeV, almost independently on HNL mass (see Fig. 1 and Appendix D 2). Plugging T min 0 = 1.50 MeV into (11), we get our final limit To obtain the bound (11), we considered exclusively meson-driven p ↔ n processes for T > T min 0 and only weak SBBN processes for T < T min 0 . We also solved numerically the equation (13) for the neutron abundance in the presence of both mesons-driven and SBBN p ↔ n conversion rates in Appendix D 2, and obtained constraints at the level of 0.019 − 0.021 s, in perfect agreement with the bound (15). We have also repeated our analysis for the case of the GeV-mass scalar that mixes with the Higgs and found an excellent agreement with [34,35].
Our analysis remains valid until τ N reaches O(40 s). HNLs with longer lifetimes survive until the onset of the nuclear reactions. Mesons from their decays may then dissociate already formed nuclei, driving Y p down again, see Appendix E. This effect has been analyzed in [46] demonstrating that such long-lived HNLs are also excluded. Conclusion. We demonstrated that HNLs with semileptonic decay channels significantly affect the primordial 4 He abundance, as mesons from their decays drive the p ↔ n conversion rates away from their SBBN values (c.f. [33][34][35]). In order to avoid 4 He overproduction, mesons should disappear from the primordial plasma by T = T min 0 1.50 MeV. The neutron abundance will then have enough time to relax down to its SBBN value before the onset of deuteron formation. These requirements severely constrain the parameter space of the HNLs with 0.023 s ≤ τ N ≤ 40 s for masses m N > 140 MeV. Accidentally, HNLs with τ N 40 s are also excluded [46,47], but we do not consider this case here.
We show our results for the case of two nearlydegenerate in mass HNLs that entered in thermal equilibrium and then froze out (which is reflected in both abundance calculation and decay width/patterns), as motivated by the Neutrino Minimal Standard Model (or νMSM) [see e.g. 6, 48,49]). The final bounds for different mixing patterns are shown in Figs. 2 and 3. Our constraints can be generalized to other HNL models, see e.g. [46].
Confronted with the bounds from accelerator searches, we ruled out HNLs with mass below 500 MeV (for electron mixing) and 350 MeV (for muon mixing). Moreover, tighter bound means that future searches at Intensity Frontier (specifically, SHiP experiment [50]) can reach the BBN bottom line and completely rule out HNLs with the masses up to 750 MeV, which was not the case before [see e.g. 8, 51].
The comparison with the previous results [22,25,32] is shown in Fig. 2 (right panel). Our bound (15) is a factor of ∼ 5 stronger than the previous result [22]. The recent reanalysis [32] did not take into account the effects of mesons, therefore their results are a factor 2 − 3 less conservative.
The clear qualitative effect discussed in this paper not only leads to a tighter bound on HNL lifetime and provides an reachable goal for experimental searches, but also allows for an analytic description, unusual in the realm of BBN predictions driven by sophisticated numerical codes. U 2 U e 2 :U μ 2 :U τ 2 = 0:0:1 Figure 3. Bounds for HNLs mixed with a particular flavor. The blue area is excluded by our present analysis combined with [43] (for HNL masses below the charged pion production threshold Over the last 6 years five works determined the primordial 4 He abundance from stellar measurements [10-14]. The formal statistical errors of Y p are at the level of 1 − 3%, however, the scatter between different groups is larger, see Fig. 4. All these works determine astrophysical Helium abundance through measurements of recombination emission lines of 4 He and H in the metal-poor extragalactic ionized regions, then linearly extrapolating the measurements to zero metallicity. Given the high precision of the results, it is important to take into account various smaller effects: including 4 He fluorescent emission, different ion temperatures, spatial temperature fluctuations, and others [58,59]. Additionally, while it is true that the metallicity and Helium abundance are positively correlated, the linear extrapolation to zero-metallicity may be prone to systematic uncertainties. The value of Y p predicted within the framework of SBBN isȲ p = 0.24709 ± 0.00019 (see, e.g., [15]). The effect of mesons leads to an increase of Y p as compared to the SBBN value. Therefore, in order to get a conservative upper bound we assume that the maximally allowed Y p is given by the 1σ deviation from the maximal value predicted by [10-14], which is Y p,max = 0.2573. Note that this upper value significantly deviates from the PDG-recommended value [44] [43]. Dilution factor is calculated, when most of HNLs has decayed and do not contribute to entropy density. Note, that we define abundance at the moment of decoupling, hence it does not change with decays.
the O(1) accuracy, the decoupling temperature T dec may be estimated as where T ν,dec ≈ 1.4 MeV is the decoupling temperature of active neutrinos, a parameter n int = where n N,eq is the number density of HNLs at equilibrium (i.e., calculated using the Fermi-Dirac distribution). Expressing the total HNL number density in terms of the abundance Y N , where s = g * After the freeze-out, the comoving number density of HNLs changes only due to HNL decays. The physical number density thus evolves as Br N→h N→π -X, e mixing N→π -X, μ mixing N→π -X, τ mixing N→K -X, e mixing N→K -X, μ mixing N→K -X, τ mixing N→K 0 L X, e mixing N→K 0 L X, μ mixing N→K 0 L X, τ mixing where a −1 sbbn (T ) ∝ g 1/3 * T is the scale factor in SBBN, and the scale factors are evaluated at times t τ N . To calculate ζ, we solve the Friedmann equation under an assumption that neutrinos are in perfect equilibrium and neglecting the mass of electrons: where the number density of HNLs is given by Eq. (B5). This is a reasonable assumption, since most of the HNLs with lifetimes τ N 0.1 s decay much earlier than neutrinos decouple.
3. Strong meson effects mean that we need to trace the number of HNLs even at times t τ N : where t(T ) is the same as in SBBN. 5 Because of the suppression, the effect of this population on the expansion of the Universe may be neglected. However, this exponential tail still may produce mesons in amounts sufficient to change the dynamics of the n/p ratio.
The values of the HNL abundance and the dilution factor versus its mass and lifetime are given in Fig. 5.

Appendix C: Hadronic decays of HNLs
In this work, we consider a pair of HNLs, degenerate in mass and having similar mixing angles. Two such HNLs form a single quasi-Dirac fermion [60,61]. The abundance of a meson h produced from such HNLs is proportional to the quantity Y N · Br N →h . The mass dependence of Br N →h for different mesons h and mixing patterns is shown in Fig. 6. We are interested only in the abundances of light mesons (pions and kaons) and for HNL masses well above pion/kaon thresholds we should account for "secondary mesons". This is discussed below, mainly following [40].
Decays into pions. In the case of the pure e/µ mixings, the charged pion production threshold corresponds to m N = m π + m l , where l = e/µ. For τ mixing, the similar charged current-mediated channel opens up only at m N = m τ + m π 1.9 GeV. However, for all types of mixings charged pions may appear as secondary particles in decays of neutral mesons, Therefore, for τ mixing charged pions may appear at masses m N ≥ m η 0 . We use the branching ratios Br η 0 →π ± X ≈ 0.27, Br ρ 0,± →π ± X ≈ 1 [44]. Above m N 1 GeV, decays of HNLs into pions cannot be approximated by single meson decays. Indeed, decays of GeV mass range HNLs are similar to decays of τ lepton [40], whereas for the latter hadronic decays are dominated by multi-pion channels [44]. We estimate the width of multi-pion decays as the difference between the total width into quarks and the width into single mesons: For multiplicities N of decays of HNLs into charged pions N → π ± (i.e., the amount of π ± per multi-hadronic decay of HNLs), we will use multiplicities for multihadronic decays of τ leptons. Namely, The effective branching into π − from multi-pion decays is Since the bound on the meson driven p ↔ n conversion is only logarithmically sensitive to the value of Br N →π ± , our results depend on these assumptions weakly.
Decays into kaons. Below m N = m φ , charged kaons may appear only through the mixing with e/µ in the process N → K − l. This decay is Cabibbo suppressed [40] and almost two orders of magnitude smaller than into pions. Neutral kaons appear only in the final states with three or more particles (such as N → K 0 +K 0 + ν α and N → K + +K 0 + − , etc). HNLs heavier than φ meson may produce both charged and neutral kaons via decays N → φν, φ → KK. We assume that K 0 contains equal admixtures of K 0 L and K 0 S , i.e. Br N →K 0 L = Br N →K 0 /2. We use the branching ratios Br φ→K − ≈ 0.5, Br φ→K 0 L ≈ 0.34 [44].
Appendix D: Changes in p ↔ n rates due to the presence of mesons In this section, we provide details on our estimate of the effect of mesons on BBN.

Processes and cross-sections
.1. Pions. The threshold-less processes with charged pions are The cross-sections at threshold are [37] σ π − p→n v ≈ 4.3 · 10 −23 F π c (T ) m 3 /s, where F h c is the Sommerfeld enhancement of the cross-section due to presence of two oppositely charged particles in the in-state: where v e ≈ T m h + T mp is the relative velocity between a nucleon and a meson. F c is of order of one at T 1 MeV.
.2. Kaons. The threshold-less n ↔ p conversions driven by kaons are K − + p → Σ ±/0 /Λ + π ∓/0 /π 0 → n + 2π, where Λ, Σ are the lightest strange hadronic resonances [33]. Their effect is similar to the one of pions, but with small differences: (i) cross-sections of above reactions are higher than the cross-sections of (1) 6 , (ii) there is no isotopic symmetry -K + mesons do not contribute to p ↔ n conversion, since there are no threshold-less processes n + K + → p + X. Indeed, the process n + K + → p + K 0 has the threshold Q ≈ 2.8 MeV, while the threshold-less processes going through s-quark resonances, similar to (D4), would require resonances with negative strangeness and positive baryon number, that do not exist, (iii) neutral kaons do not lose the energy before decaying (however, we follow [33] and approximate the cross-sections by threshold values).
The threshold cross-sections are .3. Conversion probabilities. A probability for a meson h to convert p ↔ n before decaying is given by where Γ h decay is the decay width and n B is the baryon number density. The decay widths of mesons are [44] Γ π ± decay ≈ 3.8 · 10 7 s −1 , Γ K − decay ≈ 8.3 · 10 7 s −1 , Γ K 0 L decay ≈ 2 · 10 7 s −1 (D8) Using (D2), (D5), (D8), for the p → n conversion probabilities we obtain The largeness of the probabilities is caused by the fact that the decay of mesons proceeds through weak interactions, while the p ↔ n conversion is mediated by strong interactions. In particular, at T 2 MeV kaons participate in the conversion faster than they decay.

Numeric study
To verify the analytic estimate (15), we numerically solve equation for the neutron abundance X n , where we include both weak conversion p ↔ n processes and the meson driven processes (D1)-(D4). The system of equations has the form The reason is that these reactions have higher available phase space and go through hadronic resonances.
Here the quantities dX n dt π = (1 − X n )n π − σ π − p→n v − X n n π + σ π + n→p v , are the rates of change of X n due to different mesons (K = K − /K 0 L ); n B is the baryon number density n B = η B n γ . In equations for the number density of mesons n h , the first term comes from HNLs, the second due to decays of mesons and the last term is due to p ↔ n conversion. The time-temperature relation and the scale factor dynamics are provided by the solution of Eq. (B7), and the HNL number density may be obtained using Eq. (B5).
During times t eq (Γ h decay ) −1 ∼ 10 −8 s, which are small in comparison to any other time scale in the system, the solution for n h reaches the dynamical equilibrium: where K = K − /K 0 L . Therefore, we solve a single equation where we use meson number densities given by Eqs. (D12) and (D13) in the meson-driven conversion rates (D11). The results are shown in Fig. 7. Our main result is the right panel of Fig. 7 -it shows that the value T min 0 1.50 MeV and that its variation as a function of the HNL mass is within ±1%.  With the help of Eqs. (D2), (D5), we obtain the value of the neutron abundance driven solely by a given meson h. As long as T T 0 (see Eq. (8) and left panel of Fig. 7), the weak interaction processes may be completely neglected, and the resulting X n are given by X π ± n = σ π − p→n v · n π − σ π − p→n v · n π − + σ π + n→p v · n π + ≈ 0.9F π c (T ) 1 + 0.9F π c (T ) The values of X π − /K − n grow with the decrease of the temperature due to the growth of the Coulomb factor F c , which enhances the rate of the p → n process.
The quantities (D15) provide us the qualitative estimate of the value of X n in presence of different mesons, Fig. 7. Below the kaon production threshold, X h n = X π ± n . At larger masses, in order to find X h n we need to set the whole right hand-side of Eq. (D10) to zero. Below the K 0 L production threshold (which occurs at m N = m φ ), the value of X h n grows, since charged kaons tend X n to higher values than X π − n . Above the neutral kaon production threshold, the ratio Br N →K − /Br N →π − increases (Fig. 6) and X h n grows further. However, kaons K 0 L , that are present in small amounts, somewhat diminish this growth.
The value of X h n (m N ) provide us the mass dependence of T min 0 (m N ), which is the smallest temperature allowed by observations (c.f. Fig. 1). We show it in Fig. 7 (right panel).
3. SBBN evolution at T < T min 0 If HNLs disappear from the plasma before neutrinos froze out, the evolution of the neutron abundance and subsequent nuclear reactions proceed exactly as in SBBN case (albeit with modified initial value of X n at T = T min 0 ). Indeed, the onset of nuclear reactions is determined by the dynamical balance between reactions of deuterium synthesis and dissociation. This balance depends on the value of η B . The latter gets diluted by the factor ζ due to decays of HNLs, see Section B. However, we fix η B at the beginning of nuclear reactions to be the same as measured by CMB. This of course means that η B has been ζ −1 times higher before decays of HNLs, but no observables can probe the value of η B in this epoch.
Another ingredient that affects dynamics of nuclear reactions is the time-temperature relation, traditionally encoded in the value of N eff . If HNLs have τ N 0.02 s, neutrinos are in equilibrium and therefore HNL decays do not change N eff , see detailed analysis in [62].
As a result, evolution of primordial plasma below T min