Heavy Neutral Leptons from kaon decays in the SHiP experiment

We calculate the signal rate of hypothetical heavy neutral leptons (HNL or sterile neutrinos) from kaon decays expected in the framework of the SHiP experiment. The kaons are produced in the hadronic shower initiated in the beam-dump mode by 400\,GeV protons from CERN SPS. For a sufficiently light HNL (when the decays are kinematically allowed) we find kaon decays to be a noticeably richer source of HNL as compared to $D$-meson decays adopted in previous studies of the HNL phenomenology at SHiP. In particular, SHiP is capable of fully exploring the central part of the kinematically allowed region of the HNL mass and mixing with active neutrinos down to the lower cosmological bound. The latter is associated with HNL decays in the early Universe to energetic products rescattering off and thus destroying light nuclei produced at the primordial nucleosynthesis. A consistency of the HNL model with smaller mixing would require either a non-standard cosmology or new ingredients in the HNL sector, closing the room for the minimal seesaw type I model with sterile neutrinos lighter than kaons.

1. The experiment SHiP (Search for Hidden Particles) [1,2] on the base of the CERN SPS proton beam has been proposed [3] with primary goal of searching for hypothetical feebly interacting light particles. Many well-motivated examples can be found in Ref. [4] highlighting the SHiP physics case. The SPS protons of 400 GeV hit the Molybdenum-Tungsten target, producing secondary particles, which subsequently rescatter and decay forming hadronic and electromagnetic showers degrading downstream in the target. The interesting hypothetical particles can emerge in the target and, being feebly interacting, propagate freely. On the contrary, the host of ordinary, Standard Model (SM) particles must be either absorbed in the downstream dumper or deflected (mostly muons) by specially designed magnets. The SHiP target complex is optimized to obtain a flux of hypothetical particles (if any) with as little contamination by the ordinary particles as possible. The hypothetical particles then may show up downstream in the detector system by decaying into SM particles, which make the signal events.
While this scheme is fairly universal giving a good chance to explore many extensions of the SM with feebly interacting light particles, it has been chosen in [3] based on the seemingly best performance in direct searches for sterile neutrinos -massive fermions, singlet with respect to the SM gauge group, which mix with active neutrinos -dubbed heavy neutral leptons (HNL). These particles, if sufficiently light, can be produced via mixing mostly in decays of mesons and similarly decay into lighter SM particles. If the mixing is small, a beam-dump facility is the best place to search for these particles [5]. The heavier the meson, the bigger the room for HNL mass to be investigated. However, presently available in a beam-dump mode proton beams are not sufficiently energetic to produce beauty mesons in the required amount, and most promising candidates to be the HNL sources in the beamdump experiment are D-mesons. Based on this idea a new dedicated experiment at CERN has been proposed [6], and its structure, detector geometry and position were chosen [3] to get the highest signal rate of HNL from D-mesons produced in a fixed-target experiment by protons at CERN SPS. The estimate of this signal was further revised in [4] and later in [7], motivated by both the project natural development and theoretical progress [8].
While D-mesons can produce HNLs with masses up to about 2 GeV, the HNL of smaller masses, up to about 0.5 GeV and 0.14 GeV can be also produced in decays of kaons and pions, respectively. So far these production modes have been fully neglected. The reason is natural: in the target material the kaon and pion decay lengths grossly exceed their interaction lengths. Consequently, these light mesons get absorbed and certainly degrade in energy. The 3-momenta of produced HNL are generically not pointed towards the detector, and their typical values are too small to pass the selection criteria for signal events and to suppress the background. Nevertheless, this question has not been studied before, and here we present our estimates for the signal rates of HNL appeared in kaon decays inside the SHiP target. We observe, that disadvantages of the mildly-relativistic kinematics are compensated to large extent by a high kaon production rate. We find this result encouraging, since the expected signal rate at the SHiP detector is not only noticeably higher than that for HNL (of the same mass) from D-mesons, but also high enough to fully explore a large part of the allowed by the Big Bang Nucleosynthesis region. Inside it the HNL decays in the early Universe do not spoil the production of light nuclei in the cosmic plasma (the energetic SM particles from HNL decays might be dangerous, destroying the light nuclei). This finding strengthens the physics case of SHiP, which results may be of great importance for the seesaw type I mechanism [9,10] of generating active neutrino masses.
The paper is organized as follows. We start with the model description explaining possible origin of the HNL mixing with active neutrinos and introduce the required notations. Then we outline the part of SHiP facility relevant for our study: the target structure, the position and geometry of the main detector. Further we turn to the kaon production by the beam protons on target, that we calculate by making use of the GEANT4 package [11]. Then we trace kaon decays into HNL and the latter subsequent decays inside the detector volume to the SM particles forming the signal events. Applying the cuts, presently accepted by the SHiP collaboration, as well as the zero-background hypothesis, we calculate the number of signal events and the corresponding limits on the sterile-active neutrino mixing to be placed by the SHiP experiment in case of no evidence. Finally, we summarize with discussion of possible options to improve the SHiP performance in searches for HNL.
2. In this study we follow a pure phenomenological approach, which guarantees a universal result, relevant for various models with HNL. Namely, to the Standard Model (SM) we introduce one hypothetical fermion of a mass M N smaller than the kaon mass, M N < M K . The fermion is a singlet with respect to the SM gauge group, but it mixes with active neutrinos (electron ν e , muon ν µ , and τ -neutrino ν τ ), that we parameterise with dimensionless variables U α , α = e , µ , τ . Considering the models with |U α | 1 and M N much exceeding the scale of the active neutrino mass, one can treat this fermion as the fourth massive neutrino (that explains the notation in mixing parameter and the very term HNL), but not as the neutrino of a hypothetical fourth SM generation (that explains the alternative name sterile neutrino).
From the phenomenological side we indeed have in the model a heavy neutrino, which participates in all weak processes (if kinematically allowed) but with effective coupling of the weak gauge constant multiplied by the corresponding parameter U α . This determines both the HNL production and decay. In the absence of any other interactions this fermion is an example of the feebly interacting particles.
From the theoretical side, this construction can be found in many extensions of the SM. Apparently, the simplest yet self-contained example is provided by the sterile neutrino model in the framework of the see-saw type I mechanism [9,10]. There the NHLs N I , I = 1, 2, . . . are Majorana fermions with the Lagrangian Here Y αI and M I are dimensionless Yukawa couplings and Majorana masses, L α are SM lepton doublets, andH = × H * is the conjugated SM Higgs doublet (the weak indices are omitted). The Higgs field gains a non-zero vacuum expectation value, which yields a mixing mass term for N I and active neutrinos ν α . Upon diagonalization the active neutrinos acquire masses (that explains neutrino oscillations), and HNL mass states mix with active neutrinos. If all M I 1 eV, the mixing is small and original N I are almost identical to the HNL mass states. At the same condition, the active neutrino mass scale is suppressed by mixing squared (the see-saw mechanism at work), which may explain why this scale is much lower than any other mass of the SM particles.
At least three sterile neutrinos are needed to provide masses to all three active neutrino states. The sterile neutrino masses are free parameters in between 1 eV and 10 14 GeV (the latter scale refers to the perturbativity limit of |Y αI | < 1), for details see e.g. [12]. If some of the sterile states are lighter than the kaon (but significantly heavy to decay into electrically charged SM particles), the corresponding phenomenology at SHiP is exactly the same as we consider in this paper. A particular model can predict not one, but two, three, etc HNL states in the interesting mass range. Then the mixing parameters U αI may be appropriately constrained and related. We do not study this case, but our results can be straightforwardly modified to cover it. In a similar way, one can refine our limits for the case with degenerate sterile neutrinos, where the detector fails to distinguish two (three, etc) close in mass states. This situation is realized in νMSM (for a review see [13]), where two sterile neutrinos are highly degenerated to explain both phenomena: neutrino oscillations and the baryon asymmetry of the Universe. The mixing pattern, |U αI |, and the corresponding signals are specifically constrained in this case [5].
In what follows we consider the model with one HNL of mass M N M K and small and unrelated mixing parameters |U α | 1, which for a realistic models are additionally bounded from above by direct searches for HNL, for a recent summary see e.g. [14].
3. In the SHiP experiment a 400 GeV proton beam hits a target, where secondary particles form hadronic showers dominated by pions and kaons. When a kaon decays it has a chance (proportional to the corresponding differential branching ratio) to produce a sterile neutrino of a given 3-momentum. For a heavy neutrino two-body kaon decays are not suppressed by the chirality flip, so they dominate over three-body decays because of the phase volume, see formulas and plots in Ref. [5]. If HNL is heavier than a pion the domination factor is at least 50 and typically much higher. Moreover, the three-body decay is relevant only for more narrow range of HNL masses, M N < M K − M π , due to kinematics. Therefore, we account only for charged kaons and consider only their two-body decays as the HNL sources.
The HNL lives long and interacts feebly, hence it propagates freely along its 3-momentum at production. To be observed at SHiP, its trajectory must pass through the decay volume, which is placed at a distance of L = 64 m from the target. The decay volume is of ∆l = 50 m length and has a conical form along the beam with a rectangular base [7]:  detector almost coincides with its decay volume.
To be observed, the HNL must decay inside the detector volume. The background-free signatures assume (at least) two charged particles in the final state, each with 3-momenta of ≥1 GeV/c. Both particles must cross the far end of the decay volume to be registered by the SHIP decay spectrometer. The detection efficiency is expected to be high, so we set it to 100% in our study.
4. To obtain the distribution of the HNL decay products over 3-momenta, we start with the analysis of the kaon production and decay in the SHiP target. We simulate the kaon production by 400 GeV protons by exploiting the GEANT4 framework [11]. The target design is taken from [1]. The schematic view of the Molybdenum-Tungsten target with water interlayers is shown in Fig. 1. The emerged in proton interactions light mesons propagate farther, mostly scattering in the target material and thus producing more mesons. The physics list we utilize for the processes in the target is called "QGSP BERT". It uses the quark-gluon string model for the high energy events (energy in the center-of-mass frame is 10 GeV) and the Bertini cascade model for lower energies ( 10 GeV). All the standard EM processes and decays are accounted for.
Naturally, the most energetic kaons are produced in the very first collision of the incident protons. The spectrum of these charged kaons K + is presented in Fig. 2. The secondary hadrons (and degraded protons) also produce kaons inside the target. The distribution of the charged kaon production points is shown in Fig. 3. One observes that most of the kaons are produced in the Molybdenum layers in the upstream target region. Since secondary particles degrade in energy while the hadronic shower develops, the produced by them kaons have lower energy. The spectrum of all K + at their production points inside the target is shown in Fig. 4.
After production, each koan propagates in the target, and having a strange quark in its content, can either decay via weak coupling or interact with a nuclei. In case of inelastic interactions, some charged kaons may transform to either the neutral ones or strange baryons. The former may transform via inelastic scattering back to charged kaons or decay into nonstrange particles, while the latter decay into non-strange baryons and mesons. Both neutral kaons and strange baryons can produce HNLs via semileptonic weak decays, but those are three-body decays, so the rates are significantly suppressed, and overall kinematics makes Positive kaon production kinematics Kaon production point them open only for lighter HNLs. Therefore, only charged kaons produced in the target are considered as parents of potentially interesting HNLs.
The kaon can decay at any moment, but in the laboratory frame this chance grows with kaon attenuation: most kaons decay (almost) at rest, when its 3-momentum falls below about Positive kaon production kinematics 0.1 MeV/c. However, in this case the decay products have no preferable directions, and only a very small part of them pass through the SHiP decay volume. Moreover, the kaon, decaying at rest, cannot produce HNL with the total energy >0.5 GeV, that subsequently constrains the HNL decay product energy, well below the threshold of the SHiP decay spectrometer. Therefore, we neglect these stopped kaons as well.
To conclude, the only relevant for HNL searches are charged kaons decaying in flight. Based on the simulated statistics of N stat = 1.5 × 10 6 POT we observe that each proton produces about 8.01 K + and 3.51 K − , out of which 4.52 K + and 2.85 K − are captured (while strangeness transfer mostly to neutral kaons and baryons, respectively) and 3.18 K + and 0.59 K − decay at rest. Therefore, the only profitable are 0.29 K + and 0.07 K − decaying in flight. Their distribution over total 3-momentum and polar angle (with zenith direction along the beam line) is presented in Fig. 5. Let us note in passing that about 0.1% of the produced charged kaons leave the target. We neglect their contribution. Likewise, in our simulations we ignore a possible impact of the active muon shield with magnetic field (designed primary to deflect the muons outside the target) on the propagation of charged kaons inside the target.
5. The described above GEANT4 simulation yields a set of coordinates {x 1 , y 1 , z 1 } and 3-momentum p K = {p x 1 , p y 1 , p z 1 } of the charged kaons at the moment of their decay. We take into account both K + and K − . Since each proton in the beam produces on average M p→K = 0.36 charged kaons that decay in flight, the total number of in-flight kaon decays is expected to be K tot = N P OT × M p→K .
Each charged kaon may decay into HNL with the probability equal to the corresponding branching ratio (recall, we account only for the two-body kaon decays) [5] Br (1) Here f K = 155.6 MeV and τ K = 1.238 × 10 −8 s are the kaon decay constant and lifetime, and M l is the charged lepton mass, l = e, µ. The total number of HNL produced by the relevant kaons at SHiP reads To use all the GEANT4 statistics, in our numerical calculations each decaying in flight charged kaon is assumed to produce a sterile neutrino. In the final calculation we account for that via multiplying the number of "signal events" by N N /K tot . Unless the statistics is too low, we obtain viable results. If our statistics proves to be insufficient, we improve it by calculating several (tens, hundreds) sterile neutrino production scenarios for each event of the kaon decay simulation. The HNL trajectories and chances to decay inside the SHiP decay volume are obtained following the procedure of Ref. [15]. The trajectory starts at the point of the kaon decay, {x 1 , y 1 , z 1 }, provided by the simulations. In the kaon rest frame the HNL 3-momentum p = {p x , p y , p z } is randomly distributed over the dimension-2 sphere of radius fixed for a given decay mode (1) as To obtain the NHL 3-momentum in the laboratory frame p N , we make a Lorentz boost defined by the kaon 3-momentum in the laboratory frame taken from the simulations. We choose the point where proton beam crosses the target as the origin of the coordinate system, z axis is directed towards the detector and x, y axes are chosen so that for a particular decaying kaon its longitudinal component is along z-axis, while the transverse component is along x: p K L ≡ p Kz = p z 1 , p K T ≡ p Kx , p Ky = 0. Then the components of HNL 3-momentum in this frame read: p Ny = p y , This frame is related to the laboratory frame (where the kaon 3-momentum equals p K = {p x 1 , p y 1 , p z 1 }) through rotation along the z-axis which restores a non-zero y-component of the kaon momentum. Hence, to adopt formulas (3)-(5), we perform two rotations: before the Lorentz-boost we rotate to the frame there p Ky = 0, and after the boost we rotate back to the laboratory frame. We plot the overall spectra of the obtained in this way HNL in Fig. 6. In what follows we keep for the HNL 3-momentum in the laboratory frame the same notations as in (3)-(5), that is p N = {p Nx , p Ny , p Nz }. 6. Next, we check, whether the HNL trajectory crosses the decay volume or not. The interesting HNL must have p Nz > 0, otherwise it flies backwards. Then, at the moment when its position is at z N = L (where the decay volume begins, see Tab. 1), its other coordinates are If HNL decays in the decay volume its coordinates are: where u ∈ (0, 1) is some number. Inequalities guarantee that the decay happens inside the decay volume. We count all neutrinos that enter the decay volume along the suitable trajectory and call the resulting number N detector . We present the ratios of N detector /N N of sterile neutrinos reaching the decay volume on the left panel of Fig. 7. The ratios grow with the HNL mass and finally reach about 1%. They illustrate the signal suppression due to the geometrical factor. The HNL spectra are shown on right panel of Fig. 7.
The probability for the sterile neutrino to decay at the moment it reaches coordinates (7) -(9) for any u (that means it decays inside the decay volume) reads the limit is valid for our interesting case of small mixing, when the HNL decay length exceeds the length of the decay volume (and most of HNLs safely leave the detector volume); the most relevant quantity here is the HNL life-time τ N . To obtain the number N decay of HNLs which decay inside the decay volume we count all of the HNL which reach the detector N detector with a weight P from eq. (10). We present the ratios of N decay /N detector in Fig. 8 (left panel), setting |U α | 2 = 10 −8 as a reference number in the HNL lifetime τ N , which is The sterile neutrino lifetime (right panel). We consider mixing with either electron or muon neutrinos, calculated as explained below. 7. HNLs mix with active neutrinos, and hence decay into the SM particles due to weak interactions. The formulas for HNL decay rates can be found in Ref. [5]. The branching ratios relevant for the interesting mass range are presented in Fig. 9. One observes, that in HNL decays as well as in their production, the two-body decay modes, when kinematically open, dominate. Moreover, the 2-body mode into charged particles dominates over that into neutral particles (except the threshold regions). Thus, for the signal events we consider mostly only these charged modes (π ± l ∓ , where l is electron and muon), also because they are treated presently at SHiP as the background-free ones. When 2-body decays are kinematically forbidden, we consider pure leptonic modes with two charged leptons, N → l + α l − β ν. They are also treated as background free at SHiP [7].
The HNL lifetime for the reference mixing |U α4 | 2 = 10 −8 is shown on the right panel of Fig. 8. One finds that the HNL decay length decreases with its mass, but for the reference level of mixing it always vastly exceeds length of the SHiP decay volume, which justifies the limit (10).
To be observed, the charged particles must cross the far end of the decay volume to enter the decay spectrometer installed there, which allows to identify the particle type and to measure its energy. Moreover, its presently accepted operation scheme is supposed to be almost 100% effective, if at least two charged particles from the decay have 3-momenta exceeding 1 GeV/c. These constraints on the 3-momenta of the HNL decay products we further implement in our study applying the procedure above we use for the HNL trajectories. We briefly describe it here for the 2-body decay modes. Namely, we replace N with π and K with N in (2) and (5), (4). In (6) at the moment z π = L + ∆l we replace starting point {x 1 , y 1 , z 1 } with {x N , y N , z N } from (7) - (9). In this way we obtain the pion coordinates at the moment it reaches the end of the decay volume {x π , y π , z π }. Taking into account that p N = p π + p l , similarly we obtain {x l , y l , z l }, z l = L + ∆l. Both pion and charged lepton cross the backside of the decay volume if |x π | < ∆w 2 , |x l | < ∆w 2 and |y π | < ∆h 2 , |y l | < ∆h 2 . In the same way as we do above with N detector , we count all sterile neutrinos that pass these geometrical criteria and call the resulting number N observed . Thus we find the geometrical acceptance of the detector det ≡ N observed /N decay for the HNLs decaying in the decay volume. We present the geometrical acceptance of the decaying sterile neutrinos on the left panel of Fig. 10. They only mildly depend on the HNL mass. The spectra of those HNLs, which decay products fulfill the geometrical criteria in the SHiP experiments are given on the right panel of Fig. 10. Finally, we apply the cut on the 3-momenta of the charged particles the HNL decay to, π ± and l ∓ , obeying the geometrical criteria. This is rather constraining criterion, that we illustrate on the left panel of Fig. 11 for the two chosen values of 0.5 GeV/c and in the signal statistics.
We arrive at the similar conclusions for the leptonic decay modes. 8. To obtain the number N f inal of the observed sterile neutrinos that decay inside the decay volume we count all of the observed sterile neutrinos N observed with a weight P from eq. (10). We present the final number of the observed sterile neutrinos N f inal in Fig. 12 for 2-body decay channels. Accepting the background-free scheme at SHiP, we ask N f inal ≤ 2.3 to put the 90% CL upper limit on the active-sterile mixing at a given HNL mass, assuming no evidence for the chosen HNL signature. Recall, that we accepted no suppression, i.e. |U α | 2 = 1, for the productions and use |U α | 2 = 10 −8 for HNL decays as the reference number to make numerical estimates and present the results on plots. Since the number of events depends on the product of our two reference numbers, the correct limit in our scheme is given as where α and β refer to the mixings responsible for HNL production and decay, respectively. Moreover, the limits are placed investigating a particular signature -decay mode -and we consider them separately. The corresponding SHiP sensitivity is presented in Figs.13, 14, 15 and 16 together with SHiP sensitivity based on D-meson decays [7] and previous limits Figure 13: Present limits on the HNL mixing with electron neutrino and expected sensitivity of the SHIP experiment. Also shown are the results obtained in PS191 [16] and T2K [17] (within Feldman-Cousins approach).
on the sterile-active mixing from direct searches, from cosmology (grey regions: Big Bang Nucleosynthesis constraints are scanned from Fig. 4.10 of Ref. [4] which in the interesting mass range is repeated from the restriction of τ s < 0.1 sec in Fig. 1 of Ref. [13]) and the reference line for the see-saw mechanism relation, U 2 ≈ √ m 2 atm M N ≈ 0.05 eV/M N . One observes that the SHiP experiment is capable of fully exploring the central parts of kinematically available regions above the cosmological BBN constraint (regions in grey). Remarkably, the see-saw line refers approximately to the minimal mixing at which the corresponding sterile neutrino contributes sufficiently large amount to active neutrino masses. The region below BBN line is forbidden from cosmology. Therefore, if the direct searches reveal negative results, this will close the simple seesaw scheme for the light HNL provided BBN limits are solid.
Note, that in our study we treat HNL as Majorana particles, while they are considered to be Dirac in PS191 and in E949 experimental analyses. Therefore, to present these results on our plots we have divided their limits by square root of two. Moreover, the original analysis of PS191 collaboration is based only on HNL interactions via charged currents: the neutral current contributions to all processes with HNL have been ignored. Following the analysis of Ref. [19] we include the missed contributions from the neutral currents and implement all the relevant corrections for PS191 limits presented in Figs. 13-15. Figure 14: Present limits on the HNL mixing with muon neutrino and expected sensitivity of the SHIP experiment. Also shown are the results obtained in PS191 [16], E949 [18], and T2K [17]. Figure 15: Present limits on the combination |U e U µ | (left panel) and |U µ U e | (right panel) of HNL mixing with muon and electron neutrinos and expected sensitivity of the SHIP experiment. The parameter products are the same, but the corresponding observables are different, so we present them separately for convenience.

9.
To conclude, in this paper we calculate the signal events from hypothetical heavy neutral leptons emerged from decays of charged kaons inside the Mo-Tn target of the SHiP experiment. We find this source of HNL is more promising for masses in the range of 150 - 500 MeV than previously considered D-mesons [7], so that it allows SHiP to fully explore the large part of phenomenologically and cosmologically viable region of sterile-active mixing in models with HNLs lighter than kaons. In this study we account for only the main production channels and HNL decay modes. It seems that corrections from the subdominant processes do not change much. First, consider decays. The dominant decay modes are 2-body decays, which rates are amplified by the chirality flip factor ∝ M 2 N . The 3-body decays are suppressed by the volume phase factor and are kinematically forbidden for a heavy HNL. The low-threshold three-body leptonic mode, N → e + e − ν becomes competitive only for a light HNL, M N m π , see Fig. 9, but then pions must be introduced to the analysis as a possible source of such HNL, that is beyond the scope of this paper. Second, turn to the production. There, the strange baryons weakly decaying to ordinary ones would also contribute [20] but only to light HNL of M N 100 MeV, where again, the pions must be considered. Notice, neutral kaons may decay into a heavier HNL (e.g. K 0 → πeN ), but with strongly suppressed rates, see Ref. [5].
The signal rate of HNL with m π < M N < M K can be increased by a factor of about 2 with a lower cut (say, 0.5 GeV/c) on the 3-momenta of the HNL decay products, see the left panel of Fig. 11. Further, in that case one can also investigate the decay events when one of the charged particle goes backwards (upstream, towards the target) and use the tracker system and the neutrino detector in front of the SHiP decay volume to identify the particle type and measure its momentum. An ultimate (but probably unrealistic) option is that sufficiently low cuts would allow SHiP to register the HNL from kaon decays at rest, which statistics exceeds that of decays in flight by an order of magnitude. The detector geometric acceptance, see left panel of Fig. 10, can be possibly increased by using the veto system on the surface of the decay volume to help in detecting the HNL decay event: time and position of crossing by the outgoing charged particle.
The overall SHiP sensitivity to models with HNL can possibly be increased in this way, but modestly. The number of signal events scales as |U α | 2 × |U β | 2 , and these steps can help to pull down the expected limits by a factor of 2, farther probing the region around the cosmological BBN constraint. The next major scientific goal of direct searches -testing the region of the seesaw line |U α | 2 ∼ m ν /M N , see Figs. 14-15 -requires a factor of hundred increase in the number of signal events and will be a great challenge for the next generation experiments.
The work was supported in part by the RSF Grant No. 17-12-01547.