Search for heavy neutral leptons decaying into muon-pion pairs in the MicroBooNE detector

We present upper limits on the production of heavy neutral leptons (HNLs) decaying to μπ pairs using data collected with the MicroBooNE liquid-argon time projection chamber (TPC) operating at Fermilab. This search is the first of its kind performed in a liquid-argon TPC. We use data collected in 2017 and 2018 corresponding to an exposure of 2 . 0 × 10 20 protons on target from the Fermilab Booster Neutrino Beam, which produces mainly muon neutrinos with an average energy of ≈ 800 MeV. HNLs with higher mass are expected to have a longer time of flight to the liquid-argon TPC than Standard Model neutrinos. The data are therefore recorded with a dedicated trigger configured to detect HNL decays that occur after the neutrino spill reaches the detector. We set upper limits at the 90% confidence level on the element j U μ 4 j 2 of the extended PMNS mixing matrix in the range j U μ 4 j 2 < ð 6 . 6 – 0 . 9 Þ × 10 − 7 for Dirac HNLs and j U μ 4 j 2 < ð 4 . 7 – 0 . 7 Þ × 10 − 7 for Majorana HNLs, assuming HNL masses between 260 and 385 MeV and j U e 4 j 2 ¼ j


I. INTRODUCTION
The Standard Model (SM) describes massless neutrinos as left-handed states. The observation of neutrino oscillations [1] has demonstrated, however, that neutrinos must have mass, requiring extensions of the SM, such as the neutrino minimal Standard Model (νMSM) [2,3]. The νMSM predicts additional right-handed neutral leptons that, unlike SM neutrinos, are not charged under the weak interaction and thus manifest themselves only through their mixing with SM neutrinos. According to the νMSM, a right-handed neutral lepton at the keV mass scale can provide a candidate for dark matter, while the other two right-handed leptons are expected to have masses at the GeV scale [2,3]. In general, the masses of these righthanded states and their coupling to SM neutrinos are not predicted by the model, and their allowed values can thus span many orders of magnitude. In this paper, we report results from a search for heavy neutral leptons (HNLs) with masses of Oð100Þ MeV.
The MicroBooNE detector [4] began collecting data from the Booster Neutrino Beam (BNB) [5] in 2015, making it the first fully operational detector of the three liquid-argon time projection chambers comprising the Short-Baseline Neutrino Program (SBN) [6]. The SBN program will address the short-baseline anomalies observed by the MiniBooNE and LSND Collaborations [7,8]. One possible explanation of these anomalies is the existence of light sterile neutrinos with masses of the order of eV, which would lead to short-baseline neutrino flavor oscillations not expected in the SM.
In addition to the studies of the effects of eV-scale sterile neutrinos, the energy range of the BNB allows us to extend the sensitivity of the SBN detectors to the production and decay of HNLs with masses of Oð100Þ MeV. HNLs produced by the BNB would travel along the beam line and could then decay in flight to μπ pairs inside the MicroBooNE detector, located 463 m downstream from the neutrino production target. Due to their mass, some of the HNLs are expected to arrive late compared to the arrival of the BNB spill. To suppress background from SM neutrino interactions, we therefore use data collected with a dedicated HNL trigger. This trigger was commissioned in 2017 and is used to search for late signatures occurring after the arrival of the SM neutrino beam spill. This trigger allows us to perform a search for HNLs in the mass range * microboone_info@fnal.gov 260-385 MeV using data taken in 2017 and 2018 that corresponds to 2.0 × 10 20 protons on target (POT).

II. HEAVY NEUTRAL LEPTONS
We define the HNL in terms of its relevant parameters: its mass m N , and the elements of the extended Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix jU α4 j 2 (α¼e, μ, τ). The flavor eigenstates of the left-handed neutrinos ν α are written as a linear combination of the SM neutrino mass eigenstates ν i (i ¼ 1, 2, 3) and the heavy neutral lepton state, N, in the form HNLs can be produced (and decay) via SM gauge interactions, with a rate suppressed by the relevant jU α4 j 2 element through mixing-mediated interactions with SM gauge bosons. The decay of charged kaons and pions from the BNB can thus produce a flux of HNLs, which then propagate to the MicroBooNE detector, where they are assumed to decay to SM particles. Figure 1 shows diagrams for the production and decay channels. In this paper, we only consider HNL decays to μπ final states, where the HNLs are produced through the process K þ → μ þ N.
HNL states can include both Dirac and Majorana mass terms. Majorana HNLs would decay in equal numbers into μ þ π − and μ − π þ final states. Since the BNB with positive horn polarity used for this search produces predominantly neutrinos and not antineutrinos, Dirac HNLs could only decay through the process N → μ − π þ .
Assuming jU e4 j 2 ¼ jU τ4 j 2 ¼ 0, the HNL production rates and the μπ decay width are both proportional to jU μ4 j 2 only [9], and we therefore place limits exclusively on the jU μ4 j 2 mixing matrix element. The accessible HNL masses are given by the requirement that the decay and production be kinematically allowed, i.e., m K − m μ > m N > m μ þ m π .
The angular distributions of the decay products is given by the angle between the polarization vector of the HNL and the momentum of the charged lepton in the HNL rest frame. The angular distributions differ between the charge combinations (μ þ π − and μ − π þ ) but the combined distribution describing Majorana HNLs is isotropic, and the expected rate is double the rate of Dirac HNL decays [10,11].

III. HNL FLUX IN THE BNB
The BNB impinges protons from the Fermilab Booster synchrotron on a beryllium target. The protons are delivered in a spill with a duration of 1.6 μs and an average repetition rate of up to 5 Hz [5]. The proton kinetic energy of 8 GeV limits the types of mesons produced by p-Be interactions to kaons and pions, generating a muonneutrino beam with average energies of 800 MeV. This restricts the highest mass of HNL that can be produced at the BNB to m K − m μ for jU μ4 j 2 mediated channels.
We calculate the HNL production rate from the BNB using the SM neutrino flux simulation [5]. The decay kinematics of each SM neutrino parent are calculated for an HNL of mass M N . Each event is then weighted by a kinematic factor to account for the effect of M N on the parent decay rate and by a geometric factor describing the probability of the HNL reaching the MicroBooNE detector. The geometric factor enhances the flux since the HNLs with higher mass are boosted into the beam direction.
The kinematic factor suppresses HNL production at the kinematic threshold and takes into account the smaller helicity suppression due to the mass of the HNLs [12][13][14].
Due to the Lorentz transformation into the lab frame, the decay probability becomes inversely proportional to HNL momentum. In contrast, the number of SM neutrino interactions is given by their interaction cross sections on argon, which rises with energy. The HNL flux is thus expected to be enhanced at lower momenta, leading to correspondingly longer travel times to the detector.

IV. MicroBooNE DETECTOR
The MicroBooNE detector [4] is a liquid-argon time projection chamber (TPC) situated at near-ground level at a location 463 m downstream from the target of BNB, receiving a 93.6% pure ν μ beam. The MicroBooNE TPC has an active mass of 85 t of liquid argon, in a volume 2.6 × 2.3 × 10.4 m 3 in the x, y, z coordinates, respectively. The MicroBooNE detector is described by a right-handed coordinate system. The x axis points along the negative drift direction with the origin located at the anode plane, the y axis points vertically upward with the origin at the center of the detector, and the z axis points along the direction of the beam, with the origin at the upstream edge of the detector. The polar angle is defined with respect to the z axis and the azimuthal angle ϕ with respect to the y axis. Neutrinos that cross the detector can interact with the argon nuclei, or, in the case of HNLs, decay to SM particles and produce secondary charged particles that ionize the argon atoms along their trajectories producing ionization electrons and scintillation light. An electric field of 273 V=cm causes the electrons to drift towards the anode plane, requiring 2.3 ms to drift across the width of the detector. The anode planes are positioned perpendicular to the electric field and comprise three planes of sense wires with a spacing of 3 mm between adjacent wires and the same spacing separating the wire planes. Ionization electrons induce a bipolar signal when they pass through the first two planes of wires, oriented at AE60°with respect to the vertical, before being collected on the third plane with vertically oriented wires producing a unipolar signal.
The waveforms measured by the 8192 wires are digitized in a 4.8 ms readout window. The signal processing on the raw TPC waveforms includes noise filtering and deconvolution to convert wire signals into hit information [15]. Subsequently, individual hits corresponding to a localized energy deposit are extracted for each wire. The combination of timing information and energy deposit contained in each waveform is used to create two-dimensional (2D) projective views of the event. Figure 2 shows such a 2D view for a simulated HNL decay. The Pandora [16,17] toolkit is then used to reconstruct three-dimensional (3D) tracks (produced by muons, pions and protons) and showers (produced by electrons and photons) from the 2D views.
A calibration is performed to take into account all the microphysics in the detector, including electron-ion recombination and the space charge effect (SCE) [18]. The SCE is caused by slowly drifting ions produced by cosmic rays that create variations in the electric fields. This variation impacts energy deposits and track trajectories, which appear distorted, particularly near the edges of the TPC. The individual energy deposits are corrected by the timeaveraged calibration factors obtained from data.
An array of 32 8-inch photomultiplier tubes (PMTs) with 16 ns timing resolution collects the scintillation light produced by argon ionization. We use this measurement to determine the time of the neutrino interaction and for triggering. Light flashes are reconstructed with a timing resolution of 100 ns by summing waveforms from the 32 PMTs.

V. TRIGGERS AND DATA SAMPLES
To reduce the amount of recorded data, online software triggers are deployed, processing the waveforms of the light collection system so that activities in coincidence with BNB spills are identified and stored. SM neutrinos arrive at the MicroBooNE detector 1.5 μs after they have been produced, while the time of flight of the HNLs depends on their mass and momentum. As the HNL mass increases, an increasing fraction of HNL events would thus arrive at the detector after the end of the BNB trigger window of 1.9 μs.
HNLs decaying into μπ pairs within the BNB trigger window need to be discriminated from HNL-like backgrounds due to SM charged-current neutrino interactions that include muons, pions, and also protons in the final state. Neutrino-induced charged-current coherent pion production, where no additional activity around the vertex is expected, could present a nearly irreducible background. The only background relevant for HNL decays that occur after the end of the BNB trigger window are crossing cosmic-ray muons, which have a distinctly different topology from the signal. We therefore focus on the HNLs arriving after the end of the BNB trigger window in this analysis.
In June 2017, we introduced an HNL trigger that starts concurrently with the BNB trigger-in coincidence with the beam spill arrival-but extends the trigger window from 1.9 to 2.5 μs. Figure 3 shows the relationship between the time of flight, the initial energy of HNLs, and the trigger windows, illustrating that the HNL trigger retains a larger fraction of HNL decays when the HNL is heavier. The HNL trigger requires the number of photoelectrons (N pe ) recorded by the PMTs to be N pe > 10.5, which is slightly higher than the requirement of N pe > 6.5 for the BNB trigger. This choice optimizes signal efficiency and trigger rate.
The MicroBooNE detector is exposed to a large flux of cosmic-ray muons, traversing at a rate of ≈5.5 kHz, since it is situated just below ground level at a depth of ≈6 m with no significant overburden. Events in coincidence with BNB spill ("on-beam") data contain up to ≈20 cosmic-ray muons within the readout window of 4.8 ms.
For background studies, we collect data sets based on identical trigger settings as for on-beam data, with the exception of the beam coincidence requirement. These "off-beam" data sets were taken with either the BNB or the HNL trigger requirement on the number of photoelectrons. They contain mainly cosmic rays and no SM neutrino interactions.
In summary, we make use of three data samples: (1) On-beam HNL data, taken in coincidence with a BNB neutrino spill. This data set requires an event to fulfil the HNL trigger condition with a veto on the BNB trigger to reject activity produced by neutrino interactions during the BNB spill. (2) On-beam BNB data, taken in coincidence with a BNB neutrino spill and fulfilling the BNB trigger conditions. (3) Off-beam data, taken with identical trigger settings and in a time window of the same length as for the on-beam BNB or HNL data, but at a time when no beam spills are received at the detector.

VI. MONTE CARLO SAMPLES
Simulated data sets are used to evaluate the reconstruction and selection efficiency, to train a boosted decision tree for signal discrimination and to provide a control sample of SM neutrino interactions for validation.
HNL signal samples are simulated for ten different HNL masses. All HNLs are assumed to travel collinearly with the beam axis, i.e., parallel to the longitudinal z axis of the MicroBooNE detector, such that p x ¼ p y ¼ 0. The time distribution of the HNL production is assumed to be uniform within the beam spill, neglecting the ≈19 ns wide bunch structure, which cannot be resolved with the light reconstruction used in this analysis. The calculated arrival time distributions for HNLs and SM neutrinos produced in a BNB spill is shown in Fig. 4. We simulate HNL decays into μ − π þ and μ þ π − final states with isotropic angular distributions. The μ − π þ decays are reweighted to obtain the angular distribution for Dirac HNLs.
Interactions of SM neutrinos in liquid argon, as well as in the material surrounding the detector ("dirt" events), are simulated within the LArSoft [19] framework using the GENIE [20] Monte Carlo program. For HNL signal events, cosmic rays crossing the detector are modeled by overlaying data from zero-bias off-beam data events, whereas the CORSIKA [21] program is used to simulate cosmic rays for SM neutrino interactions. The detector simulation and propagation of secondary particles in liquid argon is simulated with GEANT4 [22,23].

VII. EVENT SELECTION
As this analysis focuses on the HNLs arriving after the end of the BNB trigger window, events containing signal candidates need to pass the HNL trigger condition but not the BNB trigger condition. Since signal candidates are expected to comprise two reconstructed tracks sharing a common vertex, we select vertices reconstructed by Pandora with exactly two associated tracks. At this stage, events can contain more than one such HNL candidate. To reduce background, we apply further selections: (1) Fiducial Volume. The reconstructed vertex associated with the HNL candidate is required to be located in a fiducial volume, defined as a cuboid contained in the active TPC volume. The vertex location must be greater than 12 cm from the border of the active volume along the x axis, 35 cm along the y axis, 25 cm from the upstream edge of the z axis, and 85 cm from the downstream edge. HNL candidates with vertices located within a 1 m wide gap along the z axis in the range 675 < z < 775 cm are also rejected since the wires are not optimally performant in this region [15]. The fiducial volume thus defined corresponds to a mass of 43 t of liquid argon. (2) Vertex Track Distance. We require the starting points of the two tracks associated with the HNL candidate to lie within 5 cm from the location of the reconstructed vertex. (3) Minimum Number of Hits. The tracks associated with an HNL candidate are required to each have more than 30 associated hits in the collection plane, which corresponds to a minimum energy deposit of ≈20 MeV. (4) Flash Distance. We require the distance in the yz plane between the center of the reconstructed light flash and the yz projection of the vertex location to be less than 150 cm. (5) Track Containment. The tracks associated with an HNL candidate have to be fully contained within a volume defined by a distance of 25 cm from the edges of the active volume on the y axis and 10 cm on the x and z axes. This requirement removes tracks crossing the TPC edges that are more severely affected by the SCE. It also rejects cosmic-ray background. (6) Kinematics. A large fraction of cosmic-ray muons misidentified as HNL candidates are caused by "broken" tracks where the reconstruction algorithm has split the muon track, assigning the two sections to a common vertex at the point where the track is broken. Such candidates have a large opening angle Δα between the two tracks and are thus rejected by requiring Δα < 2.8. We also require that the mass as determined from the momenta assigned to the tracks is < 500 MeV. The momenta of the two tracks representing the μπ pair from the HNL decay are determined by the length of the track under the hypothesis that the longer particle track is from the muon and the shorter the pion. The length of the tracks associated with the HNL candidates is of order 10 cm and depends on the HNL mass. The impact of the muonpion assignment on the results is negligible. After applying these selection requirements, (45-50)% of the HNL candidates in the HNL simulation are retained, with a corresponding background selection efficiency of 1.6%. The contribution from dirt events in the background sample are found to be negligible. Table I shows the numbers of HNL candidates for the data samples used for this analysis, corresponding to 2.0 × 10 20 POT, and for the corresponding HNL signal simulation assuming a mass of 370 MeV and a mixing angle of jU μ4 j 2 ¼ 1.4 × 10 −7 . The expected number of background candidates is derived from the off-beam data by normalizing the number of time windows in the off-beam data to the number of beam spills of the on-beam data. The background expectation agrees with the data within the statistical uncertainties.

VIII. SIGNAL EXTRACTION
We train a boosted decision tree (BDT) on a set of kinematic variables to discriminate between signal candidates and background using the XGBoost framework [24]. The following variables are used as input to the BDT: (1) the 3D opening angle Δα between the two tracks associated with the HNL decay; (2) the momentum jp N j of the HNL candidate; (3) the polar angle θ of the HNL candidate; (4) the azimuthal angle ϕ of the HNL candidate; and (5) the invariant mass m N of the μπ pair, where θ and ϕ are defined in the MicroBooNE coordinate system. Figure 5 shows the distributions of these variables for signal HNL candidates with a mass of 370 MeV, and for the off-beam data set with the selection applied. Since we observe no statistically significant difference in the kinematic distributions between the BNB and HNL triggered off-beam data sets, the higher statistics off-beam BNB data set is shown here.
The distributions show good separation between signal and background. Background is constituted mostly of "broken" cosmic-ray tracks. We train separate BDT models for ten different mass hypotheses in the range 260-385 MeV. The BDT score distributions for the signal and for the off-beam background shown in Fig. 6 TABLE I. Number of candidates remaining after the selection requirements for the HNL signal with a mass of 370 MeV and for jU μ4 j 2 ¼ 1.4 × 10 −7 , the expected background rate derived from the off-beam data set, and the on-beam HNL data, corresponding to 2.0 × 10 20 POT. The uncertainty on the expected background rate is given by the statistical uncertainty of the off-beam data. demonstrate good separation between signal and background for different HNL masses. The discrimination between signal and background improves with increased HNL mass.

HNL signal Background Data
The distributions in Fig. 5 are given for HNL candidates, where each event can contain more than one HNL candidate. Since the probability to observe more than one HNL candidate is negligible for the mixing angles considered here, we retain only the candidate with the most signal-like BDT score in each event.

IX. CONTROL SAMPLES
We use a statistically independent control sample to validate the determination of the reconstruction and selection efficiencies, and the BDT performance. The control sample is the on-beam BNB data set, which is expected to contain SM neutrino interactions with similar final-state topologies as HNL decays. In addition to applying the HNL selection to the control sample, we reject HNL candidates where one of the tracks has an associated energy deposit consistent with the specific energy loss expected for a proton. This gives a better representation of the HNL selection, since candidates with protons in the final state are not expected in the signal sample.
The trigger used to record the on-beam BNB sample selects events with any kind of interaction, i.e., the data sample will contain events with cosmic-rays only and no SM neutrino interaction. We therefore subtract the distributions obtained from an off-beam data sample containing cosmic-ray events to obtain distributions statistically representing a sample of mainly SM neutrino interactions.
The numbers of HNL candidates remaining after each preselection step are given in Table II for the control sample and the Monte Carlo (MC) simulation of BNB SM neutrino events. We then apply the same BDTs trained for three HNL masses as used for the signal selection in the previous section to the control samples after the selection to obtain the BDT score distributions in Fig. 7. The distributions for the MC simulation of BNB SM neutrino events are normalized to the same number of POTs as the data. We observe good agreement between data and simulation in terms of shape and normalization. For BDTs trained with higher HNL masses, signal-like events appear at higher BDT scores. In the MC simulation, these events are mostly due to charged-current (CC) neutrino interactions. Kinematics and track length in CC neutrino interactions become more similar to the simulated kinematics for higher mass HNLs. This background is not present in the data samples used for the HNL search.

X. SYSTEMATIC UNCERTAINTIES
Systematic uncertainties on the simulation of the HNL signal originate from the simulation of the HNL flux, the trigger efficiency, and several calibration and reconstruction effects. To quantify the systematic uncertainties, we use a signal-enriched sample with a BDT score > 0.95. The following uncertainties on the signal are considered: (1) HNL flux uncertainties are estimated by simultaneously varying all parameters used for the flux simulation [25]. The relevant parameters are related to the beam-line simulation, which includes variations of the horn current and the skin effect, and to uncertainties on the kaon production cross section. The contributions to the total flux uncertainty from both types of uncertainty are approximately equal. For all HNL masses, the overall uncertainty on the flux is 8%. In addition, a constant 2% uncertainty accounts for the POT counting performed with the beam toroid. (2) Systematic uncertainties on the trigger efficiency originate from the timing resolution of the PMTs, which we use to define the HNL trigger window and vetoes of the BNB trigger window. Uncertainties caused by variations of the light yield are negligible. The trigger uncertainty is estimated to be in the range (5-10)% depending on the mass of the HNL signal. (3) Dynamically induced charge (DIC) refers to the charge induced on the wires beyond the wire closest to the ionization trail [26,27]. This effect impacts the algorithm that determines the deposited charge, especially on induction planes, as well as the pattern recognition. To estimate the impact of DIC effects, we compare samples where DIC is simulated in the region up to the adjacent 20 wires (10 on each side) with a MC simulation that does not model DIC. The resulting uncertainty on the normalization of the signal sample is ≈10%. (4) Distortions caused by the SCE are corrected by timeaveraged calibration factors obtained from data. We estimate the uncertainties originating from the SCE correction by comparing the reconstruction efficiency in the MC samples with and without the simulation of the SCE. The corresponding uncertainty is found to be negligible (< 1%).  only differ in the simulation of the detector response. The response is either obtained from data after detector calibration or from simulation. The resulting uncertainties due to the detector response are estimated to be small (< 1%). Uncertainties from nuclear interaction modeling are negligible, since the main background source is cosmic-ray muons. For Dirac HNLs, we consider only μ − π þ decays and not the sum of the charge combinations. The difference in efficiency between the charge combinations is ≈ð2-3Þ%, which is small compared to the total systematic uncertainty. This difference is thus neglected.
The contribution of the systematic uncertainties on the signal efficiency in the signal-enriched sample are summarized in Table III for a HNL mass value of 325 MeV. The systematic uncertainties grow linearly with HNL mass, from 10% to 18% in the mass range 260-385 MeV. Uncertainties on the background estimation, which is derived from data, is dominated by the statistical fluctuations of the off-beam HNL data sample.

XI. RESULTS AND DISCUSSION
The BDT score distributions for signal, background expectation (off-beam HNL data) and data are shown in Fig. 8. Signal and background are well separated, and no data excess is observed in the signal region with high BDT scores. We therefore proceed to set limits on the HNL production rate as a function of mass.
The limits are determined using the modified frequentist CL s method [28][29][30]. We calculate a log-likelihood ratio (LLR) test statistic using Poisson probabilities for estimated background events, signal yields, and the observed number of events for different HNL mass hypotheses. The confidence levels are derived by integrating the LLR distribution in pseudoexperiments using both the signalplus-background (CL sþb ) and the background-only hypotheses (CL b ). The excluded signal rate is defined by the signal strength for which the confidence level for signal, CL s ¼ CL sþb =CL b , equals 0.1.
Systematic uncertainties on both background and signal are taken into account using Gaussian priors. The total systematic uncertainty of approximately 15% on the signal is found to have a negligible impact on the sensitivity, since the uncertainty is dominated by the statistics of the background sample.
To calculate limits, we split the distribution into a signalenriched and signal-depleted region with BDT scores >0.95 and 0.5<BDTscore<0.95, respectively. In Table IV, we compare the expected number of background events to the data for the different HNL mass hypotheses, as well as for the signal-enriched and signal-depleted samples. The expected number of HNL signal events is calculated assuming jU μ4 j 2 ¼ 1.4 × 10 −7 .
The observed upper limits at the 90% confidence level as a function of mass are presented in Fig. 9 and Table V for Majorana HNLs, together with the median expected limit and the 1 and 2 standard deviation bands on the median expected limit. The bands are asymmetric since Poisson statistics are used where the expected number of events is small. The observed and expected limit agree within 1 standard deviation over the entire mass range. The decay  Table V by a factor of ffiffi ffi 2 p . The results are of similar or better sensitivity as those obtained by the NA62 [31] and NuTeV [32] Collaborations for the same mixing parameter in the overlapping HNL mass range. The E949 Collaboration [33] sets limits in a lower mass range between 175 and 300 MeV by measuring K þ meson decays at rest. The PS191 [34,35] and T2K [36] Collaborations place more stringent limits in the HNL mass range between 260 and 360 MeV, but their location at an off-axis angle of 2 degrees restricts the sensitivity to slightly lower masses. The MicroBooNE detector is located on axis, allowing it to set the most constraining limits in the mass range up to 385 MeV.

XII. CONCLUSION
We present the first search for HNLs in a liquid-argon TPC using data recorded with the MicroBooNE detector in 2017-2018 with a novel trigger recording events that arrive at the MicroBooNE detector after the BNB beam spill. The data correspond to 2.0 × 10 20 POT. We assume that the HNLs are produced in kaon decays and decay exclusively into a muon pion final state. We obtain constraints on the element jU μ4 j 2 of the extended PMNS mixing matrix of jU μ4 j 2 <ð4.7-0.7Þ×10 −7 for Majorana HNLs and jU μ4 j 2 < ð6.6-0.9Þ × 10 −7 for Dirac HNLs with masses between 260 and 385 MeV and assuming jU e4 j 2 ¼ jU τ4 j 2 ¼ 0.  9. Limits on jU μ4 j 2 at the 90% confidence level as function of mass for a Majorana and Dirac HNL decaying into μπ pairs. The observed limit is compared to the median expected limit with the 1 and 2 standard deviation (σ) bands.