Phenomenology of GeV-scale heavy neutral leptons

We review and revise phenomenology of the GeV-scale heavy neutral leptons (HNLs). We extend the previous analyses by including more channels of HNLs production and decay and provide with more refined treatment, including QCD corrections for the HNLs of masses O\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{O} $$\end{document}(1) GeV. We summarize the relevance of individual production and decay channels for different masses, resolving a few discrepancies in the literature. Our final results are directly suitable for sensitivity studies of particle physics experiments (ranging from proton beam-dump to the LHC) aiming at searches for heavy neutral leptons.


Contents
We review and revise phenomenology of the heavy neutral leptons (HNLs) with masses in the GeV range. The interest to these particles has recently increased, since it was recognized that they are capable of resolving 3 major observational BSM phenomena: neutrino oscillation, baryon asymmetry of the universe and dark matter [1,2] (for review see e.g. [3,4], [5,Chapter 4] and references therein). Several particle physics experiments, that put the searches for heavy neutral leptons among their scientific goals, have been proposed in the recent years: DUNE [6], NA62 [7][8][9] SHiP [5,10], CODEX-b [11], MATHUSLA [12][13][14][15], FASER [16][17][18]. The searches for HNLs (also often called "Majorana neutrinos" or "sterile neutrinos") have been performed and are ongoing at the experiments LHCb, CMS, ATLAS, T2K, Belle (see e.g. [19][20][21][22][23]25]) with many more proposals for novel ways to search for them [24,. This interest motivates the current revision. The information relevant for sensitivity studies of the GeV-scale HNLs is scattered around the research literature [37,39,[48][49][50][51][52][53][54][55][56][57] and is sometimes controversial. We collect all relevant phenomenological results and present them with the unified notation, discussion of the relevance of the individual channels and references to the latest values of phenomenological parameters (meson form factors) that should be used in practical application. The relevance of individual channels depending on the masses of HNLs is present in the resulting Table 5. We also discuss existing discrepancies in the literature, pointing out the way of obtaining the correct results and analyze new channels of production and new modes of decay neglected in the previous literature.

General introduction to heavy neutral leptons
Heavy neutral leptons or sterile neutrinos N are singlets with respect to the SM gauge group and couple to the gauge-invariant combination (L c α ·H) (where L α , α = 1, . . . , 3, are SM lepton doublet,H i = ε ij H * j is conjugated SM Higgs doublet) as follows L Neutrino portal = F α (L α ·H)N + h.c. , (1.1) with F α denoting dimensionless Yukawa couplings. The name "sterile neutrino" stems from the fact that the interaction (1.1) fixes SM gauge charges of N to be zero. After electroweak symmetry breaking the SM Higgs field gains nonzero vacuum expectation value v and interaction (1.1) provides heavy neutral leptons and SM (or active) neutrinos -with the mixing mass term (v = 246 GeV) The truly neutral nature of N allows one to introduce for it a Majorana mass term, consistent with the SM gauge invariance, resulting in the HNL Lagrangian at GeV scale The mass eigenstates of the active-plus-sterile sector are the mixtures of ν and N , but with small mixing angles and large splitting between mass scales of sterile and active neutrinos. The heavy mass eigenstates are "almost sterile neutrinos" while light mass eigenstates are "almost active neutrinos". In what follows we keep the same terminology for the mass states as for the gauge states. As a result of mixing, HNL couples to the SM fileds in the same way as active neutrinos, except the coupling is strongly suppressed by the small mixing angles In (1.3) α are charged leptons of the three SM generations. The number of model parameters increase with the number of NHL (see e.g. reviews [3,4]). In particular in the model with 2 sterile neutrinos there are 11 free parameters and in the case of 3 sterile neutrinos there are 18 parameters [3]. Not all of them play important role in phenomenology. The collider phenomenology is sensitive only to masses of the HNL(s) and absolute values of mixing angles, |U α |. When sterile neutrinos are not degenerate in mass, in all the processes they are produced and decay independently, without oscillations between themselves, in contrast to the behavior of active neutrinos. So, from the phenomenological point of view it is enough to describe only 1 sterile neutrino, which needs only 4 parameters: sterile neutrino mass M N and sterile neutrino mixings with all three active neutrinos U α , Eq. (1.2).
The papers is organized as follows: in Section 2 we review the different HNL production channels; in Section 3 we discuss the most relevant NHL decay channels. The summary and final discussion is given in the section 4. Appendices provide necessary technical clarifications.

HNL production in proton fixed target experiments
In fixed target experiments (such as NA62, SHiP or DUNE) the initial interaction is proton-nuclei collision. In such collisions HNLs can be produced in a number of ways: a) Production from hadron's decays; b) Production from Deep Inelastic Scattering (DIS) p-nucleus interaction; c) Production from the coherent proton-nuclear scattering. Below we provide overview of each of the channels summarizing previous results and emphasizing novel points.

Production from hadrons
The main channels of HNL production from hadrons are via decays of sufficiently long-lived hadrons, i.e. the lightest hadrons of each flavour 1 . In the framework of the Fermi theory, the decays are inferred by the weak charged currents. One can also investigate the hidden flavored mesons J/ψ(cc, 3097), Υ(bb, 9460) as sources of HNLs. These mesons are short-lived, but 1.5-2 times heavier than the corresponding open flavored mesons, giving a chance to produce heavier HNLs.
As the region of HNL masses below that of the kaon is strongly bounded by the previous experiments (see [5] for details, reproduced in Fig. 1), in what follows we concentrate on production channels for HNL masses M N > 0.5 GeV.
HNLs are produced in meson decays via either 2-body purely leptonic decays (left panel of Fig. 2) or semileptonic decays (right panel of Fig. 2) [61,62]. The branching fractions of leptonic decays have been found e.g. in [49,51]. For the semileptonic decays only the processes with a single pseudo-scalar or vector meson in the final Figure 1. Existing limits and future prospects for searches for HNLs. Only mixing with muon flavour is shown. For the list of previous experiments (gray area) see [5]. Black solid line is the recent bounds from the CMS 13 TeV run [22]. The sensitivity estimates from prospective experiments are based on [27] (FCC-ee), [9] (NA62), [58] (SHiP) and [59] (MATHUSLA). The sensitivity of SHiP below kaon mass (dashed line) is based on the number of HNLs produced in the decay of D-mesons only and does not take into account contributions from kaon decays, see [58] for details. The primordial nucleosynthesis bounds on HNL lifetime are from [60]. The Seesaw line indicates the parameters obeying the seesaw relation |U µ | 2 ∼ m ν /M N , where for active neutrino mass we substitute m ν = ∆m 2 atm ≈ 0.05 eV [5].  state have been considered so far [51] (see also [55] and [37]) (where h P is a pseudo-scalar and h V is a vector meson) and their branching ratio has been computed. We reproduce these computations in the Appendix A paying special attention to the treatment of form factors. Finally, to calculate the number of produced HNLs one should ultimately know the production fraction, f (qq → h) -the probability to get a given hadron from the corresponding heavy quark. The latter can either be determined experimentally or computed from Pythia simulations (as e.g. in [57]).

Production from light unflavored and strange mesons
Among the light unflavored and strange mesons the relevant mesons for the HNL production are: 2 π + (ud, 139.6), K + (us, 494), K 0 S (ds, 498) and K 0 L (ds, 498). The only possible production channel from the π + is the two body decay π + → + α N with = e, µ. The production from K + is possible through the two-body decay of the same type. There are also 3-body decays K + → π 0 + α N and K 0 L/S → π − + α N . The resulting branching ratios for corresponding mesons are shown in Fig. 3. For small HNL masses the largest branching ratio is that of K 0 L → π − + α N due to the helicity suppression in the two-body decays and small K 0 L decay width.

Production from charmed mesons
The following charmed mesons are most relevant for the HNL production: D 0 (cū, 1865), D + (cd, 1870), D s (cs, 1968). BR(B→X+N) Figure 4. Dominant branching ratios of HNL production from different charmed and beauty mesons. For charged mesons two-body leptonic decays are shown, while for the neutral mesons decays are necessarily semi-leptonic. For these plots we take U e = 1, D 0 is a neutral meson and therefore its decay through the charged current interaction necessarily involves a meson in a final state. The largest branching is to K meson, owing to the CKM suppression |V cd |/|V cs | ≈ 0.22. Then the mass of the resulting HNL is limited as M N < M D − M K ≈ 1.4 GeV. For the charmed baryons the same argument is applicable: they should decay into baryons and the most probable is strange baryon, hence M N < M Λc − M Λ ≈ 1.2 GeV. Therefore these channels are open only for HNL mass below ∼ 1.4 GeV.
Charged charmed mesons D ± and D s would exhibit two body decays into an HNL and a charged lepton, so they can produce HNLs almost as heavy at themselves. The branching of D s → N + X is more than a factor 10 larger than any similar of other D-mesons. The number of D s mesons is of course suppressed as compared to D ± and D 0 mesons, however only by a factor of few 3 . Indeed, at energies relevant for cc production, the fraction of strange quarks is already sizeable, χs s ∼ 1/7 [63]. As a result, the two-body decays of D s mesons dominate in the HNL production from charmed mesons.

Production from beauty mesons
The lightest beauty mesons are B − (bū, 5279), B 0 (bd, 5280), B s (bs, 5367), B c (bc, 6276). Similarly to the D 0 case, neutral B-mesons (B 0 and B s ) decay through charged current with a meson in a final state. The largest branching is to D meson because of the values of CKM matrix elements (|V cb |/|V ub | ≈ 0.1). Thus the mass of the resulting HNL is limited: Charged beauty mesons B ± and B ± c have two body decays into HNL and charged lepton, so they can produce HNLs almost as heavy at themselves. Due to the CKM suppression the branching ratio of B + → N + + is significantly smaller than that of B c → N + . However, unlike the case of D s mesons, the production fraction of f (b → Inclusive branching: l = e, µ 11.0 ± 0.
Hence 1-meson modes contribute additionally 1.09 ± 0.12 Sum of other multimeson channels, n > 1: Inclusive branching: l = τ not known Dominant one-meson channels: pseudo-scalar meson Table 1. Experimentally measured branching widths for the main semileptonic decay modes of the B + and B 0 meson [63]. Decays to pseudoscalar (D) and vector (D * ) mesons together constitute 73% (for B + ) and 69% (for B 0 ). Charmless channels are not shown because of their low contribution B c ) has only been measured at LHC energies, where it is reaching few × 10 −3 [64]. At lower energies it is not known.

Multi-hadron final states
D and especially B mesons are heavy enough to decay into HNL and multimeson final states. While any single multi-meson channel would be clearly phase-space suppressed as compared to 2-body or 3-body decays, considered above, one should check that the "inclusive" multi-hadron decay width does not give a sizeable contribution.
To estimate relative relevance of single and multi-meson decay channels, we first consider the branching ratios of the semileptonic decays of B + and B 0 (with ordinary (massless) neutrino ν in the final state) where X are one or many hadrons. The results are summarized in Table 1. Clearly, by taking into account only single meson states we would underestimate the total inclusive width of the process (2.3) by about 20%.
In case of semileptonic decays in the HNL in the final state, the available phase space shrinks considereably, see Fig. 5. The effect of the mass can also be estimated by comparing the decays involving light leptons (e/µ) and τ -lepton in the final state. A comparison with SM decay rates into τ -lepton shows that 3-body decays into heavy sterile neutrinos are suppressed with respect to decays to light neutrinos. Thus inclusive semi-leptonic decay of flavoured mesons to HNLs are dominated by single-meson final states with the contributions from other state introducing small correction.

Quarkonia decays
Next we investigate the hidden flavored mesons J/ψ(cc, 3097) and Υ(bb, 9460) as sources of HNLs. These mesons are short-lived, but 1.5-2 times heavier than the corresponding open flavored mesons, giving a chance to produce heavier HNLs. We have studied these mesons in Appendix D, here we provide the summary of the results.
The number of HNLs produced from J/ψ decays is always subdominant to the number of HNLs produced in D-meson decays (for M N < m D ). Therefore, the range of interest is 2 GeV ≤ M N ≤ m J/ψ where this number should be compared with the number of HNLs produced via B-meson decays. The resulting ratio is given by  where we have adopted f (B) × BR(B → N + X) ∼ 10 −2 (c.f. Fig. 4, right panel) and used f (J/ψ) ∼ 10 −2 . The numbers in (2.4) are normalized to the 400 GeV SPS proton beam. One sees that J/ψ can play a role only below bb production threshold (as X bb tends to zero). For experiments where sizeable number of bb pairs is produced one can use the Υ decays to produce HNLs with M N 5 GeV. The number of thus produced HNLs is given by where N Υ is the total number of Υ mesons produced and we have normalized U 2 to the current experimental limit for M N > 5 GeV (c.f. Fig. 1). It should be noted that HNLs with the mass of 5 GeV and U 2 ∼ 10 −5 have the decay length cτ ∼ cm.

Production from baryons
Semileptonic decays of heavy flavoured baryons ( Table 2) produce HNLs. Baryon number conservation implies that either proton or neutron (or other heavier baryons) must be produced in the heavy baryon decay, which shrink by about 1 GeV the kinematical window for sterile neutrino. The corresponding heavy meson decays have an obvious advantage in this respect. Moreover, since both baryons and sterile neutrinos are fermions, only the baryon decays into three and more particles in the final state can yield sterile neutrinos, which further shrinks the sterile neutrino kinematical window with respect to the meson case, where two-body, pure leptonic decays can produce sterile neutrinos. Furthermore, light flavored baryons, strange baryons (see Table 2) can only produce HNLs in the mass range where the bounds are very strong already (roughly below kaon mass, see FIG. 1). Indeed, as weak decays change the strangeness by 1 unit, there the double-strange Ξ-baryons can only decay to Λ or Σ baryons (plus electron or muon and HNL). The maximal mass of the HNL that can be produced in this process is smaller than (M Ξ − − M Λ 0 ) 200 MeV. Then, Ω − baryon decays to Ξ 0 − N with the maximal HNL mass not exceeding M Ω − − M Ξ 0 350 MeV. Finally, weak decays of Λ or Σ baryons to (p, n) can produce only HNLs lighter than ∼ 250 MeV.
The production of HNL in the decays of charmed and beauty hyperons has been investigated in Ref. [52], which results have been recently checked in [65]. The number of such baryons is of course strongly suppressed as compared to the number of mesons with the same flavour. At the same time the masses of HNLs produced in the decay of charmed (beauty) baryons are below the threshold of HNL production of the corresponding charm (beauty) mesons due to the presence of a baryon in the final state. This makes such a production channel strongly subdominant. A dedicated studies for SHiP [57] and at the LHC [65] confirm this conclusion. It should be noted that Refs. [52,57] use form factors from Ref. [66] which are about 20 years old. A lot of progress has been made since then (see e.g. [67,68], where some of these form factors were re-estimated and a factor ∼ 2 difference with the older estimates were established).

HNL production from tau lepton
At centre of mass energies well above thecc threshold τ -leptons are copiously produced mostly via D s → τ + X decays. Then HNLs can be produced in τ decay and these decays are important in the case of dominant mixing with τ flavour (which is the least constrained, see [5,Chapter 4]). The main decay channels of τ are τ → N + h P/V , τ → N ανα and τ → ν τ α N , where α = e, µ. The computations of the corresponding decays widths are similar to the processes N → α h P/V (c.f. Appendix B.2) and purely leptonic decays of HNL (see Section 3.1.1). The results are Figure 6. HNL production channels: a) Drell-Yan-type process; b) gluon fusion; c) quarkgluon fusion.
where y i = m i /m τ , V U D is an element of CKM matrix which corresponds to quark content of the meson h P , f h and g h are pseudoscalar and vector meson decay constants (see Tables 8 and 9) and λ is the Källén function [69]: The results of this section fully agree with literature [51].

HNL production via Drell-Yan and other parton-parton scatterings
The different matrix elements for HNL production in proton-proton collision are shown in Fig. 6. Here we are limited by the beam energy not high enough to produce real weak bosons on the target protons. There are three type of processes: Drell-Yan-type process a), gluon fusion b) and W γ/g fusion c). Process b) starts to play an important role for much higher centre-of-mass energies [70,71], process a) and c) should be studied more accurately. Let us start with the process a) in Fig. 6. The cross section at the parton level is [72,73] where V qq is an element of the CKM matrix, N c = 3 is a number of colors and the centre-of-mass energy of the systemqq is given by where x 1 and x 2 are fractions of the total proton's momentum carried by the quark q and anti-quarkq respectively. The total cross section therefore is written as where f q (x, Q 2 ) is parton distribution function (PDF). The corresponding integral S( √ s, M N ) as a function of M N and the production probability for this channel are shown in Fig. 7. For numerical estimates we have used LHAPDF package [74] with CT10NLO pdf set [75]. This can be roughly understood as follows: PDFs peak at x 1 (see Fig. 8) and therefore the probability that the center-of-mass energy of a parton pair exceeds the HNL mass, √ s parton M N , is small. On the other hand, the probability of a flavour meson to decay to HNL (for |U | 2 ∼ 1) is of the order of few % and therefore "wins" over the direct production, especially at the fixed-target experiments where beam energies do not exceed hundreds of GeV. In case of the quark-gluon initial state (process c) in Fig. 6) the similar considerations also work and the resulting cross section is also small, with an additional suppression due to the 3-body final state. We see that the direct production channel is strongly suppressed in comparison with the production from mesons for HNLs with masses M N 5 GeV.

Coherent proton-nuclei scattering
The coherent scattering of a proton off the nuclei as a whole could be an effective way of producing new particles in fixed target experiments. There are two reasons for The suppression of the integral as compared to M N = 0 case is due to PDFs being small at x ∼ 1 and condition x 1 x 2 s > M 2 N . Total p-p cross section is taken from [63]. this. First, parton scattering in the electromagnetic field of the nuclei is proportional to Z 2 (where Z is the nuclei charge) which can reach a factor 10 3 enhancement for heavy nuclei. Second, the centre of mass energy of proton-nucleus system is higher than for the proton-proton scattering. The coherent production of the HNLs will be discussed in the forthcoming paper [76]. Here we announce the main result: the HNL coherent production channel is subdominant to the meson decay for all HNL masses and mixing angles (for HNL masses below 5 GeV). In case of SHiP on expects less than 1 HNL produced via coherent scattering for 10 20 PoT. Figure 9. Possible Feynman diagrams for the HNL production in the proton coherent scattering off the nuclei.

Summary
In summary, production of HNL in proton fixed target experiments occurs predominantly via (semi)leptonic decays of the lightest c-and b-mesons (Fig. 4). The production from heavier mesons is suppressed by the strong force mediated SM decays, while production from baryons is kinematically suppressed. Other production channels are subdominant for all masses 0.5 GeV ≤ M N ≤ 5 GeV as discussed in Sections 2.3-2.4.

HNL decay modes
All HNL decays are mediated by charged current or neutral current interactions (1.3).
In this Section we systematically revisit the most relevant decay channels. Most of the results for sufficiently light HNLs exist in the literature [37, 49-51, 53, 54]. For a few modes there are discrepancies by factors of few between different works, we comment on these discrepancies in due course. All the results presented below do not take into account charge conjugated channels which are possible for the Majorana HNL; to account for the Majorana nature one should multiply by 2 all the decay widths. The branching ratios are the same for Majorana case and for the case considered here.

3-body basic channels
Two basic diagrams, presented in the Fig. 10, contribute to all decays. For the charged current-mediated decay ( Fig. 10(a)) the final particles (U, D) could be either a lepton pair (ν α , α ) or a pair of up and down quarks (u i , d j ). For the neutral currentmediated decay f is any fermion. The tree-level decay width into free quarks, while unphysical by itself for the interesting mass range, is important in estimates of the full hadronic width at M N Λ QCD , see Section 3.2.2 below. Figure 10. Diagram for the HNL decays mediated by charged a) and neutral b) currents. For the decays N → ν α − α + α and N → ν α ν ανα both diagrams contribute, which leads to the interference (see Section 3.1.2).

Charged current-mediated decays
The general formula for the charged current-mediated processes [50,53,54,77] The factor N W = 1 for the case of the final leptons and N W = N c |V ij | 2 in the case of the final quarks, where N c = 3 is the number of colors, and V ij is the corresponding matrix element of the CKM matrix. The function I(x u , x d , x l ) that describes corrections due to finite masses of final-state fermions is given by  Fig. 11. It decreases with each argument.

Decays mediated by neutral current interaction and the interference case
Decay width for neutral current-mediated decay N → ν α ff depends on the type of the final fermion. For charged lepton pair l βlβ the results are different for the case α = β and α = β, because of the existence of the charge current mediated diagrams in the latter case. Nevertheless, the decay width can be written in the unified way, and N Z = 1 for the case of leptons in the final state or N Z = N c for the case of quarks. The values of C f 1 and C f 2 are given in the Table 3. This result agrees with [51,53,54].
In the case of pure neutrino final state only neutral currents contribute and the decays width reads Table 3. Coefficients C 1 and C 2 for the neutral current-mediated decay width.

Decay into hadrons
In this Section we consider hadronic final states for M N both below and above Λ QCD scale and discuss the range of validity of our results.

Single meson in the final state
At M N Λ QCD the quark pair predominantly binds into a single meson. There are charged current-and neutral current-mediated processes with a meson in the final state: and g h are the corresponding meson decay constants (see Appendix C.1), θ W is a Weinberg angle and the function λ is given by eq. (2.10). The details of the calculations are given in the Appendix B.2.
The decay width to the charged pseudo-scalar mesons (π ± , K ± , D ± , D s , B ± , B c ) is given by (3.6) in full agreement with the literature [51,53,54]. The decay width to the pseudo-scalar neutral meson (π 0 , η, η , η c ) is given by Our answer agrees with [51], but is twice larger than [53,54]. The source of the difference is unknown. 5 The HNL decay width into charged vector mesons (ρ ± , a ± 1 , D ± * , D ± * s ) is given by that agrees with the literature [51,53,54]. However, there is a disagreement regarding the numerical value of the meson constant g ρ between [51] and [53,54]. We extract the value of this constant from the decay τ → ν τ ρ and obtain the result that numerically agrees with the latter works, see discussion in Appendix C.1.3.
For the decay into neutral vector meson (ρ 0 , a 0 1 , ω, φ, J/ψ) we found that the result depends on the quark content of meson. To take it into account we introduce dimensionless κ h factor to the meson decay constant (B.36). The decay width is given by Our result for ρ 0 and results in [51] and [53] are all different. The source of the difference is unknown. For decays into ω, φ and J/ψ mesons we agree with [53]. The result for the a 0 1 meson appears for the first time. 6 The branching ratios for the one-meson and lepton channels below 1 GeV are given on the left panel of Fig. 12.

Full hadronic width vs. decay into single meson final state
Decays into multi-hadron final states become kinematically accessible as soon as M N > 2m π . To estimate their branching fractions and their contribution to the total decay width, we can compute the total hadronic decay width of HNLs, Γ had and compare it with the combined with of all single-meson states, Γ 1 meson . The total hadronic decay width can be estimated via decay width into quarks (Sections 3.1.1-3.1.2) times the additional loop corrections. 6 Refs. [53,54] quote also two-body decays N → ν α h 0 V , h 0 V = K * 0 ,K * 0 , D * 0 ,D * 0 , with the rate given by (3.9) (with a different κ). This is not justified, since the weak neutral current does not couple to the corresponding vector meson h 0 V at tree level. The QCD loop corrections to the tree-level decay into quarks have been estimated in case of τ lepton hadronic decays. In this case the tree level computation of the τ decay to two quarks plus neutrino underestimates the full hadronic decay width by 20% [78][79][80]. The loop corrections, ∆ QCD , defined via have been computed up to three loops [80] and is given by: where α s = α s (m τ ). 7 We use (3.11) with α s = α s (M N ) as an estimation for the QCD correction for the HNL decay, for both charged and neutral current processes. We expect therefore that QCD correction to the NHL decay width into quarks is smaller than 30% for M N 1 GeV (Fig. 13). Full hadronic decay width dominates the HNL lifetime for masses M N 1 GeV (see Fig. 12). The latter is important to define the upper bound of sensitivity for the experiments like SHiP or MATHUSLA (see Fig. 1). This upper bound is defined by the requirements that HNLs can reach the detector.

Multi-meson final states
When discussing "single-meson channels" above, we have also included there decays with the ρ-meson. By doing so, we have essentially incorporated all the two-pion decays N → π + π 0 − for M N > m ρ . Indeed, we have verified by direct computation The same for neutral current channels. of N → π + π 0 − that they coincide with N → ρ + − for all relevant masses (Fig. 14).
Of course the decay channel to two pions is also open for 2m π < M N < m ρ , but its contribution there is completely negligible and we ignore this in what follows. Figs. 12 and 15 demonstrate that one-meson channels are definitely enough for all the hadronic modes if sterile neutrino mass does not exceed 1 GeV. The ratio between the combined decay width into single-meson final states (π ± , π 0 , η, η , ρ ± , ρ 0 , ω+φ, D s ) and into quarks is shown in Fig. 15. 8 One sees that the decay width into quarks is larger for M N 2 GeV, which means that multi-meson final states are important in this region.  The main expected 3-and 4-body decays channels of HNL and decay are presented in Table 4. In this table we also add information about multimeson decays of τ because they give us information about decay through charged current of the HNL of the same mass as τ -lepton. The main difference between HNL and τ -lepton comes from the possibility of the HNL decay through the neutral current, which we discuss below.

3-body decays
As one observes, the main hadronic channels of the τ are n-pions channels. Decay channel into 2 pions is the most probable, but there is a large contribution from the 3 pions channels and still appreciable contribution from the 4 pions ones. For bigger masses the contribution from the channels with higher multiplicity become more important as Fig. 15 demonstrates.
The decay into kaons is suppressed for the τ -lepton. For some channels like τ → ν τ K or τ → ν τ Kπ this suppression comes from the Cabibbo angle between s and u quarks. The same argument holds for HNL decays into lepton and D meson, but not in D s . The decays like τ → ν τ K − K 0 are not suppressed by CKM matrix and still are small. We think that this is because for such decays the probability of the QCD string fragmentation into strange quarks is much smaller than into u and d quarks for the given τ -lepton mass (see diagram a) in Fig. 16). At higher masses the probability of such fragmentation should be higher, but still too small to take it into account. On the other hand, the HNL decay into two kaons can give a noticeable contribution, because of existence of the neutral current decay (see diagram b) in Fig. 16), for which the previous arguments do not apply.

Summary
In this paper we revise the phenomenology of the heavy neutral leptons, including both their production and decays. We concentrated on the HNL masses up to O(10) GeV.
The mechanisms of the HNL production are secondary decays of the hadrons produced in the initial collision (Section 2.1), production in proton-nucleon collision (Section 2.3) and coherent scattering of the proton off nuclei (Section 2.4). Of these mechanisms the production from the lightest flavored mesons dominate at all masses of interest. Production from baryons is not efficient at any HNL mass (see discussion in Section 2.1.6). The main production channels above the kaon threshold are production from D mesons for M N 2 GeV and production from the beauty mesons for M N m Υ . For leptonic decays and two body semileptonic decays, the calculations are performed in Appendix A. Our results agree with [51], for the case of the pseudoscalar and vector mesons we present the simplified version of the final formulas. We additionally analyzed the HNL production in B meson decays including multimeson final states, that were not previously discussed. We estimate that contribution of the multimeson final state give not more than 20% of production from B mesons (Section 2.1.4).
The HNL are unstable and decay into light SM particles which can be detected. The HNL decay channels with branching ratio above 1% in the region M N < 5 GeV are summarized in the Table 5. For each channel we present the mass at which it opens, mass range where it is relevant and maximal branching ratio. The total decay width and the lifetime are summarized in Fig. 17.
All HNL decay channels can be divided into purely leptonic and semileptonic (hardonic) ones. The decay widths into leptons are given by (3.1), (3.4), (3.5) and are in full agreement with the literature [51,53].
For HNL masses above m π semileptonic decay channels quickly start to dominate, the hardonic branching ratio reaches ∼ 70% at M N 1 GeV. Single-meson final states (including decay into on-shell ρ mesons) saturate hadronic decay width till about 1.5 GeV (Fig. 15). In the HNL mass region 2 − 5 GeV from 50% to 80% of the semileptonic decay width is saturated by multimeson states. For completeness we summarize all relevant hadronic form factors in Appendices.
Our final results are directly suitable for sensitivity studies of particle physics experiments (ranging from proton beam-dump to the LHC) aiming at searches for heavy neutral leptons. Table 5: Relevant HNL decay channels. Only channels with the branching ratio above 1% covering the HNL mass range up to 5 GeV are shown. The numbers are provided for |U e | = |U µ | = |U τ |. For neutral current channels (with ν a in the final state) the sum over neutrino flavors is taken, otherwise the lepton flavor is shown explicitly. Columns: (1) the HNL decay channel. Notation (nπ) a means a system of n pions with the total charge a. (2) The HNL mass at which the channel opens; (3) The HNL mass starting from which the channel becomes relevant. For multimeson final states we provide our "best-guess estimates"; (4) HNL mass above which the channel contributes less than 1%; "-" means that the channel is still relevant at M N = 5 GeV, "?" means that we could not estimate the relevance of the channel; (5) The maximal branching ratio of the channel for M N < 5 GeV; (6)

A HNL production from hadrons
The calculation of weak decays, involving hadrons is summarized in [81]. In the absence of QED and QCD corrections the effective weak interaction Lagrangian at low energies can be written as where the charged current terms have the form and V U D is the CKM element which corresponds to quark flavour transition in hadronic current. For neutral current the interaction has the same form where summation goes over all fermions, and I 3f is the fermion isospin projection and Q f is its electric charge (Q e = −1). In the following Sections we describe different processes with HNL and hadrons.

A.1 Leptonic decay of a pseudoscalar meson
Consider a decay of pseudoscalar meson h into charged lepton and HNL: Fig. 2. The corresponding matrix element is given by where the corresponding quark contents of meson h is |h = Ū D . In order to fix the notations we remind that the charged meson coupling constant, f h , for a pseudoscalar meson constructed from up (U ) and down (D) type quarks is defined as where p µ is 4-momentum of the pseudo-scalar meson h. The numerical values of the decay constants for different mesons are summarized in Table 8. After standard calculation one finds the decay width of this reaction where y = m /m h , y N = M N /m h and λ is given by (2.10).

A.2 Semileptonic decay of a pseudoscalar meson
The process with pseudoscalar or vector meson h P/V in the final state is mediated by the current that has V − A form (see right diagram in Fig. 2). Properties of the hadronic matrix element h P/V J µ hadron h depend on the type of final meson h [82]. In the case of pseudoscalar meson only vector part of the current plays role: where q µ = (p − p ) µ is a transferred momentum and For the case of a vector meson h V in the final state both vector and axial part of the current contribute. The standard parametrization with form factors is 16) where µ is a polarization vector of the vector meson h V .
For the decay into vector meson the expression is more bulky, where I V,F G are parts of the decay width that depend on the F G form factors combination 9 and C K is a Clebsh-Gordan coefficient, C K = 1/ √ 2 for decays into ρ 0 and C K = 1 for all other cases in this paper. It turns out that I V,gf = I V,ga + = I V,ga − = 0, the other terms are given by (A.26) 9 In this computation we take all form factors as real-valued functions.
I V,f a + = 1 12 (A.28) where the notation is the same as in Eqs. (A. 18-A.20) and An important symmetry of the low-energy theory of strong interactions is the socalled G-symmetry which is a combination of the charge conjugationĈ and rotation of 180 • around the y axis in the isotopic spaceR y . 10 The operation of charge conjugation acts on bilinear combinations of fermions f 1 , f 2 as follows: R y acts on the isospin doublet asR Acting on pion states, that are pseudoscalar isovectors, one getŝ so any pion is an odd state under G-symmetry. As a consequence, for the system of n pionsĜ For ρ mesons, which are vector isovectors, G-parity is positive, while for a 1 mesons, which are pseudovector isovectors, G-parity is negative,

B.1.2 Classification of currents
Unflavoured quarks system interacts with electromagnetic field, W -and Z-bosons through currents Table 6. Properties of axial and vector currents under G-symmetry.
To divide currents (B.13)-(B.15) into G-odd and G-even parts let us introduce isoscalar and isovector vector currents 20) and isoscalar and isovector axial currents Currents (B.18)-(B.23) have certain G-parity presented in Table 6. Using these currents one can rewrite physical currents as

B.1.3 Connection between matrix elements
G-even part of currents (B.24)-(B.26) belongs to one isovector family, therefore there is an approximate connection between matrix elements for the system of even number of pions or ρ-meson, The special case to mention here is |2π 0 state. In V π 0 π 0 vertex, where V = γ/Z, system of 2 pions should have total angular momentum J = 1. Pions are spinless particles, so their coordinate wavefuction has negative parity, which is forbidden by the Bose-Einstein statistics. Therefore This result is equivalent to the prohibition of the ρ 0 → 2π 0 decay. G-odd parts of the currents (B.24)-(B.26), see Table 6, belong to one isoscalar and one isovector families, so there is only one relation between matrix elements for the system of odd number of pions or for a 1 -mesons, (B.29) The last formula can be simplified in the case of the one-pion or a 1 state. The direct interaction between photon and π 0 is forbidden because of the C symmetry, while photon-to-a 1 interaction violates both P and C symmetry. Therefore, the matrix element 0| J EM µ |π/a 1 = 0 and All the approximate relations discussed above hold up to isospin violating terms of order (m π + − m π 0 )/m π ∼ 3.4%.

B.2 HNL decays to a meson and a lepton
There are 4 types of this decay: N → α + h P/V and N → ν α + h P/V , where h P and h V are pseudoscalar and vector mesons respectively. Reaction N → α + h P is closely related to the process calculated in Section A.1. It utilizes the same matrix element and differs only by kinematics. Using the same notation, the decay width is and function λ is given by eq. (2.10).
In case of the neutral current-mediated decay N → ν α + h P the hadronic matrix element reads (see Section C.1.1 for details) where p µ is the 4-momentum of the pseudo-scalar meson h, J Z µ current is given by Eq. (B.15). The decay width is where x h = m h /M N and f h are neutral meson decay constants presented in the right part of Table 8. Consider the process N → α + h V . For the vector meson the hadronic matrix element of the charged current is defined as  Figure 19. Diagram for the HNL decay into 2 pions.
where ε µ (p) is polarization vector of the meson and g h is the vector meson decay constant. The values of the g h are given in Table 9. Using previous notations, the decay width of this process is Finally, to calculate HNL decay into neutral vector meson N → ν α + h V we define the hadronic matrix element as where g h is the vector meson decay constant and κ h is dimensionless correction factor, their values are given in Table 9. For the decay width one obtains

B.3 HNL decays to a lepton and two pions
For the case of 2 pions the matrix element of the axial current is equal to zero, so the general expression for matrix element is (c.f. (A.13)) where J µ is one of the currents (B.13)-(B.15) and q µ = (p − p ) µ . Because of isospin symmetry (B.27) the form factors are related as Electromagnetic current conservation q µ J µ = 0 implies f EM − (q 2 ) = 0. Therefore all the matrix elements could be expressed via only one form factor, called pion electromagnetic form factor, Pion electromagnetic form factor is related to the cross section of reaction e + e − → 2π as where β π (s) = 1 − 4m 2 π /s, so it is well-measured experimentally. There are a lot of data on electromagnetic form factor [85][86][87][88][89][90], which agree with each other. Good description of the data is given by the vector-dominance model (VDM), see Fig. 18 [85] and Appendix F for model description.
The decay into 2 pions is significantly enhanced by the ρ-resonance. It turns out, that this is the dominant channel, see Fig. 14 with comparison of the decay width of HNL into 2 pions and into ρ-meson. Therefore, one can replace the decay into 2 pions with 2-body decay into ρ-meson.

C Phenomenological parameters
In this Section we summarize parameters used in this work. Values of the CKM matrix elements are given in Table 7.  Table 7. CKM matrix elements [63] adopted in this work.

C.1 Meson decay constants
The decay constants for charged pseudoscalar mesons are defined by Eq. (A.10), Values of f h (Table 8) are measured experimentally and/or obtained by lattice calculations [91].
Meson decay constants for the mesons with the same-flavour quarks are defined by Eq. (B.32). There is a discrepancy regarding their values in the literature, therefore we have computed them directly (see Appendix C.1.1). The results of these computations are given in the right column of Table 8. The meson decay constants for neutral mesons consisted of quarks of different flavors (such as K 0 , D 0 , B 0 , B s ) are not needed in computing HNL production or decay, we do not provide them here.
For vector charged mesons the decay constants g h are defined by Eq. (B.34). In the literature they often appear as f h , connected to our prescription by mass of the meson g h = f h m h . Their values are presented in Table 9. For vector neutral mesons the decay constants g h and dimensionless factors κ h are defined by Eq. (B.36). Their values are presented in Table 9 as well.    ). Decay constants for D * (s) mesons in [94] show large theoretical uncertainty, we quote only the average value here.

C.1.1 Decay constants of η and η mesons
To describe HNL decays into η and η mesons we need to know the corresponding neutral current decay constants, that we define as (B.32) where where p µ is the 4-momentum of the pseudo-scalar meson h, J Z µ current is given by Eq. (B.15). The choice of the additional factor (−1/ √ 2) is introduced in order to obtain f π 0 = f π ± and f π 0 > 0, see discussion below. Taking into account that for pseudoscalar mesons only axial part of the current contributes to this matrix element we can write the matrix element as The relevant decay constants are f 0 and f 8 , they come from the set of extracted from experiments decay constants defined as [96] 0| J a µ |h = if a h p µ , where J a µ =qγ µ γ 5 λ a √ 2 q with λ a being the Gell-Mann matrices for a = 1 . . . 8 and The overall factor in λ 0 is chosen to obey normalization condition Tr(λ a λ b ) = 2δ ab . Within the chiral perturbation theory (χPT) (see [97] and references therein), the lightest mesons corresponds to pseudogoldstone bosons φ a , that appear after the spontaneous breaking of U L (3) × U R (3) symmetry to group U V (3). States φ a are orthogonal in the sense where and f a b are corresponding decay constants. Using we can rewrite the axial part of the weak neutral current (B.32) as a linear combination of the J 0 µ , J 3 µ and J 8 and f h is given by For example, π 0 meson corresponds to φ 3 state in χPT, so f 0 π 0 = f 8 π 0 = 0 and Eq. (C.8) gives f π 0 = f 3 π 0 = f π + because of isospin symmetry, in full agreement with Eq. (B.30).
For η and η application of Eq. (C.8) is not so straightforward. These mesons are neutral unflavoured mesons with zero isospin and they can oscillate between each other. So η and η do not coincide with any single φ a state. Rather they are mixtures of φ 0 and φ 8 states. In real world isospin is not a conserved quantum number, so φ 3 state also should be taken into account, but its contribution is negligible [98], so we use f 3 η = f 3 η = 0. Another complication is U (1) QCD anomaly for J 0 µ current that not only shifts masses of corresponding mesons but also contributes to f 0 h meson constant. To phenomenologically take into account the effect of anomaly it was proposed to use two mixing angles scheme [99], Taking parameter values from the recent phenomenological analysis [96], (2)f π , f 0 = 1.14(5)f π , θ 8 = −21.2(1.9) • , θ 0 = −6.9(2.4) • , (C. 10) we find These numbers should be confronted with the values quoted in [51] and [53].

C.1.3 Decay constant of ρ meson
There are 2 parametrizations of the ρ charged current matrix element using g ρ , defined by (B.34), or f ρ , which is related to g ρ are f ρ = g ρ /m ρ . The value of the decay constant can be obtained by 2 methods: from ρ → e + e − using approximate symmetry (B.27) or from the τ -lepton decay. Results, obtained in Ref. [93] by these two method differ by about 5%, f ρ,ee = 220(2) MeV and f ρ,τ = 209(4) MeV. We calculate and get g ρ,τ = 0.162 GeV 2 and g ρ,ee = 0.171 GeV 2 , which corresponds to f ρ,τ = 209 MeV and f ρ,ee = 221 MeV in full agreement with [93]. The difference between these results can be explained by the approximaty of the relation (B.27). So we use g ρ,τ value as more directly measured one. The results of our analysis agrees with f ρ value in [53] (within about 10%), but differ from the value adopted in [51] by ∼ 25%.

C.2 Meson form factors of decay into pseudoscalar meson
To describe the semileptonic decays of the pseudoscalar meson into another pseudoscalar meson one should know the form factors f + (q 2 ), f 0 (q 2 ), f − (q 2 ) defined by Eq. (A.13), only two of which are independent. We use f + (q 2 ), f 0 (q 2 ) pair for the decay parametrization.
In turn, there are many different parametrizations of meson form factors. One popular parametrization is the Bourrely-Caprini-Lellouch (BCL) parametrization [102] that takes into account the analytic properties of form factors (see e.g. [103,104]), where the function z(q 2 ) is defined via The choice of t 0 and of the pole mass M pole varies from group to group that performs the analysis. In this work we follow FLAG collaboration [104] and take The coefficients a + n and a 0 n are then fitted to the experimental data or lattice results.

C.2.2 D meson form factors
In the recent paper [107] the form factors for D → K and D → π transitions are given in the form where z 0 = z(0). The best fit parameter values are given in Table 11.  Table 11. Best fit parameters for the form factors (C.21) of D → π and D → K transitions [107].

C.2.3 B meson form factors
Most of B meson form factors are available in literature in the form (C.16), their best fit parameter values are given in Table 12. The form factors for B s → D s are almost the same as for B → D transition [109], so we use the same expressions for both cases.

C.3 Meson form factors for decay into vector meson
One of the relevant HNL production channel is pseudoscalar meson decay h P → h V α N . To compute decay width of this decay one needs to know the form factors g(q 2 ), f (q 2 ), a ± (q 2 ), defined by Eqs. (A.15, A.16 [110][111][112].   [110][111][112]. Masses of B c , D s and D * s are taken from [63], while for B * c theoretical prediction [113] is used.
combinations are introduced as For these linear combinations the following ansatz is used

(C.29)
A hh 1/2 (q 2 ) = f hh (C.30) Best fit values of parameters are adopted from papers [110][111][112]. f , σ parameters are given in Table 13, while ξ and the pole masses M V and M P are given in Table 14.
D Production from J/ψ and Υ mesons D.1 Production from J/ψ The process J/ψ → Nν allows to creates HNLs with mass up to M J/ψ 3.1 GeV and therefore contribute to the production above the D-meson threshold.

E Production of heavy flavour at SHiP
For a particular application of the obtained results we revise the HNL production at the SHiP experiment. The number of mesons produced by E p = 400 GeV proton beam at the SHiP target can be estimated as where X qq represents the qq production rate, f (h) is the meson h production fraction 16 and expected number of protons on target is N P oT = 2 · 10 20 . The following cross sections have been used for the estimates: • the proton-nucleon cross section is σ(pN ) 10.7 mbarn.
• σ(cc) = 18 µbarn [115] and the fraction X cc = 1.7 × 10 −3 • σ(bb) = 1.7 nbarn [116] and the fraction X bb = 1.6 × 10 −7 To calculate the meson production fractions the dedicated simulation is needed. It should take into account the properties of the target (materials, geometry) and the cascade processes (birth of the excited meson states like D * and its decay into D). The values of f (h) for the case of SHiP were calculated in the paper [57]. These values with the number of different mesons are presented in Table 15. For kaons we do not divide them for species. Taking into account production fractions of different mesons the main production channels from charm and beauty quarks for SHiP are shown in Fig. 20.

F Vector-dominance model
Here we provide F π (s) formula, given by vector-dominance model [85]  and h (s) is a derivative of h(s).