Improved big bang nucleosynthesis constraints on heavy neutral leptons

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 overproduction of primordial helium-4. This constrains the lifetime of HNLs to be τ N < 0 . 02 s for masses above the mass of the pion (as compared to 0.1 s reported previously). In combination with accelerator searches, this allows us to put a new lower bound on the HNL masses and define the “ bottom line ” for HNL searches at the future Intensity Frontier experiments.


I. INTRODUCTION
Heavy neutral leptons (HNLs or right-handed neutrinos) are hypothetical particles capable of explaining neutrino masses and oscillations [1] and resolving other beyond-the-Standard-Model phenomena, i.e., the origin of the baryon asymmetry of the Universe (see, e.g., Ref. [2]) and the nature of dark matter [3]. Extending the Standard Model (SM) with exactly three HNLs with masses below the electroweak scale and no other heavy physics allows to explain 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., Refs. [1,7] and references therein) and future intensity frontier experiments (see, e.g., Ref. [8]).
The single requirement that HNLs contribute sizably to the masses of neutrinos combined with existing accelerator exclusions does not allow us to put 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., Ref. [20] for a 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), and (ii) the effective number of relativistic species N eff . This was used to put an upper bound on the 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 the 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 Þ times 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 Ref. [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 us to derive the bounds purely analytically, avoiding any computations of complicated Boltzmann equations.

II. MESON-DRIVEN p ↔ n CONVERSION
Sufficiently heavy HNLs can decay into mesons h ¼ π, K, etc. (see Refs. [39,40] or Appendix C). Charged pions drive the p ↔ n conversion via [33] The cross section of these reactions is very large: The large cross section, absence of a threshold, and isotopic symmetry of these processes mean that if pions are present in the plasma in amounts at least comparable to that of baryons, they drive the number densities of protons and neutrons to equal values, n n =n p ≃ hσ π p→n vi=hσ π n→p vi ≃ 1. 1 The effect of kaons is qualitatively similar, but leads to a slightly different neutron-to-proton ratio (Appendix D 2).
The impact of this effect on the 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, or (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 Ref. [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 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.
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 baryonto-photon ratio, and tðTÞ is the time-temperature relation. tðTÞ is given by the Standard Model relation tðTÞ is the reduced Planck mass, and g Ã ðTÞ ≃ 10.6 for T ≃ 1-2 MeV. 3 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 nonequilibrium, and the corresponding value of n n =n p is not given by the usual Boltzmann exponent. 2 Besides thermal production, out-of-equilibrium production mechanisms of HNLs exist (see, e.g., Refs. [41,42]), which we leave for future works (see, however, Ref. [31] where a part of this parameter space was explored). 3 This is indeed the case for short-lived HNLs with τ N ≪ 0.1 s. Let us rewrite the logarithmic factor in Eq. (8) as a dec a 0 3 n N;dec 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: a dec T dec a 0 T min 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 Using values of Br N→h , P conv , and the scale factors ratio (Appendixes C, D 1, and B, respectively), we conclude that the logarithm term in Eq. (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 [43]. The smallest error bars come from measuring Y p in low-metallicity interstellar regions and extrapolating its value to zero metallicity (pioneered in Ref. [44]). Several groups [10][11][12][13][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, albeit with the altered initial condition p↔n ðtÞ are SBBN rates; see Ref. [15].] A non-SBBN value of X n ðT 0 Þ is the dominant effect of short-lived 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 X n ðT 0 Þ is translated into ΔY p via 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 of HNL mass (see Fig. 1 and Appendix D 2). Plugging T min 0 ¼ 1.50 MeV into Eq. (11), we get our final limit To obtain the bound (11), we consider exclusively mesondriven p ↔ n processes for T > T min 0 and only weak SBBN processes for T < T min 0 . We also numerically solve Eq. (13) for the neutron abundance in the presence of both mesondriven and SBBN p ↔ n conversion rates in Appendix D 2, and obtain constraints at the level of 0.019-0.021 s, in perfect agreement with the bound (15). We also repeat our analysis for the case of the GeV-mass scalar that mixes with the Higgs and find an excellent agreement with Refs. [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 Ref. [45], demonstrating that such long-lived HNLs are also excluded.

III. 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 (cf. Refs. [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 [45,46], but we did not consider this case here.
We show our results for the case of two nearly-degenerate 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., Refs. [6,47,48]). 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., Ref. [45].
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, a tighter bound means that future searches at Intensity Frontier (specifically, the SHiP experiment [51]) can reach the BBN bottom line and completely rule out HNLs with masses up to 750 MeV, which was not the case before (see, e.g., Refs. [8,57]).
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, and therefore their results are a factor of 2-3 less conservative. The gray region is excluded as a result of this work (for masses below the pion threshold we use the results of Ref. [49]). The magenta shaded region corresponds to the domain excluded in Ref. [31]. We consider 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 adopted different values for the maximally admissible 4 He abundance when deriving their bounds: Y p;max ¼ 0.2696 in Refs. [22,25] and Y p;max ¼ 0.253 in Ref. [32], as compared to Y p;max ¼ 0.2573 in this work (see text for details).
The clear qualitative effect discussed in this paper not only leads to a tighter bound on the HNL lifetime and provides a reachable goal for experimental searches, but also allows for an analytic description, unusual in the realm of BBN predictions driven by sophisticated numerical codes.

ACKNOWLEDGMENTS
We thank K. Bondarenko, N. Sabti for useful discussions, and Yu. Izotov for discussions related to the uncertainties of primordial helium determination. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (GA 694896) and from the Carlsberg Foundation.
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 of these works determined the 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 AE 0.00019 (see, e.g., Ref. [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 in Refs. [10][11][12][13][14], which is Y p;max ¼ 0. In this work we consider HNLs that have entered thermal equilibrium in the early Universe and then froze out. Below we provide necessary details of their evolution, to make our FIG. 4. Measurements of Y p from recent works [10][11][12][13][14]. The green shaded region is the PDG recommended value [43] (with AE1σ). The gray dashed line denotes the SBBN predictionȲ p ¼ 0.247 from Ref. [15]. The red dash-dotted line is the maximal admissible value Y p;max on which we base our analysis.
FIG. 3. Bounds for HNLs mixed with a particular flavor. The blue area is excluded by our present analysis combined with Ref. [49] (for HNL masses below the charged pion production threshold). The dark gray area denotes the excluded HNL parameter space from previous searches [1], including the latest NA62 search [50]. The red and greed dashed lines show the sensitivity of several future Intensity Frontier experiments with the highest sensitivity in the regions of interest: SHiP [51,52] and DUNE [53][54][55] (see Ref. [8]). Finally, the black dashed line denotes the seesaw bound applicable if two degenerate-in-mass HNLs are responsible for neutrino oscillations (as in the νMSM) [1,56]. The lower bound of the blue region shows the limit of the applicability of our approach: we only consider HNLs that are produced thermally and are sufficiently short lived such that they do not change the nuclear reaction framework by their meson decay products. presentation more self-contained [49]. The evolution of HNLs with masses above the pion production threshold proceeds in three stages: (1) The freeze-out of HNLs occurs at some temperature T dec defined via where Γ int N is a rate of processes that keep HNLs in thermal equilibrium with the primordial plasma, 4 and the Hubble rate H corresponds to its Standard Model value (at the decoupling, HNLs with masses m N ≳ m π and lifetimes in the range of interest contribute just a small fraction of the energy density of the Universe). Within 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, and the parameter n int ¼ MeVÞ to ≃10 at GeV temperatures (for Dirac HNLs). According to this estimate, for lifetimes τ N ∼ 0.01-0.1 s, HNLs decouple while being nonrelativistic for m N ≲ 200 MeV and relativistic otherwise. To improve this estimate, we solve the following equation: 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 Ã 2π 2 45 T 3 is the entropy density, we find Y N ≃ 0.6 g Ã ðT − Þ in the ultrarelativistic regime and a factor of Oð2Þ larger otherwise; see Fig. 5.
After the freeze-out, the comoving number density of HNLs changes only due to HNL decays. The physical number density thus evolves as (2) Decays of HNLs inject energy into the primordial plasma. This changes the time-temperature relation and the scale factor evolution as compared to SBBN. The HNL decays provide additional dilution of any decoupled relics (including themselves) in comparison to the SBBN case: 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 the assumption that neutrinos are in perfect equilibrium and we neglect 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 may still 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 Ref. [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-currentmediated 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 →π AE X ≈ 0.27, Br ρ 0;AE →π AE X ≈ 1 [43]. 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 [43]. 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 → π AE (i.e., the amount of π AE 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 Br multi-pion Since the bound on the meson driven p ↔ n conversion is only logarithmically sensitive to the value of Br N→π AE , 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 2 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 þ l − , etc.).

APPENDIX D: CHANGES IN p ↔ n RATES DUE TO THE PRESENCE OF MESONS
In this Appendix, we provide details on our estimate of the effect of mesons on BBN.

Processes and cross sections a. Pions
The thresholdless processes with charged pions are The cross sections at threshold are [37] hσ π − p→n vi ≈ 4.3 × 10 −23 F π c ðTÞ m 3 =s; where F h c is the Sommerfeld enhancement of the cross section due to the presence of two oppositely charged particles in the in-state: where v e ≈ ffiffiffiffi ffi is the relative velocity between a nucleon and a meson. F c is of order of one at T ≃ 1 MeV.

b. Kaons
The thresholdless n ↔ p conversions driven by kaons are where Λ, Σ are the lightest strange hadronic resonances [33]. Their effect is similar to that of pions, but with small differences. (i) The cross sections of the 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 thresholdless processes n þ K þ → p þ X. Indeed, the process n þ K þ → p þ K 0 has the threshold Q ≈ 2.8 MeV, while the thresholdless processes going through s-quark resonances, similar to Eq. (D4), would require resonances with negative strangeness and positive baryon number, which do not exist. (iii) Neutral kaons do not lose energy before decaying (however, we follow Ref. [33] and approximate the cross sections by threshold values).
The threshold cross sections are The reason is that these reactions have higher available phase space and go through hadronic resonances.
c. 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 [43] Γ π AE decay ≈ 3.8 × 10 7 s −1 ; Using Eqs. (D2), (D5), and (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 the 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 8 > > > > > > > > > > > > < > > > > > > > > > > > > : Here the quantities dX n dt π ¼ ð1 − X n Þn π − hσ π − p→n vi − X n n π þ hσ π þ n→p vi and dX n dt 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 is 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 the 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, which shows that the value T min 0 ≃ 1.50 MeV and that its variation as a function of the HNL mass is within AE1%.
With the help of Eqs. (D2) and (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 the left panel of Fig. 7], the weak interaction processes may be completely neglected, and the resulting X n are given by 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 with the qualitative estimate of the value of X n in the presence of different mesons (Fig. 7). Below the kaon production threshold, X h n ¼ X π AE 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 push 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 , which are present in small amounts, somewhat diminish this growth.
The value of X h n ðm N Þ provides us with the mass dependence of T min 0 ðm N Þ, which is the smallest temperature allowed by observations (cf. Fig. 1). We show it in Fig. 7 (right panel).

SBBN evolution at T < T min 0
If HNLs disappear from the plasma before neutrinos freeze out, the evolution of the neutron abundance and subsequent nuclear reactions proceed exactly as in the SBBN case (albeit with a 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 Appendix B. However, we fix η B to be the same as that measured by cosmic microwave background observations by the beginning of nuclear reactions. 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 the 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 the detailed analysis in Ref. [62].
As a result, the evolution of primordial plasma below T min 0 is governed by the SBBN equations.

APPENDIX E: DISSOCIATION OF LIGHT ELEMENTS BY MESONS AND DOMAIN OF VALIDITY OF OUR TREATMENT
The analysis presented in this paper is valid as long as HNLs decay before neutrino decoupling times. This puts an upper limit on the HNL lifetime [Eq. (15)]. It is clear, however, that once the lifetime is long enough and HNLs (their decay products) can survive until the onset of nuclear reactions, our treatment needs to be changed.
Indeed, pions, if present in the plasma at temperatures T ≲ T BBN , dissociate light nuclei. The 4 He thresholdless dissociation processes are (see Ref. [33]) π − þ 4 He → T þ n; π − þ 4 He → D þ 2n; To estimate the lifetimes at which the processes (E1) can be neglected, we compare the number density of mesons available for the dissociation with the number density of 4 He nuclei: n h He diss ðT BBN Þ ≪ n He ðT BBN Þ; ðE2Þ cf. Eq. (7). Here, n π He diss is defined via n π He diss ðT BBN Þ ¼ n N · Br N→π − · P He diss ðE3Þ and P4 He diss is the probability for a single meson to dissociate 4 He nuclei before decaying: where we used the total cross section of the dissociation processes (E1), hσ π He diss vi ≃ 6.5 · F Heπ − mb [a factor F Heπ − ≃ 3.5 accounts for the Coulomb attraction (D3)].
In our estimates we use T BBN ¼ 84 keV, assuming that all free nucleons become bounded in 4 He nuclei at this temperature. We also do not take into account that after the dissociation of helium the abundance of lighter elements will also be increased significantly. Assuming n He ≃ n B =4 in Eq. (E2), and using Eqs. (E3) and (E4) with the HNL number density given by Eq. (B8), we arrive at the upper bound on HNL lifetimes for which our analysis is applicable: τ N ≲ 40 s. Notice that the dependence on ζ is only logarithmic, and therefore its exact value and its mass dependence play no role. The analysis for τ N > 40 s should be performed separately; see, e.g., Ref. [45].