Search for heavy neutral lepton production in $K^+$ decays

A search for heavy neutral lepton production in $K^+$ decays using a data sample collected with a minimum bias trigger by the NA62 experiment at CERN in 2015 is reported. Upper limits at the $10^{-7}$ to $10^{-6}$ level are established on the elements of the extended neutrino mixing matrix $|U_{\ell 4}|^2$ ($\ell=e,\mu$) for heavy neutral lepton mass in the range $170-448~{\rm MeV}/c^2$. This improves on the results from previous production searches in $K^+$ decays, setting more stringent limits and extending the mass range.


Introduction
Non-zero masses and mixing of the Standard Model (SM) neutrinos are now firmly established. However many SM extensions have been proposed, involving massive "sterile" neutrinos, also called heavy neutral leptons (HNLs), which mix with the ordinary light "active" neutrinos. For example, the Neutrino Minimal Standard Model (νMSM) postulates three HNLs, explaining dark matter and baryon asymmetry of the universe in a way consistent with the results of neutrino oscillation experiments [1]. One of these HNLs with the expected mass of O(10 keV/c 2 ) is a dark matter candidate, while the others are expected to have masses of O(1 GeV/c 2 ).
Mixing between HNLs (denoted N below) and active light neutrinos gives rise to HNL production in meson decays, including K + → ℓ + N (ℓ = e, µ). The branching fraction of the latter decay is determined by the HNL mass m N and mixing parameter |U ℓ4 | 2 as follows [2,3]: (1) Here B(K + → ℓ + ν) is the measured branching fraction of the SM leptonic decay (including inner bremsstrahlung), and ρ ℓ (m N ) is a kinematic factor: with x = (m ℓ /m K ) 2 , y = (m N /m K ) 2 and λ(a, b, c) = a 2 + b 2 + c 2 − 2(ab + bc + ac). By definition, ρ ℓ (0) = 1. Numerically, the product B(K + → ℓ + ν) · ρ ℓ (m N ) is O(1) over most of the allowed m N range. However it drops to zero at the kinematic limit m N = m K − m ℓ and, in the positron case, reduces to B(K + → e + ν) = 1.582(7) × 10 −5 [4] for m N → 0 due to helicity suppression. A search for K + → ℓ + N decays in HNL mass range 170-448 MeV/c 2 using a data sample collected with a minimum bias trigger by the NA62 experiment at CERN during the first physics data-taking in 2015 is reported here. The obtained upper limits on |U ℓ4 | 2 complement, and improve on, those obtained in earlier HNL production searches in pion and kaon decays [5][6][7][8][9].

Beam, detector and data sample
The layout of the NA62 beamline and detector [10] is shown schematically in Fig. 1. A secondary positive hadron beam with a central momentum of 75 GeV/c and 1% momentum spread (rms) is derived from primary 400 GeV/c protons extracted from the CERN SPS and interacting with a beryllium target in spills of 3 s effective duration at nominal intensity of 1.1 × 10 12 protons/s. Beam kaons are tagged with a 70 ps time resolution by a differential Cherenkov counter (KTAG) with nitrogen radiator at 1.73 bar pressure contained in a 5 m long vessel. Beam particle momenta are measured by a silicon pixel detector (GTK, under commissioning in 2015 and not used for this analysis). Inelastic interactions of beam particles with the last of the three GTK stations are detected by an array of scintillator hodoscopes (CHANTI). The beam is delivered into a vacuum tank containing a 75 m long fiducial decay volume (FV) starting 2.6 m downstream of the last GTK station. The beam transverse size at the FV entrance is 53 × 24 mm 2 , and the beam divergence in 2015 was 0.22 (0.11) mrad in the horizontal (vertical) plane. The nominal instantaneous particle rate at the FV entrance is 750 MHz, mainly due to π + (70%), protons (23%) and K + (6%). The fraction of kaons decaying in the FV is 13%, leading to 6 MHz nominal K + decay rate. The beam is accompanied by a flux of muons produced by K + and π + decays upstream of the vacuum tank (the beam halo), with 3 MHz nominal rate in the detector acceptance. Central holes in detectors downstream of the FV and a beam pipe traversing most of these detectors allow the undecayed beam particles to continue their path in vacuum. The beam intensity during the 2015 run was typically O(1%) of the nominal value. The momenta of charged K + decay products are measured by a magnetic spectrometer (STRAW) located in the vacuum tank downstream of the FV. The spectrometer consists of four tracking chambers made of straw tubes, and a dipole magnet located between the second and the third chamber providing a horizontal momentum kick of approximately 270 MeV/c. The spectrometer momentum resolution is σ p /p = (0.30 ⊕ 0.005 · p)%, where the momentum p is expressed in GeV/c. A 27X 0 thick quasi-homogeneous liquid krypton (LKr) electromagnetic calorimeter, built for the earlier NA48 experiment [11] and equipped with a new readout system, is used for photon detection. The calorimeter has an active volume of 7 m 3 , and is segmented transversally into 13248 projective ∼ 2 × 2 cm 2 cells. Its energy resolution in the NA62 conditions is σ E /E = (4.8/ √ E ⊕ 11/E ⊕ 0.9)%, where E is expressed in GeV. To achieve hermetic acceptance for photons emitted in K + decays in the FV at angles up to 50 mrad, the LKr calorimeter is supplemented by annular lead glass large-angle veto (LAV) detectors installed in 12 positions along and downstream of the FV, and two lead/scintillator sampling calorimeters (intermediatering calorimeter, IRC, and small-angle calorimeter, SAC) located close to the beam axis.
A ring-imaging Cherenkov detector (RICH) consisting of a 17.5 m long vessel filled with neon at atmospheric pressure is used for identification of charged K + decay products. Its Cherenkov threshold for muons is 9.5 GeV/c, and it provides timing measurement for tracks above threshold to better than 100 ps precision. The LKr calorimeter, a hadronic iron/scintillator sampling calorimeter formed of two modules (MUV1,2) and a scintillator-tile muon detector (MUV3) located behind an 80 cm thick iron wall are also used for particle identification. A plastic scintillator hodoscope (CHOD) built for the NA48 experiment, located in front of the calorimeters, provides a fast trigger with efficiency above 99% and track-timing measurement to 200 ps precision.
The data sample used for this analysis is obtained from 1.2 × 10 4 SPS spills recorded in 5 days of operation in 2015 at beam intensity varying from 0.4% to 1.3% of the nominal value with a minimum-bias trigger scheme. The low-level hardware trigger required a CHOD signal (downscaled typically by a factor of 3) to collect K + decays to muons (which account for 67% of the decay rate), and a CHOD signal in anti-coincidence with a MUV3 signal (not downscaled) to collect decays with no muons in the final state. The high-level software trigger required a kaon signal in the KTAG detector within ±10 ns of the low-level trigger signal. Loose timing conditions are used in this analysis because the accidental rates are small, due to the low beam intensity.

Event selection
Assuming |U ℓ4 | 2 < 10 −4 and considering HNL decays into SM particles [12], the smallest possible average decay length of a HNL produced in the K + → ℓ + N decays in NA62 conditions exceeds 10 km. Under the above assumption, HNL decays in flight in the 156 m long volume from the start of the FV to the last detector (SAC) can be neglected, and the K + → ℓ + N decay is characterized by a single detected track in the final state, similarly to the SM K + → ℓ + ν decay. The principal selection criteria are listed below.
• A single positively charged track reconstructed in the spectrometer with momentum in the range 5-70 GeV/c is required. Additional spectrometer tracks and LKr energy deposition clusters not geometrically compatible with the track are not allowed within ±100 ns of the track time measured by the CHOD. Activity in the large-angle and small-angle photon veto detectors and the CHANTI detector within ±10 ns of the track time is not allowed. Track impact points in the straw chambers, LKr calorimeter, CHOD and MUV1-3 detectors should be within their fiducial geometrical acceptances.
• The kaon decay vertex is reconstructed as the point of closest approach of the track and the beam axis (the latter is monitored with fully reconstructed K + → π + π + π − decays), taking into account the measured stray magnetic field map in the vacuum tank. The reconstructed closest distance of approach (CDA) between the track and beam axis should be less than 25 mm, as determined by the beam transverse size.
• To suppress beam halo background from K + decays upstream of the KTAG and beam π + decays, the presence of a kaon signal in the KTAG is required within ±10 ns of the track time measured by the CHOD.
• Beam halo background from K + → µ + ν decays over the approximately 30 m long path between the KTAG and the last GTK station (with the muon deflected by magnetic fields and scattered in magnet yokes and collimators before reaching the vacuum tank) is suppressed by geometrical conditions established by studies of upstream K + decays. For the K + → e + N selection, the reconstructed vertex position is required to be at least 10 m downstream of the start of the FV. For the K + → µ + N selection, the minimal required distance between the decay vertex and the start of the FV depends on the muon emission angle with respect to the beam axis and lies in the range 10-33 m.
• Positrons and muons are identified by the ratio of energy deposit, E, in the LKr calorimeter to momentum, p, measured by the spectrometer: 0.9 < E/p < 1.15 and E/p < 0.2, respectively. No signals in MUV1-3 detectors within ±20 ns of the track time and geometrically consistent with e + candidate tracks (accounting for detector granularity and multiple scattering) are allowed; MUV1-3 signals are required for µ + candidate tracks. Additionally, an identification algorithm based on the RICH hit pattern is applied for tracks with p < 40 GeV/c.
The squared missing mass is computed as m 2 miss = (P K − P ℓ ) 2 , where P K and P ℓ are the kaon and lepton 4-momenta, respectively. P K is obtained from the beam average 3-momentum (monitored with K + → π + π + π − decays) in the K + mass hypothesis, while P ℓ is evaluated from the reconstructed track 3-momentum in the corresponding ℓ + mass hypothesis.
Simulation of particle interactions with the detector and its response is performed with a Monte Carlo (MC) simulation package based on the Geant4 toolkit [13]. The m 2 miss spectra of the selected events from both data and simulation are displayed in Fig. 2. Signals from the SM leptonic decays K + → ℓ + ν are observed as peaks at m 2 miss = 0 with m 2 miss resolutions of 4.9 (4.7) × 10 −3 GeV 2 /c 4 in the e + (µ + ) case. These resolutions are dominated by the

Measurement principle
The K + → ℓ + N decay rates are measured with respect to the rates of the normalization SM K + → ℓ + ν decays with similar topologies and known branching fractions. The expected numbers of K + → ℓ + N signal events N ℓ S are related to the assumed branching fractions B(K + → ℓ + N ) and acceptances A N ℓ of the K + → ℓ + N selections as Here N ℓ K are the numbers of K + decays in the FV, computed from the numbers N ℓ of selected data events with m 2 miss in the SM signal region: where A ℓ 2 ℓ 1 is the acceptance of the K + → ℓ + 1 ν selection (with m 2 miss in the SM signal region) for the K + → ℓ + 2 ν decay evaluated with MC simulations, and B(K + → ℓ + ν) is the branching fraction of the K + → ℓ + ν decay [4]. The inputs to the computation of N ℓ K are summarized in Table 1: Inputs to the computation of the numbers N ℓ K of kaon decays in the FV: numbers of selected data events in the SM signal region, acceptances evaluated with MC simulations and their statistical errors (notation is defined in the text), and K + → ℓ + ν branching fractions [4].
(1.582 ± 0.007) × 10 −5 0.6356 ± 0.0011 Table 1. The number of K + decays in the µ + case is smaller than that in the e + case due to the downscaling factor of typically 3 applied to the muon trigger chain. The above approach relies on first-order cancellation between signal, normalization and background yields of the effects of residual detector inefficiencies, trigger efficiencies and random veto not fully accounted for by the MC simulation.
The background in the K + → e + ν sample from K + → µ + ν decays, due to both µ + misidentification and decay in flight, is taken into account in the computation of N e K . This background is dominated by µ + mis-identification due to 'catastrophic' bremsstrahlung in the LKr calorimeter at track momenta p > 40 GeV/c, where identification relies on calorimetry only as the RICH does not provide useful information. The probability of a muon having E/p > 0.90 in the LKr calorimeter has been measured in a dedicated study to be P µe ∼ 10 −5 , and found to be reproduced by simulation to 10% relative precision [14]. The background in the K + → µ + ν sample is negligible.
The quoted uncertainty on N e K receives contributions from the statistical error (2.4%), precision on the simulation (evaluated by stability checks versus variation of the selection conditions and considering the precision on P µe simulation, 2.0%), MC statistical precision on the acceptance for the K + → µ + ν background (1.9%) and the external parameter B(K + → e + ν) (0.4%), combined in quadrature to obtain a total relative error of 3.7%. The uncertainty on N µ K receives two contributions of similar size: due to the precision of the simulation (evaluated by variation of the selection conditions) and due to the external input B(K + → µ + ν), combined in quadrature to obtain a total relative error of 1.9%.

Background estimates with MC simulations
The HNL search procedure, presented in Section 5, is based on a data-driven background estimation method, but this is only valid provided there are no peaking background structures in the HNL mass region. Backgrounds to HNL production have been estimated by MC simulations (Fig. 2) to understand qualitatively their contributions and to optimize the event selection. The results of these simulation studies, reported below, justify the adopted procedure.

Backgrounds to K + → e + N
The principal background to K + → e + N decays comes from the K + → µ + ν decay followed by muon decay in flight µ + → e + νν. It is characterized by a broader CDA distribution than the signal, and is suppressed by the CDA and vertex position selection criteria (Section 2). The CDA selection criterion and therefore the background level are determined by the beam transverse size. The background due to K + → µ + ν decays with muon mis-identification (Section 3) is constrained to low m 2 miss values outside the HNL signal region. Beam pion decays π + → e + ν, as well as π + → µ + ν followed by muon decay in flight, contribute to the background via π + mis-identification by the KTAG due to accidental coincidence with a beam kaon not decaying in the FV. The contribution from direct π + mis-identification by the KTAG is negligible. Pion mis-identification probability for the employed KTAG-CHOD timing condition, averaged over the data sample, is computed to be (0.9 ± 0.1 syst )% from the beam K + rate measured via the rate of out-of-time K + signals in the KTAG. This estimate is consistent with the number of observed π + → e + ν decays in the (P π − P e ) 2 spectrum, where P π is the beam pion 4-momentum.
Backgrounds from all other major K + decays with branching fractions above 1%, and all K + decays to positrons and branching fractions above 10 −5 [4] have been considered. The m 2 miss spectra of the estimated background components are displayed in Fig. 2a, showing good agreement with the data spectrum.

Backgrounds to
The largest component of the background to K + → µ + N decays comes from the K + → µ + νγ decay, mainly due to photons emitted at angles greater than 50 mrad with respect to the beam axis and escaping the LAV geometrical acceptance. It is simulated including inner bremsstrahlung and structure-dependent processes as well as their interference [15]; decays with the photon energy in the kaon rest frame E γ below and above 10 MeV are simulated separately to increase the MC statistics in the latter case.
Residual background due to K + decays between the KTAG and the last GTK station, which is suppressed by the cut on the vertex longitudinal position (Section 2), is estimated from a dedicated simulation. Backgrounds from all other major K + decays are also considered: the largest of them is due to the K + → π + π + π − decay.
The m 2 miss spectra of the estimated background components are displayed in Fig. 2b; current agreement with the data spectrum in the HNL signal region is marginal. Observation in the data of a background component (not reproduced with MC simulation) with muons propagating close to the yz plane (which is the bending plane of the GTK dipole magnets) suggests that the data/MC disagreement in the HNL signal region is due to the limited precision on the beamline simulation affecting the estimated background from upstream K + decays. On the other hand, the disagreement at negative m 2 miss is due to the limited precision of the description of the resolution, affected by the simulation of the beam momentum spectrum and divergence.

Search for HNL production
Mass scans are performed in the HNL signal regions with a step size of 1 MeV/c 2 . The event selection employed for each HNL mass hypothesis involves an additional condition: the reconstructed missing mass should be within ±1.5σ ℓ m of the assumed HNL mass, where σ ℓ m is the mass resolution evaluated with MC simulations (Fig. 3a). The above width of the signal mass window leads to near-optimal expected upper limits on B(K + → ℓ + N ) in the absence of signals across the whole HNL signal regions. A loose selection with a relaxed vertex longitudinal position constraint (requiring the vertex to be in the FV) is applied in the K + → e + N case for mass hypotheses of 350 MeV/c 2 and higher, reflecting the fact that the beam halo background does not populate this mass range. Acceptances, A N ℓ , of the selections (including the ±1.5σ ℓ m mass cut) as functions of HNL mass obtained with MC simulations are shown in Fig. 3b.
To certify that the missing mass resolution and therefore the signal acceptance are simulated correctly outside the SM K + → ℓ + ν peaks, the resolution on ∆m 2 3π = (P K − P 3 ) 2 − (P 1 + P 2 ) 2 in fully reconstructed K + → π + π + π − decays, where P i (i = 1, 2, 3) are the pion 4-momenta reconstructed from the spectrometer information and P K is the kaon 4-momentum defined as for the m 2 miss computation (Section 2), has been studied as a function of (P 1 + P 2 ) 2 . Given that ∆m 2 3π = 0 by construction, the resolution on ∆m 2 3π can be measured for both data and MC. Data and MC resolutions have been found to agree within 1%. For the adopted ±1.5σ ℓ m mass window, 1% change on σ ℓ m translates into 0.4% relative change on the signal acceptance. In each HNL mass hypothesis considered, the background is evaluated from sidebands of the data m miss distribution. The number of expected background events N exp within the ±1.5σ ℓ m HNL search window is estimated from a least-squares fit to the data m miss spectrum with a bin width of 5 MeV/c 2 using third order polynomial functions in the 100-460 (200-385) MeV/c 2 range for the e + (µ + ) case. Mass bins overlapping with the ±1.5σ ℓ m wide HNL search window are excluded from the fit to avoid bias caused by possible HNL signals. Statistical uncertainties δN exp on the background estimates N exp are computed by propagation of statistical errors on the fit function parameters: they are typically about 10% in relative terms. Systematic uncertainties on N exp due to the choice of background fit function, estimated by using fourth order polynomials for the fits, are negligible (typically 1%).
In each HNL mass hypothesis, the total number of observed events N obs within the ±1.5σ ℓ m HNL search window, the number of expected background events N exp and its uncertainty δN exp are used to compute confidence intervals for the number of observed K + → ℓ + N decays N ℓ S . The Rolke-López method [16] assuming Poissonian (Gaussian) distributions for the numbers of observed (expected) events is used. The procedure has been tested and found to be unbiased in the presence of artificially injected statistically significant K + → ℓ + N signals. The values of N exp , δN exp and N obs in each HNL mass hypothesis considered are shown in Fig. 4. The maximum value of the local signal significance computed as z = (N obs − N exp )/ N obs + (δN exp ) 2 is 2.2, for the e + case with m N = 283 MeV/c 2 . In the absence of statistically significant HNL production signals, upper limits on N ℓ S are established; the expected and observed limits at 90% CL are shown in Fig. 4. Perfect knowledge of the background (δN exp = 0) would improve these limits typically by 30%.      Figure 6: Upper limits at 90% CL on |U ℓ4 | 2 obtained for each assumed HNL mass compared to the limits established by earlier HNL production searches in π + decays: TRIUMF (1992) [5], PIENU (2017) [6] and K + decays: KEK (1984) [7], E949 (2015) [8], NA62-2007 (2017) [9].
Single event sensitivities (SES) defined as the values of B(K + → ℓ + N ) and the mixing parameter |U ℓ4 | 2 corresponding to the observation of one signal event, are displayed as functions of HNL mass in Fig. 5a. They are O(10 −8 ), and those in the positron case are smaller than those in the muon case due to N e K being larger than N µ K . Upper limits on the branching fraction B(K + → ℓ + N ) in each HNL mass hypothesis are computed from those on N ℓ S using eq. (2); the expected and observed limits at 90% CL are shown in Fig. 5b. Upper limits on the mixing parameter |U ℓ4 | 2 in each HNL mass hypothesis are computed from those on B(K + → ℓ + N ) according to eq. (1). These limits depend on the external inputs B(K + → ℓ + ν) only in the e + case due to the background subtraction in the N e K computation. Systematic uncertainties on the limits are, in relative terms, of the same magnitude as those on N ℓ K (Section 3). The obtained upper limits on |U ℓ4 | 2 at 90% CL together with the limits from previous HNL production searches in π + [5,6] and K + [7][8][9] decays are shown in Fig. 6. The reported result improves the existing limits on both |U e4 | 2 (over the whole mass range considered) and |U µ4 | 2 (above 300 MeV/c 2 ).

Summary
A search for HNL production in K + → ℓ + N decays has been performed with NA62 data recorded in 2015 at ∼ 1% of the nominal beam intensity with a minimum bias trigger. Upper limits on the HNL mixing parameters |U e4 | 2 and |U µ4 | 2 in the ranges 170-448 MeV/c 2 and 250-373 MeV/c 2 , respectively, have been established at the level between 10 −7 and 10 −6 . This improves on the previous limits from HNL production searches over the whole mass range considered for |U e4 | 2 (and extends the mass range in which the limits exist), and above m N = 300 MeV/c 2 for |U µ4 | 2 .