MeV-scale seesaw and leptogenesis

We study the type-I seesaw model with three right-handed neutrinos and Majorana masses below the pion mass. In this mass range, the model parameter space is not only strongly constrained by the requirement to explain the light neutrino masses, but also by experimental searches and cosmological considerations. In the existing literature, three disjoint regions of potentially viable parameter space have been identified. In one of them, all heavy neutrinos decay shortly before big bang nucleosynthesis. In the other two regions, one of the heavy neutrinos either decays between BBN and the CMB decoupling or is quasi-stable. We show that previously unaccounted constraints from photodisintegration of nuclei practically rule out all relevant decays that happen between BBN and the CMB decoupling. Quite remarkably, if all heavy neutrinos decay before BBN, the baryon asymmetry of the universe can be quite generically explained by low-scale leptogenesis, i.e. without further tuning in addition to what is needed to avoid experimental and cosmological constraints. This motivates searches for heavy neutrinos in pion decay experiments.


Introduction
The observation of neutrino-flavour oscillations is one of the few hints for new physics beyond the Standard Model that have been discovered to date. In fact, it is the only one found in the laboratory that has been established beyond doubt. One way of explaining these oscillations is by adding right-handed neutrinos ν R to the Standard Model (SM) of particle physics, thus giving mass to the light neutrinos [1][2][3][4][5][6]. Such a type-I seesaw model with right-handed neutrino masses below the electroweak scale is a minimal and testable extension of the SM that can simultaneously explain the generation of the observed neutrino masses as well as the baryon asymmetry of the Universe. With the absolute mass scale of the SM neutrinos being bounded only from above, most studies in this context merely consider two right-handed neutrinos, thus leaving one SM neutrino massless. This drastically reduces the complexity of the problem, and in many cases serves as a good proxy for the relevant dynamics. However, based on both theoretical and experimental considerations, it is necessary to go beyond this simplification. Firstly, all other fermions in the SM come in three generations and overarching concepts such as gauging the difference of baryon and lepton number -with a possible embedding in a grand unified theory -mandate the introduction of three generations of right-handed neutrinos. Secondly, the combination of a high-dimensional parameter space together with neutrino oscillation data, constrains the theory to highly non-trivial sub-manifolds of the parameter space, where the naive intuition gained from the simplified model with only two neutrinos may fail. Finally, predicting particles at an energy scale within the reach of collider experiments, a selling point of this model is its falsifiability. To guide future experimental efforts, it is thus mandatory to map out the full range of potential observables, especially since it is well known that the inclusion of three right-handed neutrinos can significantly change the JHEP01(2021)200 experimentally viable parameter space [7], cosmological constraints [8], and the perspectives for leptogenesis [9]. In this work we study a comparably unexplored region of parameter space in which all heavy neutrinos have masses below the pion mass, kinematically limiting their decay products to SM neutrinos, electrons, positrons, and photons.
Suppressing SU (2) indices for brevity, the most general renormalisable Lagrangian including SM fields and the right-handed neutrinos ν R reads Here L and φ are the left-handed lepton and Higgs doublet of the SM, respectively, F is the matrix of Yukawa couplings, and ε denotes the totally antisymmetric SU(2) tensor. The Majorana mass matrix M M introduces a new fundamental scale in nature, which is usually referred to as the seesaw scale. More precisely, for n flavours of ν R , eq. (1.1) contains n new dimensionful parameters that can be identified with the eigenvalues of M M , which roughly coincide with the physical masses M i of the heavy neutrino mass eigenstates N i (see eq. (2.3) below). The phenomenological and cosmological implications of the ν R 's existence strongly depend on the choice of the seesaw scale(s) (cf. e.g. [10] for a review).
A particularly intriguing feature of this model is the fact that the same ν R that give masses to the light neutrinos can also explain the observed matter-antimatter asymmetry in the early universe, which is believed to be the origin of all baryonic matter that is present today. 1 This is realized via the process of leptogenesis [12], which is feasible for a very wide range of possible M i (see ref. [13] for a recent review). For M i above the electroweak scale, the asymmetry is typically generated during the freeze-out and decay of the heavy neutrinos [12] ("freeze-out scenario"), while for M i below the electroweak scale, it is instead generated during their production [14][15][16] ("freeze-in scenario"). 2 It is well-known that leptogenesis is in principle feasible with M i in the range of a few MeV [18]. However, in this mass range, the model parameter space is strongly constrained by laboratory experiments, cosmology, and astrophysics.
Constraints on the properties of heavy neutrinos are conveniently expressed in terms of the mixing angles θ ai , (cf. eq. (2.2) below). In fact, for given M i , the values of θ ai determine the thermal N i production rate in the early universe, the N i lifetime, the N i contribution to the generation of light neutrino masses, and the N i production cross-section in experiments. For masses below ∼ 100 MeV and values of θ ai that are small enough to satisfy exclusion bounds from various laboratory experiments, the heavy neutrinos tend to have lifetimes larger than 0.1 s. This means that their presence in the primordial plasma and their decay may affect cosmological observables, such as the abundances of light elements that are produced during big bang nucleosynthesis (BBN), or the anisotropies in the cosmic microwave background (CMB). The resulting constraints on M i and θ ai JHEP01(2021)200 have e.g. been summarised in [19]. 3 Usually, these limits can be avoided for sufficiently small values of θ ai . However, since the mixing angles also govern the size of the light neutrino masses, there exist additional lower bounds on different combinations of θ ai from the requirement to explain the observed light neutrino oscillation parameters. These lower bounds depend on the number n of right-handed neutrino flavours and the mass m lightest of the lightest neutrino (cf. [22] for a recent discussion). In the minimal model with n = 2 and m lightest = 0, the seesaw mechanism necessarily enforces that all N i reach thermal equilibrium if their masses are below ∼ 100 MeV [23]. In combination with bounds from direct searches, this practically rules out the entire mass range below ∼ 100 MeV (∼ 350 MeV) for normal (inverted) ordering of the light neutrino masses [24]. 4 In the next-to-minimal model with n = 3 considered here, one of the N i -which we may call N 1 without loss of generality -can have small enough mixings θ a1 to avoid equilibration in the early universe and hence is no longer constrained by the lower bound on M i if m lightest 10 −3 eV [8]. 5 This leaves three distinct regions of parameter space for models with n = 3 and all M i below the pion mass. In scenario I) all three N i decay before BBN. A global fit of direct and indirect experimental constraints in this region has recently been performed in [7]. In scenario II) two of the N i decay before BBN. The third one never reaches thermal equilibrium and decays between BBN and the decoupling of the CMB. In scenario III) two of the N i decay before BBN. The third one is quasi-stable and contributes to the Dark Matter [27,28]. This scenario corresponds to the well-known Neutrino Minimal Standard Model (νMSM) [15,29].
In the present work, we present two new results regarding these scenarios. Firstly, we demonstrate that scenario II) is ruled out when combining previously unaccounted constraints from photodisintegration after BBN with constraints from CMB anisotropies [21] and the ionisation of the intergalactic medium [20]. Secondly, we find that the baryon asymmetry generated in scenario I) generically is of the right order of magnitude to explain the observed matter-antimatter asymmetry. This surprising result indicates that within the highly constrained region of parameter space where all experimental constraints are satisfied, no or little additional tuning is needed for successful leptogenesis. These results extend the previous parameter scan of scenario I) in [9] to smaller masses, M i 50 MeV. Finally, let us note that we do not consider baryogenesis in scenario III) and instead refer the reader to [18,30] for a comprehensive overview and to [17,24,[31][32][33][34] for recent updates 3 The authors in [19] ruled out lifetimes longer than the CMB decoupling time by rescaling the CMB bounds on decaying Dark Matter particles found in [20]. These were obtained under the assumption that the particles have a lifetime that exceeds the age of the universe and can therefore strictly speaking not be applied in all of the parameter space considered here. However, it turns out that the parameter region where this rescaling is not applicable is ruled out by the results obtained in [21], so that we can safely apply the bounds presented in [19] here. 4 The authors of [24] assumed a mass degeneracy among the Ni. However, since both, the lifetime bound from BBN and constraints from direct searches in good approximation apply to each Ni individually, this can at most introduce a factor 2 in the upper bound on the mixing (if the two Ni cannot be distinguished kinematically), which will not change these conclusions. 5 We do not consider the small window of M1 in the eV range that was reported in [8] because the scenario of a eV seesaw [25] is meanwhile even more disfavoured by cosmological date [26].

JHEP01(2021)200
on the viable parameter space in this model. The remainder of this article is organised as follows. In section 2, we summarise the existing laboratory and cosmological constraints, before introducing our new bound from photodisintegration after BBN in section 2.3. We comment on the supernova bound in section 2.4, which could rule out the entire scenario I) but comes with some uncertainties. We summarise all constraints in section 2.5, demonstrating that the neutrino oscillation data can be accounted for in the remaining parameter space. Section 3 is dedicated to the study of leptogenesis in scenario I), followed by a brief conclusion in section 4.

Laboratory constraints
The strongest experimental constraints on the heavy neutrino properties come from the requirement to explain the light neutrino oscillation data. If the eigenvalues of the Majorana mass matrix M M are at least a few eV in magnitude, there exist two distinct sets of mass eigenstates after electroweak symmetry breaking, which can be represented by the flavour vectors of Majorana spinors Here c.c. denotes the c-conjugation which e.g. acts as ν c is the standard light neutrino mixing matrix, and U N is its equivalent among the heavy neutrinos. The mixing between left-and right-handed neutrinos is quantified by the entries of the matrix with the Higgs field expectation value v, and the mass matrices for ν and N are given by The requirement to explain the observed light neutrino mass splittings m 2 i − m 2 j as well as the mixing angles in the matrix U ν imposes constraints on the matrix m ν , and therefore on F and M M . At low energies, this leads to restrictions on the relative size of the heavy neutrino mixing with individual SM flavours [7,24,31,[35][36][37][38][39], i.e. on the quantities U 2

JHEP01(2021)200
There also is a lower bound on the different U 2 i from neutrino oscillation data [36,37] which roughly reads U 2 i > m lightest /M i (cf. [22] for a recent discussion). The presence of weak interactions implies that a wide range of experiments is sensitive to the existence of the heavy neutrinos. An updated overview of the existing constraints that we are aware of can be found in [7]. Broadly speaking, one can distinguish between direct and indirect searches. Direct searches are experiments in which the N i appear as real particles. If kinematically allowed, the N i production cross-section is roughly given by σ N i ∼ a U 2 ai σ νa , with σ νa being the production cross-section for a SM neutrino ν a . Hence, direct searches always impose upper bounds on the different U 2 ai . For sub-GeV masses this mainly includes beam dump experiments and peak searches. Indirect searches include precision tests or searches for rare processes in the SM that are indirectly affected by the existence of the heavy neutrinos, e.g. through the modification of the light neutrinos' interactions via the mixing θ. In the mass range considered here, direct searches strongly dominate, 6 in particular from PIENU [41,42], KEK [43], LBL [44], SIN [45], TRIUMF [46] and CHARM [47] (cf. also [48][49][50]). All of these constraints are summarized in the grey regions in figure 1. The only indirect constraint that is relevant in this region comes from neutrinoless double β-decay (0νββ). However, the rate of the 0νββ-decay can be suppressed even for mixing angles that are orders of magnitude larger than the ones considered here if one requires that the Lagrangian in eq. (1.1) approximately conserves a generalisation of the SM lepton number L (more precisely, the difference between baryon number B and L) under which the heavy neutrinos are charged [51,52], 7 and the current bound on the 0νββ-lifetime rules out only a small fraction of the leptogenesis parameter space [9,31,54].

Summary of previously known cosmological constraints
In the following we assume that the N i have masses of at least a few MeV. For lower M i the analysis would have to take into account the fact that the N i are still partially relativistic when the SM neutrinos decouple from the primordial plasma, which would require a more careful analysis of BBN. This will generally lead to even stronger bounds than the ones considered here because the N i would e.g. contribute to N eff . Sterile neutrinos in the O(10 − 100) MeV mass range can alter our cosmological history and are hence strongly constrained by observations related to BBN and the CMB. Here we distinguish three cases, depending on the lifetime τ i of the sterile neutrino N i [56], If N i decays significantly before BBN, its decay products are fully thermalised and merely lead to a shift in the overall temperature of the thermal bath, which 6 For a more complete listing see the pdgLive page on HNLs [40]. 7 In ref. [52] it was pointed out that imposing a generalised B − L symmetry can lead to a parametric suppression of all lepton number violating observables. This suppression indeed happens for the 0νββ-decay and for the (Majorana) masses mi of the light neutrinos (where it is necessarily needed to allow for mixings Such short lifetimes require a sizeable mixing with SM neutrinos (above the region labeled 'cosmological history between BBN and CMB' in figure 1), which leads to a non-trivial interplay with the laboratory constraints discussed above (gray region in figure 1). 8 (ii) Long-lived N i . Heavy neutrino decays during BBN would directly alter the formation of light elements. If N i decays after BBN (but before CMB decoupling) it can impact the post-BBN cosmological history. In extreme cases, the non-relativistic sterile neutrinos can even come to dominate the energy budget of the Universe. Moreover, their decay leads to an entropy injection into the SM thermal bath. This leads to an upper bound on the mixing α |U αi | 2 of N i with the SM neutrinos ν α [19] (region labeled 'cosmological history

JHEP01(2021)200
between BBN and CMB' in figure 1). 9 A further constraint arises from the effective number of relativistic degrees-of-freedom N eff at the time of BBN. As was demonstrated in ref. [8], at least two out of the three sterile neutrinos temporarily reach thermal equilibriumand consequently a sizeable abundance -in the early Universe. This leads to a significant contribution to N eff during BBN if the sterile neutrino is relativistic at decoupling. We find the resulting upper bound on the mixing between active and sterile neutrinos to be weaker than the constraint derived in [19] in the parameter space of interest.
(iii) Quasi-stable N i . N i lifetimes beyond the time of CMB decoupling (∼ 10 12 s) are highly constrained by CMB observations [21,62] (region labeled 'CMB constraints' in figure 1), the impact of their decays on the intergalactic medium (IGM) [20,63] (region labeled 'IGM constraints' in figure 1), and the produced X-rays [64]. The constraints can be avoided for sufficiently small θ ai , which suppresses both the thermally produced abundance and the N i decay rate. In the mass range considered here, such a long lifetime requires mixing angles that are so tiny that the amount of thermally produced N i is negligible for all practical purposes (white region at the bottom of figure 1). 10 This leads to the three distinct regions I)-III) of the parameter space which survive both the laboratory and cosmological constraints. In scenario I) all three sterile neutrinos belong to population (i). They have relatively large mixing angles, and thermalise and decay before BBN. This region is found by applying the bound on the lifetime of sterile neutrinos from refs. [59][60][61]67]. In scenario II), two of the heavy neutrinos N 2 and N 3 belong to population (i). The third heavy neutrino N 1 features significantly smaller mixings θ a1 with the SM states and belongs to population (ii). N 1 avoids thermalisation [8] and obeys the bounds derived in ref. [19]. More precisely, we use the bound depicted in figure 2 of ref. [19], which leaves open a window for M 1 50 MeV and |θ ai | 2 < 10 −14 . However, as we will see in the following section 2.3, this window is closed if the effect of N i decays on photodisintegration of nuclei is taken into account. Scenario III) is similar to scenario II), but N 1 has even smaller mixings and is part of population (iii).

Additional constraints from photodisintegration of nuclei
Further, the N i decay can also disintegrate nuclei in the primordial plasma after BBN. The resulting bound strongly depends on the hadronic branching ratio of the decay, which vanishes for the mass range considered here, and does not affect any of the points in our sample [68]. Hence, we only have to take into account electromagnetic decay channels, which we incorporate via the procedure described in [69,70] by running the public code ACROPOLIS [71]. 11 On that note, we first determine the non-thermal photon/electronspectra by solving the full cascade equation [72] with the appropriate source terms ∝ 9 Note that the setup of ref. [19] contains only one sterile neutrino which couples exclusively to νe. Taking into account the actual flavour structure in the couplings, we impose the bound derived in ref. [19] on the mixing summed over all SM flavours. 10 For masses in the keV range, eq. (2.5) permits mixing angles that are large enough that thermally produced Ni can make up a considerable fraction of the DM and the bounds summarised in [65,66] should be applied. 11 Decays into SM neutrinos do not lead to photodisintegration and therefore can be neglected.

JHEP01(2021)200
n i /τ i . The resulting spectra are then used to determine the late-time modifications of the nuclear abundances via photodisintegration by solving the appropriate non-thermal Boltzmann equation. Finally, we compare the resulting abundances with the most recent set of observations [40,73]. Specifically, we use (2.8) The resulting constraints are shown in the lower right panel of figure 1 (orange) and we find that these additional limits are particularly important for closing the region of parameter space between the solid blue and purple line, i.e. the region that is otherwise neither excluded by CMB observations nor by a modified cosmological history between BBN and CMB.

Supernovae bound
The detection of SN 1987A neutrinos arriving over an interval of about 10 s, in agreement with the predictions of a core-collapse supernova with the standard cooling scenario, imposes constraints on the existence of light BSM particles which would constitute an additional channel of energy-loss, shortening the duration of the neutrino burst [74]. This has in particular been used to constrain axions [74], dark photons [75], and sterile neutrinos of different mass ranges [55,57,58,[76][77][78][79]. The constraints are particularly relevant for the mixing with ν τ , since the laboratory constraints are weakest in this case. We indicate the constraints found in [55] by the yellow shaded area in the bottom left panel of figure 1. However, as it has been recently pointed out in ref. [80], these bounds rely on the standard core-collapse supernova model. If instead the supernova is modelled by a collapseinduced thermo-nuclear explosion [81], the observed neutrino signal could stem from the accretion disk and would be insensitive to the cooling rates. With this in mind, we do not apply the supernova bounds of ref. [55] in our main analysis, but emphasize that this region of parameter space can be fully probed in the near future -both by laboratory and astrophysical observations.

Viable parameter space
The photodisintegration bound introduced in section 2.3 excludes all points of type (ii) in the mass range considered here and therefore rules out scenario II). As already stated in the introduction, the phenomenology of scenario III) corresponds to that of the much-studied νMSM and shall not be further investigated here. This leaves us with scenario I). A priori it is not clear whether there are any viable parameter choices for which all bounds can be fulfilled simultaneously. This is non-trivial because neutrino oscillation data restricts the flavour mixing pattern, i.e. the range of allowed values for U 2 ai /U 2 i , meaning that it may not be possible to fit all three N i into the allowed (white) parameter regions in figure 1. It is well-known that this considerably constrains the range of allowed masses M i below the kaon mass in scenario III), as the constraints on N 2 and N 3 in this scenario are practically

JHEP01(2021)200
identical to those in the model with only two heavy neutrinos because N 1 cannot make a measurable contribution to the seesaw mechanism. 12 When combined with direct search data and bounds from BBN, this practically rules out most heavy neutrino masses below ∼ 350 MeV, with a few small windows between 100 MeV and 300 MeV left open if the light neutrino mass ordering is normal [24]. 13 For scenario I) this question has been studied in [7], where it was found that the combination of all experimental and cosmological bounds indeed leaves a sizeable region of viable parameter space with M i below the pion mass. The reason is that neutrino oscillation data in this case allows for larger values of U 2 τ i /U 2 i than in scenario III); this permits the heavy neutrinos to decay through their mixing with the third SM generation before BBN while respecting the direct search bounds on U 2 e and U 2 µ , which are much stronger than those on U 2 τ . Hence, HNLs with masses well below the kaon mass should be "tau-philic". However, the analysis in [7] did not include the supernova bound discussed in section 2.4. A complete scan of the allowed parameter region is numerically extremely expensive because of the high dimensionality of the parameter space (18 free parameters) and the complicated shape that the sub-manifolds defined by the various experimental constraints in the mass region considered here form in this space. Instead, we perform a limited scan with randomised parameter choices. We use the radiatively corrected [83] Casas-Ibarra parameterisation [84]. For the mass splittings and the complex angles in the Casas-Ibarra parameterisation, we alternate between drawing our parameters from a linear versus a logarithmic distribution, as in ref. [9]. We apply all experimental and cosmological constraints summarised above. For the experimental bounds, we use the simple strategy adapted in refs. [24,85] and interpret the exclusion regions published by the experimental collaborations as hard cuts (rather than using full likelihood functions as in ref. [7]), which is sufficient for the purpose of this work. For the lifetime constraints from BBN we use the results from [59]. 14 We show a representative set of viable parameter points (indicated by green dots) that are consistent with all experimental and cosmological constraints in figure 1. These all correspond to the normal ordering of the SM neutrinos and require U 2 τ i > U 2 ei , U 2 µi . Taking into account that each parameter point is represented by a triplet of points in figure 1 (one for each sterile flavour), applying the supernova bound from ref. [55] would exclude all points shown. However, as pointed out in section 2.4, this bound strongly relies on the underlying model for the supernova explosion. 12 This can be seen by inserting the largest mixing angels for population (iii) in figure 1 into the seesaw formula (2.3). An estimate m lightest ∼ U 2 1 M1 yields values m lightest < 10 −12 eV for Mi in the MeV range, in which case the constraints on N2 and N3 are practically identical to those in the minimal model with two heavy neutrinos, cf. figure 11 in [7]. Bearing in mind the caveat already pointed out in footnote 4 the results found in section 2 of [24] can therefore be applied to scenario III). In the νMSM larger values of U 2 1 are allowed because M1 is in the keV range, but a similar conclusion can be drawn [82]. 13 The authors of [24] used the BBN bounds from [59] in their global analysis. The updated bounds on τi from [61] are stronger for Mi between the pion and kaon masses, and including them in a global analysis is likely to close some of these windows. 14 The more recent bounds from [60] agree with those. In [61] it was pointed out that the bound on τi becomes considerably stronger when one includes the effect of mesons that are produced in heavy neutrino decays on the primordial light element abundances. However, this does not apply to the mass range considered here because the heavy neutrinos can only decay into purely leptonic final states.

JHEP01(2021)200 3 Baryogenesis
We now proceed to compute the baryon asymmetry for all viable parameter points found in our scan, using the set of quantum kinetic equations given in section 2 of ref. [9] to describe the evolution of the heavy neutrino abundances and lepton asymmetries in the early universe. 15 We assume that the radiation dominated epoch of the cosmic history started with a matter-antimatter symmetric primordial plasma in which all SM particles were in thermal equilibrium at a temperature T R that was much hotter than the temperature T sph 131 GeV [89], above which electroweak sphalerons efficiently convert L and B into each other [90]. In inflationary cosmology, this is expected because pre-inflationary asymmetries would be diluted very efficiently by the cosmic expansion. We moreover take the initial abundance of the heavy neutrinos to be negligible. The Lagrangian in eq. (1.1) then contains all the necessary ingredients to generate the baryon asymmetry of our Universe: the heavy neutrinos are generated from thermal interactions in the plasma through their Yukawa couplings. In this out-of-equilibrium situation, the interplay of coherent neutrino oscillations and decoherent scatterings mediated by the CP -violating Yukawa couplings can generate a lepton asymmetry that is partially converted into a baryon asymmetry by the sphalerons. For the M i under consideration here, and in view of the experimental constraints on the U i , this process happens very slowly. If at least one heavy neutrino has not reached thermal equilibrium at T ∼ T sph , then the baryon asymmetry is preserved ("frozen in") at lower temperatures. This freeze-in leptogenesis mechanism, also known as Akhmedov-Rubakov-Smirnov (ARS) leptogenesis [14], has been studied by many authors, a review is e.g. given in ref. [91]. Our goal is to study the question of whether the observed baryon asymmetry of the Universe can be explained in scenario I) while respecting the constraints discussed in section 2 if all M i are smaller than the pion mass.
Our results are shown in figure 2. Remarkably, if we consider the population where all three sterile neutrinos decay before BBN, the predicted baryon asymmetry is generically in the correct ball-park to explain the observed value or larger. This is far from trivial since it is well known that marginal changes in the model parameters can lead to drastic changes in the resulting baryon asymmetry, due to the fine balance between generation and wash-out of the asymmetries.

Discussion and conclusion
We study the type-I seesaw model with three heavy neutrinos N i with masses between a few MeV and the pion mass. This part of the parameter space is relatively little studied because it is ruled out by the combination of direct searches for heavy neutrinos, cosmological constraints and light neutrino oscillation data if one considers the minimal seesaw with only two heavy neutrinos or the νMSM, where neutrino oscillation data prohibits large hierarchies between the mixings U 2 ai with individual SM generations a = e, µ, τ . However, in a general framework with three right-handed neutrinos, neutrino oscillation data permits a comparably large mixings with the third SM generation that allows the N i to decay before BBN while respecting the tighter direct search constraints on the mixings with the first two generations. There are three scenarios with all heavy neutrino masses M i below the pion mass that are potentially allowed by previously published constraints: I) all three N i decay before BBN, II) two N i decay before BBN and the third one decays between BBN and the CMB decoupling, and III) two N i decay before BBN and the third one has a lifetime that greatly exceeds the age of the universe.
Scenario III) can be ruled out when all M i are in the MeV range by noticing that the allowed mixing angles for the quasi-stable N 1 are so tiny that this particle effectively decouples, and the constraints on N 2 and N 3 are practically identically to those in the minimal model with two heavy neutrinos, where this mass range is known to be ruled out. We further find that scenario II) is ruled out by the effect that photons produced in a cascade from the long-lived N i decay would have on the disintegration of light elements in the IGM. In scenario I) a representative randomised parameter scan shows that there are viable parameter values for which the N i can avoid all constraints from laboratory experiments and cosmology for M i > 50 MeV. All these points can potentially be ruled out by the observed neutrino flux from the supernova 1987a, but this conclusion depends on the modelling of the supernova explosion. Quite surprisingly most of the viable points give a final baryon asymmetry in the correct ball-park to explain the observed value. Given the well-known strong sensitivity of the relevant Boltzmann equations to small changes in the parameters, this is a highly non-trivial result.

JHEP01(2021)200
Our results show that heavy neutrinos masses down to a few tens of MeV are allowed by experimental and cosmological constraints if all of them decay well before the onset of BBN in the early universe, and if they primarily mix with the third SM generation. N i in this mass range can be found in the decay of light mesons, in particular pions, which are produced in large numbers in accelerator based searches. While the LHC main detectors have no sensitivity in this mass range due to both, the low transverse momentum of the events and the long lifetime of the N i , dedicated detectors in the forward directions, such as FASER [92], instrumentation in the beam pipe [93] or the Forward Physics Facility [94], could potentially search for such N i . Moreover, the mass range below the pion mass can be probed by fixed target experiments. For instance, the neutrino beams of DUNE [95] and T2K [96] contain a fraction of N i , 16 the decay of which (e.g. N i → e + e − ν a ) could be observed in their near detectors [97,98] or with additional detectors [99]. The channel with a = τ can be used to probe the coupling U 2 τ i to the third generation, which tends to be the largest to ensure decay BBN while respecting the stronger experimental bounds on U 2 ei and U 2 µi . It would also be worth to investigate searches at the existing SBL [100] and NuMi [101] beam lines at FNAL.
This work provides a proof-of-existence for a viable leptogenesis scenario in this mass range, but this is by no means an exhaustive study. Due to the high-dimensional parameter space, this requires more sophisticated numerical techniques, but we hope that the results presented here will trigger further work in this direction. This will be crucial in guiding experimental effort in fully testing freeze-in leptogenesis as the mechanism to generate the baryon asymmetry of our Universe.