We present the first search for heavy neutral leptons (HNL) decaying into νe + e -or νπ 0 final states in a liquid-argon time projection chamber using data collected with the MicroBooNE detector. The data were recorded synchronously with the NuMI neutrino beam from Fermilab's Main Injector corresponding to a total exposure of 7.01 × 10 20 protons on target. We set upper limits at the 90% confidence level on the mixing parameter |Uµ4| 2 in the mass ranges 10 ≤ mHNL ≤ 150 MeV for the νe + e -channel and 150 ≤ mHNL ≤ 245 MeV for the νπ 0 channel, assuming |Ue4| 2 = |Uτ4| 2 = 0. These limits represent the most stringent constraints in the mass range 35 < mHNL < 175 MeV and the first constraints from a direct search for νπ 0 decays.

Heavy neutral leptons (HNL) appear in minimal extensions of the standard model (SM) that can explain the origin of neutrino masses, the generation of the baryon asymmetry through leptogenesis, and the nature of dark matter [1]. They are introduced through an extension of the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix by adding heavy mass eigenstates that mix very weakly with the three active neutrino states. For a single HNL state, the extended PMNS matrix has the dimension 4 × 4, which leads to four new parameters: the HNL mass m HNL and three mixing parameters, |U α4 | 2 with α = e, µ, or τ . The HNL production and decay rates are suppressed by the elements |U α4 | 2 through mixingmediated interactions with SM gauge bosons. A vibrant experimental program is dedicated to searching for HNLs and other feebly-interacting particles [2].

Here, we use data recorded with the MicroBooNE detector to perform a search for HNLs decaying to νe + e - or νπ 0 final states. The MicroBooNE detector [3] is one of the three liquid-argon time projection chambers (LArTPC) comprising the Fermilab short-baseline neutrino program [4]. The liquid-argon technology provides a powerful tool to search for signatures of physics beyond the SM as it allows us to fully reconstruct decays through its precision imaging capability. Since the two final states, νe + e -and νπ 0 (π 0 → γγ), are topologically very similar, leading to two electromagnetic showers in the LArTPC, the search is performed within a single analysis framework using boosted decision trees (BDTs). This analysis strategy is based on our previous searches for HNL decays to µπ final states [5] and decays of Higgs portal scalars into e + e -pairs [6].

The MicroBooNE detector recorded data between 2015 and 2021. It was simultaneously exposed on-axis to the booster neutrino beam (BNB) [7] and off-axis to the neutrino beam from the main injector (NuMI) [8]. Only NuMI data are used for this search, since the higher average beam energy compared to the BNB leads to a higher kaon rate and therefore potentially more HNL production. We assume that the HNLs are produced in the absorber, made from aluminium, steel, and concrete, which is located ≈ 725 m downstream from the NuMI beam's graphite target and ≈ 104 m from the MicroBooNE detector [8]. The absorber is located downstream of the MicroBooNE detector at the end of the NuMI decay pipe.

HNLs produced in the absorber would approach the detector in almost the opposite direction to the neutrinos that originate from the NuMI beam target, which significantly improves background rejection [5]. Approximately 13% of the beam protons reach the absorber and produce K + mesons that can decay at rest into HNLs through the process K + → µ + N , while most of the K -mesons are absorbed. If the HNL lifetime is sufficiently long, the HNLs could reach the MicroBooNE detector and decay into SM particles within the argon. The sensitivity of MicroBooNE to this production mechanism has previously been studied in Ref. [10].

The kinematic distributions of the final state particles in the HNL decay depend on m HNL , the kinetic energy of the HNL, and whether the HNL is assumed to be a Dirac or Majorana particle. We encode the production and decay properties of the HNL using the equations in Ref. [9] in a simulation code developed for MicroBooNE's previous HNL search [5]. We also validated the simulation with a recent implementation of HNL kinematics in the GENIE generator [11]. The branching ratios for m HNL < 300 MeV are shown in Fig. 1. We assume

already severely constrained [2] and |U τ 4 | 2 is not kinematically accessible. Neglecting the "invisible" decay N → 3ν, the decays into νe + e -final states dominate for m HNL values below the mass of the π 0 meson, and the νπ 0 final states dominate above. We generate samples for different m HNL , and in the νe + e -and νπ 0 final states, to cover the full range of accessible model parameters.

We use NuMI data corresponding to 7.01 × 10 20 protons on target (POT), which were taken in two operating modes, forward horn current (FHC) with 2.00 × 10 20 POT (Run 1) and reverse horn current (RHC) with 5.01 × 10 20 POT (Run 3). The two data sets are analyzed separately to account for differences in neutrino flux and detector configuration. We assume equal rates of K + production for the two horn polarities [12].

We select a "beam-on" data sample to search for an HNL signal where the event triggers coincide with the NuMI beam. Such beam-on events are frequently triggered by a cosmic ray and not a neutrino interaction. This type of event is modeled by selecting a "beam-off" sample collected under identical trigger conditions but when no neutrino beam is present. The "beam-off" sample is normalized to the number of triggers recorded in the beam-on data. Neutrino-induced background from the NuMI beam is modeled using a Monte Carlo (MC) simulation [13], with cosmic rays and noise from data overlaid on the simulation. The "in-cryostat ν" sample contains interactions of neutrinos with the argon inside the cryostat, and the "out-of-cryostat ν" sample describes interactions with the detector structure and surrounding material. Both samples are normalized to the numbers of POT of the data sample. An additional data-driven scaling factor is applied to the out-of-cryostat ν sample.

We reconstruct neutrino interactions and cosmic rays within the argon with a chain of pattern-recognition algorithms, implemented using the Pandora Software Development Kit (SDK) [14,15]. Hits are formed from the waveforms read out by three anode wire planes -two induction planes and one charge collection plane. We then group hits into slices to isolate neutrino interactions and cosmic rays. Slices are reconstructed under both hypotheses, and a support vector machine then calculates a "topological score" to classify slices as either a neutrino interaction or cosmic ray. We select events with exactly one neutrino slice to examine them for candidate HNLs.

"Objects" are then reconstructed either as a track, as expected for a minimum ionizing particle, or a shower, consistent with being an electron or photon. The distinction between tracks and showers is performed using a "track score" that mainly relies on the profile of the charge deposition, the range, and topological information.

The start and end points of all objects associated with the slice must lie within the TPC's fiducial volume [16], and the fraction of reconstructed hits in the slice contained within the fiducial volume has to be > 0.9. The energy of all the objects in the slice, E sl , is reconstructed from the charge read-out on the TPC's charge collection plane. We require E sl < 500 MeV as the energy deposited from the decays of HNLs with a mass m HNL < 245 MeV is expected to be lower than for most beam or cosmicray events. We require E sl < 500 MeV as the decays of HNLs with a mass m HNL < 245 MeV are expected to deposit less energy than most neutrino or cosmic-ray interactions.

Light flashes are reconstructed from the waveforms of an array of 32 photomultiplier tubes. We require that the time of the largest flash in a 23 µs window surrounding the NuMI beam trigger coincides with the NuMI beam spill of ≈ 10µs. A cosmic ray tagger (CRT) surrounding the cryostat was installed about midway through Micro-BooNE operations [17]. If there is a hit recorded by the CRT within 1 µs of the flash for the Run 3 (RHC) sample, the event is identified as a cosmic ray and is rejected [18]. The CRT is not used for the Run 1 (FHC) data set as it was not yet operational at that time. To further reduce cosmic-ray background, we require the "flash match score" to be < 15. This is calculated as a χ 2 value by comparing the light signals in the PMTs to the expected PMT signals assuming the recorded charge is due to a neutrino interaction.

Table I shows the effect of the preselection requirements for the background samples. The signal efficiency after preselection is ≈ 35% for the different m HNL and final states, while we retain ≈ 4% of the in-cryostat neutrino interactions. The contribution of the beam-off and out-of-cryostat ν events to the background sample is significantly smaller for Run 3 compared to Run 1, since the CRT improves the rejection of these classes of events. (a) (b) (c)

TABLE I. Numbers of events that remain after preselection normalized to the POT for the two data samples. The percentages are the contributions of each sample to the sum of the background predictions. Sample Run 1 (FHC) Run 3 (RHC) POT 2.00 × 10 20 5.01 × 10 20 Beam-off 3548 (46%) 3597 (33%) In-cryostat ν 3607 (47%) 6805 (63%) Out-of-cryostat ν 567 (7%) 464 (4%) Sum of predictions 7722 10866 Beam-on (data) 7598 11282 Data over prediction 0.98 1.04

The numbers of data events after the preselection agree well, within (2-4)%, with the sum of the predictions from the three main background sources. At this stage of the selection, we use XGBoost [19] to train BDTs that optimize the discrimination between signal and background in this selected sample. A separate training is performed for each mass point, final state, and data set (FHC/RHC) using subsets of the signal sample and the three background samples. The training samples are excluded from the subsequent analysis. We reduce the number of input variables to 20 from a potential set of several hundred variables by training BDT models on the full set of variables and then identifying the variables with the largest impact that are common to all tested HNL model parameters.

As variables defined for the event (slice), we use the multiplicity of objects, the track multiplicity, the total energy measured from all tracks, the total energy mea-sured from all showers, the total energy, the energy of the highest energy track, the topological and flash match scores, the energy deposited in the first 4 cm of the highest energy shower, and the shower angle θ yz , which is the average direction of all showers calculated with respect to the y axis projected onto the yz plane.

We also use the number of hits on each of the three wire planes associated with the highest-energy object and the object's track score, where the highest energy object is determined by the associated number of hits on the wires of the collection plane.

The final set of variables use angular information for the highest energy object: the polar angles, the azimuthal angles, and the z-momentum fractions, calculated as the ratio of the momentum component in the z coordinate over the total object momentum. These angular variables are calculated twice for each object by fitting it as track and as a shower.

If there is no track or shower, some variables are treated as missing in the BDT by using a placeholder. The distributions in Fig. 2 show the shower angle θ yz , the track-fit z momentum fraction, and the total shower energy for data and the background prediction in Run 3. These variables were found, during BDT training procedures, to be amongst the most sensitive to an HNL signal. Momenta of particles produced in neutrino interactions predominantly point in the +z direction, whereas signal is more clustered around -z. The shower angle for signal has two peaks, depending on whether the start and end point of the shower are correctly identified. Since the BDT uses the information of all variables, such incorrectly reconstructed events can still be identified as signal.

The background contributing to the BDT score distri-bution shown in Fig. 3 is expected to be dominated by incryostat ν interactions. The BDT identifies and rejects most charged-current ν µ interactions. For BDT scores > 3, about 40% of the simulated in-cryostat ν events are neutral-current interactions producing π 0 mesons, as this topology resembles the N → νe + e -and N → νπ 0 decays. To determine the sensitivity to a possible HNL signal, we evaluate systematic uncertainties that could modify the BDT score distributions for signal and background [20]. For the in-cryostat ν background, we consider the impact of the flux simulation, the neutrinoargon cross-section modeling, hadron interactions with argon, and detector modeling. The beam-off sample is taken from data and therefore has no associated systematic uncertainties other than the statistical fluctuations in the sample. The impact of the normalization uncertainty on the out-of-cryostat ν sample is negligible, as the contribution to the final sample is small [5].

The dominant uncertainty on the background in the signal region at high BDT scores is due to the statistical uncertainty of the samples, since most of the background has been rejected. We therefore extrapolate systematic uncertainties from higher-statistics regions of the BDT score distribution to the signal region. The quadrature sum of the background detector modelling uncertainty is taken to be 30%.

The dominant contribution to the systematic uncertainty on the signal sample arises from the rate of kaon production at rest in the NuMI absorber. It is taken to be ±30% based on the evaluation by the MiniBooNE collaboration [12]. The sum of the detector-related systematic uncertainties is (10 -20)%. The systematic uncertainties are separately evaluated for all signal parameters used in the BDT training, with consistent results. Due to the higher number of POTs and the better cosmic-ray rejection of the CRT, the signal sensitivity is dominated by the Run 3 data set.

The BDT score distributions are used to derive limits on |U µ4 | 2 for the different model parameters. We use the pyhf algorithm [21], which is an implementation of a statistical model to estimate confidence intervals for multi-bin histograms, based on the asymptotic formulas of Ref. [22]. The formalism allows for the treatment of systematic uncertainties through the use of profile likelihood ratios. The results are validated with the modified frequentist CL s calculation of the COLLIE program [23]. The pyhf code scans over a range of scaling parameters for the signal normalization and returns an interpolated value of the scaling parameter that corresponds to the 90% confidence level (CL s = 0.1). The BDT distributions for each run period (Run 1 and Run 3) enter the limit setting as separate channels before their likelihoods are combined. The statistical uncertainties on signal and background are uncorrelated, whereas the systematic uncertainties on the flux and out-of-cryostat ν normalization are taken as fully correlated between the run periods. We studied the impact of the other systematic uncertainties with different assumptions about their correlations and determined that this choice has only a small impact on the result. Using BDT models trained with neighboring mass points shows no significant deterioration of sensitivity.

TABLE II. The 90% CL observed and median expected limits on |Uµ4| 2 as a function of mHNL assuming a Majorana state. mHNL Limit |Uµ4| 2 (MeV) Observed Median Standard deviation νe + e - 10 3.32 × 10 -3 3.59 × 10 -3 2.92-4.52 × 10 -3 20 3.82 × 10 -4 3.94 × 10 -4 3.18-5.02 × 10 -4 50 2.96 × 10 -5 2.86 × 10 -5 2.36-3.55 × 10 -5 100 2.94 × 10 -6 3.16 × 10 -6 2.61-3.92 × 10 -6 150 5.99 × 10 -7 7.71 × 10 -7 6.30-9.71 × 10 -7 νπ 0 150 2.15 × 10 -7 2.25 × 10 -7 1.86-2.79 × 10 -7 180 7.01 × 10 -8 6.87 × 10 -8 5.64-8.56 × 10 -8 200 3.95 × 10 -8 4.43 × 10 -8 3.65-5.51 × 10 -8 220 3.97 × 10 -8 3.62 × 10 -8 3.00-4.46 × 10 -8 240 2.67 × 10 -8 2.71 × 10 -8 2.23-3.36 × 10 -8 245 2.26 × 10 -8 2.29 × 10 -8 1.87-2.85 × 10 -8

The observed limits for the model parameters tested, given in Table II, are all within 2 standard deviations of the median expected limit. A linear interpolation is performed between the tested m HNL hypotheses, which slightly underestimates the sensitivity in the interpolation regions.

We also consider Dirac HNL states. The Dirac HNL decay rate is a factor of 2 smaller than for Majorana

0 50 100 150 200 245 HNL mass [MeV] 10 9 10 8 10 7 10 6 10 5 10 4 10 3 10 2 |U 4 | 2 limit at 90% C.L. SIN PIENU E949 KEK E89 KEK E89/E104 PS191 e NA62 MicroBooNE FIG. 4. Limits on |Uµ4| 2 at the 90% CL assuming |Ue4| 2 = |Uτ4|

2 = 0 as a function of mass for a Majorana HNL in the ranges 10 < mHNL < 150 MeV (N → νe + e -) and 150 < mHNL < 245 MeV (N → νπ 0 ). We use a linear interpolation between the tested mHNL hypotheses. The discontinuity at mHNL = 150 MeV is due to the change in decay channel from νe + e -to νπ 0 . The constraints are compared to the published results of the SIN [24], PIENU [25], KEK-E89 [26], BNL-E949 [27], NA62 [28], and PS191 [29] collaborations. The unpublished limit using the KEK-E89/E104 data [30] is shown as a dashed line.

states at the same value of |U µ4 | 2 . Since the effect of differing decay kinematics on the sensitivity is found to be negligible, the limits for Dirac states can be obtained by scaling the results in Table II by a factor of √ 2.

We compare our results for Majorana HNLs with existing constraints on |U µ4 | 2 in Fig. 4 for |U e4 | 2 = |U τ 4 | 2 = 0. Decays of stopped pions in the process π + → µ + ν are sensitive to the range m HNL < m π -m µ . Such searches have been performed at the Swiss National Institute (SIN) in the mass range 1 < m HNL < 16 MeV [24] and by the PIENU Collaboration for 15.7 < m HNL < 33.8 MeV [25]. The muon spectrum measured in stopped K + → µ + ν decays (K 2µ ) has been used to set limits in the mass range 70 < m HNL < 300 MeV with the E89 experiment at KEK [26], in the range 175 < m HNL < 300 MeV with the E949 experiment at the Brookhaven National Laboratory (BNL) [27], and 200 < m HNL < 384 MeV by the NA62 experiment at CERN [28]. An update of the KEK-E89 result was reported in proceedings [30], extending the sensitivity to m HNL > 40 MeV. The PS191 experiment [29] was specifically designed to search for massive decaying neutrinos in the CERN-PS proton beam. Its search for N → νµe decays constrains |U µ4 | 2 in the range m HNL > m µ . Critical discussions of the PS191 results can be found in Refs. [31][32][33][34]. The T2K Collaboration has published combined limits on |U µ4 | 2 , |U e4 | 2 , and |U τ 4 | 2 for m HNL > 150 MeV in a Bayesian approach [35].

In this letter, we use NuMI beam data recorded with the MicroBooNE detector to derive the most stringent constraints on |U µ4 | 2 for the mass range 34 < m HNL < 175 MeV. It is also the first search for HNLs in νπ 0 and νe + e -final states using a LArTPC and the first ever result reported for HNL decays into νπ 0 . When combined with the MicroBooNE limits on HNL decays into µπ pairs [5], our results now cover the full mass range 10 < m HNL < 385 MeV that is kinematically accessible from kaons produced by the NuMI beam. struction, and operation of the MicroBooNE detector as well as the contributions of past collaborators to the development of MicroBooNE analyses, without whom this work would not have been possible. For the purpose of open access, the authors have applied a Creative Commons Attribution (CC BY) public copyright license to any Author Accepted Manuscript version arising from this submission.