Sensitivity of the SHiP experiment to Heavy Neutral Leptons

Heavy Neutral Leptons (HNLs) are hypothetical particles predicted by many extensions of the Standard Model. These particles can, among other things, explain the origin of neutrino masses, generate the observed matter-antimatter asymmetry in the Universe and provide a dark matter candidate. The SHiP experiment will be able to search for HNLs produced in decays of heavy mesons and travelling distances ranging between $\mathcal{O}(50\text{ m})$ and tens of kilometers before decaying. We present the sensitivity of the SHiP experiment to a number of HNL's benchmark models and provide a way to calculate the SHiP's sensitivity to HNLs for arbitrary patterns of flavour mixings. The corresponding tools and data files are also made publicly available.

Abstract: Heavy Neutral Leptons (HNLs) are hypothetical particles predicted by many extensions of the Standard Model. These particles can, among other things, explain the origin of neutrino masses, generate the observed matter-antimatter asymmetry in the Universe and provide a dark matter candidate.
The SHiP experiment will be able to search for HNLs produced in decays of heavy mesons and travelling distances ranging between O(50 m) and tens of kilometers before decaying. We present the sensitivity of the SHiP experiment to a number of HNL's benchmark models and provide a way to calculate the SHiP's sensitivity to HNLs for arbitrary patterns of flavour mixings. The corresponding tools and data files are also made publicly available.  The SHiP experiment. The Search for Hidden Particles (SHiP) experiment [1][2][3][4] is a new general purpose fixed target facility proposed at the CERN Super Proton Synchrotron (SPS) accelerator to search for long-lived exotic particles with masses between few hundred MeV and few GeV. These particles are expected to be predominantly produced in the decays of heavy hadrons. The facility is therefore designed to maximise the production and detector acceptance of charm and beauty mesons, while providing the cleanest possible environment. The 400 GeV proton beam extracted from the SPS will be dumped on a high density target with the aim of accumulating 2 × 10 20 protons on target during 5 years of operation. The charm production at SHiP exceeds that of any existing and planned facility. A dedicated detector, based on a long vacuum tank followed by a spectrometer and by particle identification detectors, will allow probing a variety of models with light long-lived exotic particles. Since particles originating in charm and beauty meson decays are produced with a significant transverse momentum with respect to the beam axis, the detector should be placed as close as possible to the target. A critical component of SHiP is therefore the muon shield [5], which deflects away from the detector the high flux of muons produced in the target, that would otherwise represent a very serious background for hidden particle searches. To suppress the background from neutrinos interacting in the fiducial volume, the decay volume is maintained under vacuum [3]. The detector is designed to reconstruct the exclusive decays of hidden particles and to reduce the background to less than 0.1 events in the sample of 2×10 20 protons on target [4]. The detector consists of a large magnetic spectrometer located downstream of a 50 m long and 5 × 10 m wide decay volume. The spectrometer is designed to accurately reconstruct the decay vertex, mass and impact parameter of the decaying particle with respect to the target. A set of calorimeters followed by muon chambers provide identification of electrons, photons, muons and charged hadrons. A dedicated timing detector measures the coincidence of the decay products, which allows the rejection of combinatorial background.
The decay volume is surrounded by background taggers to tag neutrino and muon inelastic scattering in the surrounding structures, which may produce longlived neutral Standard Model particles, such as K L , that have similar topologies to the expected signal.
The experimental facility is also ideally suited for studying the interactions of tau neutrinos. It will therefore host an emulsion cloud chamber based on the Opera concept, upstream of the hidden particle decay volume, followed by a muon spectrometer. The SHiP facility layout is shown in Fig. 1. Recent progress report [4] outlines the up-to-date experimental design as well as describes changes since the initial technical proposal [2].
Heavy Neutral Leptons. Among hypothetical long-lived particles that can be probed by the SHiP experiment are Heavy Neutral Leptons (or HNLs) [6]. The idea that HNLs -also known as right-handed, Majorana or sterile neutrinos -can be responsible for the smallness of neutrino masses goes back to the 1970s [7][8][9][10][11][12]. It has subsequently been understood that the same particles could be responsible for the generation of the matter-antimatter asymmetry of the Universe [13]. The idea of this scenario, called leptogenesis, was developed since the 1980s (see reviews [14][15][16][17][18][19] and references therein). In particular, it was found that the Majorana mass scale of pN cross-sectioncc fractionbb fraction Cascade enhancement f cascade σ pN [2] Xc c [62] Xb b [63] charm [64] beauty [64] 10.7 mb 1.7 × 10 −3 1.6 × 10 −7 2.3 1.7 Table 1. Charm and beauty production fractions and cascade enhancement factors for the SHiP experiment. Cross-section σ pN is an average proton-nucleon inelastic cross-section for the molybdenum target [2].
right-handed neutrinos can be as low as O(GeV) [20][21][22], thus providing a possibility for a leptogenesis scenario to be probed at a particle physics laboratory in the near future.
It was demonstrated in 2005 that by adding just three HNLs to the Standard Model one could not only explain neutrino oscillations and the origin of the baryon asymmetry of the Universe, but also provide a dark matter candidate [21,23]. Two of the HNLs should have masses in the GeV range, see [24] for a review. This model, dubbed Neutrino Minimal Standard Model (or νMSM), is compatible with all the measurements so far performed by accelerator experiments and at the same time provides a solution for the puzzles of modern physics [24,25]. This made models with GeV scale HNLs a subject of intensive theoretical studies in the recent years [19,[26][27][28][29][30][31][32][33][34][35][36][37][38][39][40][41][42][43][44][45].
HNLs are massive Majorana particles that possess neutrino-like interactions with W and Z bosons (the interaction with the Higgs boson does not play a role in our analysis and will be ignored). The interaction strength is suppressed compared to that of ordinary neutrinos by flavour dependent mixing angles U α 1 (α = {e, µ, τ }). Thus, even the simplest HNL model contains 4 parameters: the HNL mass M N and 3 mixing angles U 2 α . 1 The idea of experimental searches for such particles goes back to the 1980s (see e.g. [46][47][48][49][50]) and a large number of experiments have searched for them in the past (see review of the past searches in [51][52][53]). HNLs are being searched at currently running experiments, including LHCb, CMS, ATLAS, T2K, Belle and NA62 [54][55][56][57][58][59][60][61].
The sensitivity of the SHiP experiment to HNLs was previously explored for several benchmark models [2,65,66] assuming particular ratios between the three HNL mixing angles [51]. This paper updates the previous results in a number of important ways. A recent work [67] revised the branching ratios of HNL production and decay channels. In addition, the estimates of the numbers of D-and B-mesons now include cascade production [64]. We update the lower limit of the SHiP sensitivity region and also evaluate the upper bound for the first time. We discuss potential impact of HNL production from B c mesons. Moreover, our current sensitivity estimates are  Figure 2. HNL production branching ratios multiplied with the production fraction of the meson decaying into HNL, for charm (left) and beauty (right) mesons [67]. The mixing angles have been set to U 2 e = 1, U 2 µ = U 2 τ = 0. The production from D + and B + remains relevant for higher masses for D 0 and B 0 because of the fully leptonic decays h + → N + + . The B c production fraction is unknown (see text for details) and we show two examples: not limited to a set of benchmark models. Rather, we compute a sensitivity matrix -a model-independent tool to calculate the SHiP sensitivity for any model of HNL flavour mixings.
The paper is organised as follows. Section 2 describes the simulation of HNL events. The resulting sensitivity curves for mixing with each individual flavour, for the benchmark models of Ref. [2] as well as the sensitivity matrix -are discussed in Section 3. We present our method to evaluate the SHiP sensitivity to HNLs in a model-independent way in Section 4 and conclude in Section 5.

Monte Carlo simulation of heavy neutral leptons at SHiP
A detailed Monte Carlo simulation suite for the SHiP experiment, FairShip, was developed based on the FairRoot software framework [69]. In FairShip simulations primary collisions of protons are generated with Pythia 8 [70] and the subsequent propagation and interactions of particles simulated with GEANT4 [71]. Neutrino interactions are simulated with GENIE [72]; heavy flavour production and inelastic muon interactions with Pythia 6 [73] and GEANT4. Secondary heavy flavour production in cascade interactions of hadrons originated by the initial proton collision [64] is also taken into account, which leads to an increase of the overall HNL production fraction (see Table 1). The SHiP detector response is simulated using GEANT4. The pattern recognition algorithms applied to the hits on the straw spectrometer are described in [74], and the algorithms for particle identification are presented in [75].
The simulation takes the HNL mass M N and its three flavour mixings U 2 e , U 2 µ , U 2 τ as input parameters. For the pure HNLs mixing to a single SM flavour, the number Table 2. Production fraction and expected number of different mesons in SHiP taking into account cascade production [68]. For f (b → B c ) see text for details.
of detected HNL events N events is estimated as 2 where N prod is the number of produced HNLs that fly in the direction of the fiducial volume and P det is the probability of HNL detection in the Hidden Sector detector. The number of produced HNLs is is the h meson production fraction 3 at SHiP (see Table 2), BR(h → N + X) is the mass dependent inclusive branching ratios for h mesons decays with HNL in the final state and decay is the geometrical acceptance -the fraction of produced HNLs that fly into direction of the fiducial volume. Fig. 2 shows the product between the meson production fraction and its inclusive decay branching fraction into sterile neutrinos. Finally, N q is the total number of produced quarks and antiquarks of the given flavour q taking into account the quark-antiquark production fraction Xq q and the cascade enhancement factor f cascade given in Table 1, The HNL detection probability is given by where BR(N → visible) is the total HNL decay branching ratio into visible channels (see HNL decay channels in Appendix A), P decay is the probability that the HNL decays inside the fiducial volume, where l ini is the distance travelled by HNL before it entered the decay vessel; l fin is the distance to the end of the decay vessel along the HNL trajectory; l decay = cγτ N is the HNL decay length (γ and τ N being HNL gamma factor and proper lifetime). Finally, det is the efficiency of detecting the charged daughters of the decaying HNL. It takes into account the track reconstruction efficiency and the selection efficiency, further described in [2,65,75]. In order to distinguish the signal candidates from possible SM background, we put a criteria that at least two charged tracks reconstructed to the decay point are present. The reconstruction efficiencies for the decay channels N → µµν and N → µπ are given in e.g. [2, Section 5.2.2.2]. Using FairShip, a scan was done over the HNL parameter space. For each set of HNL parameters we ran a simulation with 300 HNL events, produced randomly from decay of mesons. We determined P decay , decay and det in each of them and average over simulations to find the expected number of detected events,N events . For HNLs with masses M N 500 MeV kaon decays are the dominant production channel. While O(10 20 ) kaons are expected at SHiP, most of them are stopped in the target or hadron stopper before decaying. As a consequence, only HNLs originating from charm and beauty mesons are included in the estimation of the sensitivity. SHiP can however explore the νMSM parameter space down to the constraints given by Big Bang nucleosynthesis observations [76,77], even with this conservative assumption. It is expected that the NA62 experiment will also probe the region below the kaon mass [78].
For HNL masses M N 3 GeV the contribution of B c mesons to the HNL production can be relevant because the B + c → N + + decay width is proportional to the CKM matrix element |V cb | 2 , while the decays of B + are proportional to |V ub | 2 [51,67]. The ratio |V cb | 2 /|V ub | 2 ∼ 10 2 , which explains the relative importance of B c channels even for small production fraction f (b → B c ). This production fraction has not been measured at the SHiP center of mass energy. If the B c production fraction at SHiP is at the LHC level, its contribution will be dominant. However, at some unknown energy close to the B c mass this production fraction becomes negligible. The existing Tevatron measurement place f (b → B c ) = 2.08 +1.06 −0.95 × 10 −3 at √ s = 1.8 TeV [79]. More recent LHCb measurement at √ s = 7 and 8 TeV gave [80]. Using f (b → B + ) = 0.33 from the LHCb measurement performed at √ s = 7 TeV [81], one obtains f (b → B c ) = 2.6 × 10 −3 . Theoretical evaluations have mostly been performed for TeV energies (see e.g. [82][83][84][85]) with the exception of the works [86,87] that computed the production fraction down to energies of tens of GeV (where they found the fraction to be negligible). However, by comparing predictions of [87] with LHCb or Tevatron measurements, we see that (i) it underpredicts the value of f (b → B c ) by about an order of magnitude at these energies and (ii) it predicts stronger than observed change of the production fraction between LHC and Tevatron energies. Therefore we have to treat f (b → B c ) as an unknown parameter somewhere between its LHC value and zero and provide two estimates: an optimistic estimate for which f (b → B c ) is at the LHC level and a pessimistic estimate where we do not include B c mesons at all. In the simulation we take the angular distribution of B c mesons to be the same as that of B + mesons, based on comparisons performed with the BCVEGPY [88] and FONLL [89,90] packages, while we rescale the energy distribution according to the meson mass.
Detailed background studies have proven that the yield of background events passing the online and offline event selections is negligible [2]. Therefore, the 90% confidence region is defined as the region of the parameter space where one expects on averageN events ≥ 2.3 reconstructed HNL events, corresponding to the discovery threshold with an expected background yield of 0.1 events. Figure 3 presents the 90% C.L. sensitivity curves for HNLs mixing to only one SM flavour. The sensitivity curves have a characteristic "cigar-like shape" for masses M N > 2 GeV. The upper boundary is determined by the condition that the decay length of a produced particle becomes comparable with the distance between the target and the decay volume, and therefore the HNLs produced at the target may not reach the decay volume, see Eq. (2.5). For masses M N < 2 GeV such an upper boundary also exists, but it is outside the plot range, owing to a much larger number of parent D mesons. The lower boundary of the sensitivity region is determined by the parameters at which decays become too rare (decay length much larger than the detector size). The intersection of the upper and lower boundaries defines the maximal mass which can be probed at the experiment.

SHiP sensitivity for benchmark HNL models
We also provide updated sensitivity estimates for the three benchmark models I-III presented in the Technical Proposal [2,65]. These models allow to explain neutrino flavour oscillations while at the same time maximizing the mixing to one particular flavour, and are defined by the following ratios of flavour couplings [51]: -I.

Model independent SHiP sensitivity
In this Section we provide an efficient way to estimate the SHiP sensitivity to an HNL model with an arbitrary ratio U 2 e : U 2 µ : U 2 τ . It is based on the observation that the dependence of the number of events, N events , on the mass and mixing angles of HNL factorizes, and therefore all relevant information can be extracted from a handful of simulations, rather than from a scan over an entire 4-dimensional HNL parameter space (M N , U 2 e , U 2 µ , U 2 τ ). All information about the HNL production in a particular experiment is contained in N α (M N ) -the number of HNLs that would be produced through all possible channels with the mixings U 2 α = 1 and U 2 β =α = 0: Here N h is the number of hadrons of a given type h, BR(h→N +X α ) is the branching ratio for their decay into an HNL plus any number of other particles X α with total lepton flavour number L α = 1 and decay,α is the geometrical acceptance of HNL that in general depends not only on the mass but also on the flavour. The overall number of HNLs (given by Eq. (2.2)) produced via the mixing with the flavour α and flying in the direction of the decay vessel is given by The decay probability P decay should be treated differently, depending on the ratio of the decay length and the distance from the target to the decay vessel. It also depends on the production channel through the mean gamma factor γ α entering the decay length.
In the limit when the decay length much larger than the distance between the beam target and the exit lid of the SHiP decay volume, the U 2 β dependence of the decay probability can be accounted for similarly to Eq. (4.2): where Γ β is a decay width of the HNL of mass M N that has mixing angles U 2 β = 1, U 2 α =β = 0, the definitions of lengths l ini , l fin are given after Eq. (2.5). The index α in Eq. (4.3) indicates that the HNL was produced via mixing U 2 α (although can decay through the mixing with any flavour), so γ α is the mean gamma factor of HNLs produced through the mixing with the flavour α.
In the general case, when the decay length l decay is not necessarily larger than l fin , the analogous decay probability P decay,α can be expressed via (4.3) as follows: (4.4) where BR(N → visible) is the probability that the HNL decays into the final states detectable by SHiP. Finally, we define the HNL detection efficiency as where BR(N → X β ) is the branching ratio of a decay through the mixing angle β and det,β is the probability that the HNL decay products are successfully detected. As a result, the number of detected events is given by We see that it is sufficient to know 9 functions of the HNL mass -N α (M N ), P linear decay,α (M N ) and det,α (M N ) -to determine the number of detected events for any combination of the mixing angles.
To determine these numbers we ran 9 Monte Carlo simulations for each mass. We first ran 3 simulations with vectors − → , where x is any sufficiently small number such that l decay l det . We then ran a set of 6 non-physical simulations, where a particle is produced solely via channel α and decays solely through the channel β = α. Using results of these simulations we extract N α , P α and det,α values that allow us to generate the expected number of detected events for any values of masses and couplings.
The results are available at Zenodo platform [91] with instructions for reading the file and generating sensitivity curves at different confidence levels.

Conclusion
Using a detailed Monte Carlo simulation of HNL production in decays of charm and beauty mesons, and of the detector response to the signal generated by a decaying HNL, we calculated the sensitivity of the SHiP experiment to HNLs, updating the results presented in the Technical Proposal [2]. In particular, we assess the potential impact of HNL production from B c mesons decay, showing its influence on the extent of the probed HNL mass range. We take into account cascade production of B and D mesons as well as revised estimates of branching ratios of HNL production and decay, and we extend our calculation to masses below ∼ 500 MeV, where SHiP has a potential to fully explore the allowed region. Finally, we present our results as a publicly available dataset, providing a model-independent way to calculate the SHiP sensitivity for any pattern of HNL flavour mixings.
The SHiP experiment offers an increase of up to 3 orders of magnitude in the sensitivity to heavy neutral leptons, Fig. 5. It is capable of probing cosmologicaly interesting region of the HNL parameter space, and of potentially discovering the origin of neutrino masses and of the matter-antimatter asymmetry of the Universe.
��� ����� Figure 5. Parameter space of HNLs and potential reach of the SHiP experiment for the mixing with muon flavour. Dark gray area is excluded from previous experiments, see e.g. [6]. Black solid line is the recent bound from the CMS 13 TeV run [57]. Solid and dashed-dotted red lines indicate the uncertainty, related to the production fraction of B c mesons at SHiP energies that has not been measured experimentally or reliably calculated (see Section 2 for details). 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 HNL production from kaon decays. The primordial nucleosynthesis bounds on HNL lifetime are from [76]. 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 [6].

A HNL decays
For completeness we list the relevant HNL decay channels in Table 3 (reproduced from [67]). Table 3: List of the relevant HNL decay channels with branching ratio above 1% covering the HNL mass range up to 5 GeV implemented in FairShip. The numbers are provided for |U e | 2 = |U µ | 2 = |U τ | 2 . For neutral current channels (with neutrinos in the final state) the sum over neutrino flavours is taken, otherwise the lepton flavour is shown explicitly. Columns: (1) the HNL decay channel. (2) The HNL mass at which the channel opens. (3) The HNL mass starting from which the channel becomes relevant (branching ratio of this channel exceeds 1%). For multimeson final states we provide our best-guess estimates. (4) HNL mass above which the channel contributes less than 1%, with "-" indicating that the channel is still relevant at M N = 5 GeV. (5) The maximum branching ratio of the channel for M N < 5 GeV. (6) Reference to the appropriate formula for decay width in ref. [67].

Channel
Opens