Search for long-lived heavy neutral leptons and Higgs portal scalars decaying in the MicroBooNE detector

We present a search for long-lived Higgs portal scalars (HPS) and heavy neutral leptons (HNL) decaying in the MicroBooNE liquid-argon time projection chamber. The measurement is performed using data collected synchronously with the NuMI neutrino beam from Fermilab's Main Injector with a total exposure corresponding to $7.01 \times 10^{20}$ protons on target. We set upper limits at the $90\%$ confidence level on the mixing parameter $\lvert U_{\mu 4}\rvert^2$ ranging from $\lvert U_{\mu 4}\rvert^2<12.9\times 10^{-8}$ for Majorana HNLs with a mass of $m_{\rm HNL}=246$ MeV to $\lvert U_{\mu 4}\rvert^2<0.92 \times 10^{-8}$ for $m_{\rm HNL}=385$ MeV, assuming $\lvert U_{e 4}\rvert^2 = \lvert U_{\tau 4}\rvert^2 = 0$ and HNL decays into $\mu^\pm\pi^\mp$ pairs. These limits on $\lvert U_{\mu 4}\rvert^2$ represent an order of magnitude improvement in sensitivity compared to the previous MicroBooNE result. We also constrain the scalar-Higgs mixing angle $\theta$ by searching for HPS decays into $\mu^+\mu^-$ final states, excluding a contour in the parameter space with lower bounds of $\theta^2<31.3 \times 10^{-9}$ for $m_{\rm HPS}=212$ GeV and $\theta^2<1.09 \times 10^{-9}$ for $m_{\rm HPS}=275$ GeV. These are the first constraints on the scalar-Higgs mixing angle $\theta$ from a dedicated experimental search in this mass range.


I. INTRODUCTION
The MicroBooNE detector [1] began collecting data in 2015, making it the first fully operational detector of the three liquid-argon time projection chambers comprising the Short-Baseline Neutrino (SBN) program [2] at Fermilab. The MicroBooNE detector was exposed to both the booster neutrino beam (BNB) [3] and the neutrino beam from the main injector (NuMI) [4].
We can use the NuMI beam to study beyond-the-Standard Model (BSM) phenomena such as the production and decay of heavy neutral leptons (HNL) or Higgs portal scalars (HPS), jointly referred to as long-lived particles (LLPs). In addition, we have also used it to measure electron neutrino cross sections on argon [5,6].
In this paper, we present searches for both types of LLPs originating from kaons decaying at rest in the NuMI hadron absorber, which is located downstream of the MicroBooNE detector at the end of the NuMI decay pipe. The LLP would travel ≈ 104 m to the MicroBooNE detector where it can be detected through its decay. LLPs produced in the absorber would approach the detector in almost the opposite direction to the vast majority of neutrinos, which originate from near the NuMI beam target (Fig. 1).
Event displays of simulated HNL and HPS decays in the MicroBooNE detector are shown in Fig.2. The signal topology is characterized by exactly two tracks emerging from a common vertex. Since we search for LLPs produced by two-body decays of kaons at rest, these LLPs have a fixed energy for a given mass. These two properties, the direction and energy of the signal LLPs, help to * microboone info@fnal.gov discriminate them from neutrino and cosmic-ray induced background processes. In addition, the kinematics and topologies of HPS decays to µµ pairs and HNL decays to µπ pairs are similar, which allows us to develop a single signal analysis strategy [7]. The MicroBooNE collaboration has published upper limits on the production of HNLs decaying to µπ pairs for an exposure of 2.0 × 10 20 protons on target (POT) from the BNB, using a dedicated trigger configured to detect HNL decays that occur after the neutrino spill reaches the detector. That search yielded upper limits at the 90% confidence level (CL) on the element |U µ4 | 2 of the extended Pontecorvo-Maki-Nakagawa-Sakata (PMNS) mixing matrix |U µ4 | 2 for Dirac and Majorana HNLs in the HNL mass range 260 ≤ m HNL ≤ 385 MeV and assuming |U e4 | 2 = |U τ 4 | 2 = 0 [8]. Several other collaborations have also set limits on |U µ4 | 2 [9][10][11][12][13][14][15][16]. Also using a liquid-argon detector exposed to the NuMI beam, the ArgoNeuT collaboration has published a search for HNL decays into νµ + µ − final states [17]. They derive limits on the mixing matrix element |U τ 4 | 2 .
The MicroBooNE collaboration has also published a search for HPS decaying to e + e − pairs assuming the HPS originate from kaons decaying at rest after having been produced at the NuMI absorber [18]. A data set corresponding to 1.93 × 10 20 POT is used to set limits on θ 2 in the range 10 −6 -10 −7 at the 95% CL for the mass range directly below the range considered in this paper (0 < m HPS < 211 MeV). Other direct experimental searches for HPS have been published in Refs. [19][20][21][22][23], and reinterpretations of experimental data as HPS limits have been derived in Refs. [24][25][26][27].

II. HEAVY NEUTRAL LEPTONS
The HNL is introduced through an extension of the PMNS matrix by adding a single heavy mass eigenstate that mixes very weakly with the three active neutrino states. This minimal extension adds four parameters to the model comprising the HNL mass m HNL and the ele-ments of the extended PMNS matrix, |U α4 | 2 (α = e, µ, τ ). The flavor eigenstates are with a heavy neutral lepton state N . The HNL production and decay rates are suppressed by the relevant |U α4 | 2 element through mixing-mediated interactions with SM gauge bosons.
Production of a Majorana HNL state N via mixing in a K + meson decay and its subsequent decay into a µ ∓ π ± pair. A Dirac HNL will only decay into µ − π + pairs. HNLs would be produced in the decays of charged kaons and pions originating from the proton interactions on the targets of the BNB or NuMI neutrino beams. If the HNL lifetime is sufficiently long to allow the HNL to reach the MicroBooNE detector, they can decay into Standard Model (SM) particles within the argon volume.
We consider the production channel K + → µ + N with the decay N → µ ∓ π ± as shown in Fig. 3. The HNL production rate and decay width into µπ are each proportional to |U µ4 | 2 , and the total rate therefore to |U µ4 | 4 , assuming |U e4 | 2 = |U τ 4 | 2 = 0 [28]. We thus place limits exclusively on the |U µ4 | 2 mixing matrix element. The accessible HNL masses are constrained by the condition m K − m µ > m HNL > m µ + m π .
FIG. 4. Production of an HPS boson S via mixing in a kaon meson decay and its subsequent decay. Only K + mesons contribute to this analysis, although this decay mode can also occur for neutral or K − mesons.
HNL states can include Dirac and Majorana mass terms, where Majorana HNLs would decay in equal numbers into µ + π − and µ − π + final states. Dirac HNLs from K + decays could only decay to the charge combination µ − π + to conserve lepton number.

III. HIGGS PORTAL SCALARS
The Higgs portal model [29] is an extension to the SM, where an electrically-neutral singlet scalar boson mixes with the Higgs boson with a mixing angle θ. Through this mixing, this HPS boson acquires a coupling to SM fermions via their Yukawa couplings to the Higgs boson, which is proportional to sin θ. The phenomenology of the Higgs portal model at the SBN program, including the equations describing production and decay of the scalar boson, are discussed in Ref. [30]. At the absorber the dominant production channel will be through the twobody decay K + → π + S (where the HPS is denoted by S). The dominant decay mechanism is a penguin diagram with a top quark contributing in the loop (Fig. 4).
The partial width for decays to charged leptons with mass m is proportional to m 2 [30]. If we assume that there are no new dark sector particles with masses < m HPS /2, the branching fraction into µ + µ − pair is therefore close to 100% for scalar masses in the range m µ + µ − < m HPS < m π 0 π 0 . The decays into π 0 π 0 pairs become accessible for m HPS > 269 MeV and the decays into π + π − pairs at m HPS > 279.1 MeV.
We do not consider decays to neutral pion pairs in this search, since the detector signature differs significantly from muon and charged pions. The π + π − decay signatures appear very similar to µ + µ − decays. However, the analysis is not sensitive to HPS decays in this decay channel, as the HPS will decay before it reaches the detector due its short lifetime. We set limits as a function of the mixing angle θ as both HPS production and decay rate are proportional to θ 2 .

IV. MICROBOONE DETECTOR
The MicroBooNE detector [1] is a liquid-argon time projection chamber (LArTPC) situated at near-ground level at Fermilab. It is located at an angle of 8 • relative to the NuMI beamline [4] and at a distance of 680 m from the target. The MicroBooNE LArTPC 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, and z coordinates of the MicroBooNE coordinate system 1 .
Charged particles produced in neutrino interactions with argon or in decays of LLPs will 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 plane is oriented perpendicular to the electric field and comprises 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 ±60 • 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. This is longer than the 2.3 ms drift time to allow reconstructing out-of-time cosmic rays. The signal processing on the raw Time Projection Chamber (TPC) waveforms includes noise filtering and deconvolution to convert wire signals into hit information [31,32]. 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 2D projective views of the event.
An array of 32 8-inch photomultiplier tubes (PMTs) collects the scintillation light produced by argon ionization. Light flashes are reconstructed from the waveforms of the 32 PMTs. To record an event, the NuMI online trigger requires a scintillation light signal in time with the accelerator beam spill window above a configurable threshold for the number of photoelectrons, which was 9.5 for Run 1 and then lowered to 5.75 during Run 3.
A cosmic ray tagger (CRT) surrounding the cyrostat was installed about midway through MicroBooNE operations [33]. It comprises four panels made up of interleaved plastic scintillator strips placed above, below and on the two sides parallel to the BNB beam direction. The CRT provides both fast timing and positional information of cosmic rays entering the TPC.

V. NUMI BEAMLINE
Protons with an energy of 120 GeV from the Main Injector hit a graphite target, producing particles in the NuMI beamline. The position of the MicroBooNE detector relative to the components of the NuMI beamline is shown in Fig. 5. A system of electromagnetic horns focuses the charged particles either towards or away from the beam axis, depending on the horn polarity. In forward horn current (FHC) mode, a positive (+200 kA) current is applied to the horns, which focuses positively charged particles in the beam direction. In reverse horn current (RHC) mode, a negative current (−200 kA) is applied to focus negatively charged particles. In this paper, we use NuMI data collected in both modes.
The focused particles then travel down a 675 m long decay pipe filled with helium where they decay to neutrinos or anti-neutrinos. The proton beam structure determines the intensity and timing structure of the neutrino beam. The neutrino beam has six batches which together form a spill. Each spill is 9.6 µs long.
Immediately downstream of the decay pipe, at a distance of ≈ 104 m from the MicroBooNE detector, is an absorber (5.5 m wide, 5.6 m tall, and 8.5 m deep) made of an aluminium core surrounded by steel and then concrete, designed to absorb the remaining hadrons. In the MicroBooNE coordinate system, the direction from the absorber corresponds to θ = 2.20 and φ = 1.15 (in radians). Approximately 13% of the beam protons pass through the target without interacting, travelling along the decay pipe before colliding with the absorber at a distance of 725 m downstream from the target [4]. These collisions produce a large number of K + mesons which then decay at rest, while most of the K − mesons are absorbed. The LLPs studied in this paper would be produced in this absorber from K + decays. We assume equal rates of K + production for the two horn polarities [34].

VI. FLUX GENERATION
The flux of neutrinos in the NuMI beam is simulated in several steps as described in Ref. [5]. The NuMI beamline simulation uses the g4Numi code [35], which is based on a GEANT4 description of the geometry. The PPFX software package provides a neutrino flux prediction and uncertainties on the flux [36]. We need to estimate the number of kaons decaying at rest in the absorber, in order to simulate the LLP signal with the g4Numi beamline simulation.
The NuMI flux from the absorber leads to a negligible rate of neutrino interactions contributing to, e.g., MicroBooNE cross-section results [5,6]. The modeling of the flux of neutrinos from kaons decaying at rest in the absorber is however relevant for the LLP signal simulation. We therefore improve on this flux simulation by incorporating previous work from the MiniBooNE collaboration.
The MiniBooNE detector is located downstream from the MicroBooNE detector in the NuMI beam. The Mini-BooNE collaboration measured the ν µ flux from kaons decaying at rest at the NuMI absorber [34] using several methods and compared it to several predictions. They adopted the GEANT4 prediction of 0.085 ν µ produced per POT as their central value, with a 30% uncertainty taken from the range between simulations. The flux predicted by g4Numi is almost an order of magnitude smaller than the MiniBooNE central value and is not consistent with the 30% uncertainty. The g4Numi flux would yield a neutrino cross section measurement with the MiniBooNE detector that is inconsistent with expectations. We therefore use the MiniBooNE central flux value for the signal simulation, consistent with the procedure in Ref. [18].

VII. SIGNAL KINEMATICS
We generate the LLP signal using the flux of charged kaons that produce neutrinos, decaying them instead into an HNL or an HPS through the processes K → µN or K → πS. The two-body decay is isotropic in the kaon's rest frame, and the energy E LLP of the LLP is given by where m LLP is the LLP mass, m K the kaon mass, and m µ,π the mass of the daughter particle, which is either a muon or a pion. The HPS decays into the µ + µ − final state are simulated to reproduce the branching fraction as a function of HPS mass, with an isotropic decay distribution in the HPS rest frame. We simulate both charge conjugations of HNL decays, i.e., µ − π + and µ + π − final states, again with isotropic angular distributions. The angular distributions are then re-weighted for Dirac and Majorana HNL decays [37]. We assume a uniform time distribution of the proton beam interacting in the target, with a NuMI beam window of 9.6 µs. The two daughter particles are boosted to the lab frame using the LLP momentum vector. The time of the LLP decay is calculated from the kaon decay time and the time of flight of the LLP to the decay point. All times of flight are taken into account in the simulation.
The LLP is only retained if its momentum vector intersects the detector volume. A decay vertex within the detector volume is then chosen along the trajectory of the LLP. The exponential decay of the HPS flux is accounted for when selecting a decay vertex. The HPS lifetime is proportional to θ −2 [30]. The decay length of an HPS with m HPS > 2m µ and a mixing angle in the region of interest of 10 −7 < θ < 10 −9 is similar to the distance between the absorber and the MicroBooNE detector. Therefore some HPS will decay before reaching the detector, reducing the flux in the MicroBooNE detector.
For large values of θ 2 , only a small fraction of the HPS reach the detector before decaying, which restricts the upper reach of exclusion contours as a function of θ 2 . To derive these contours, we first define an effective angle θ eff which neglects the impact of decays before the detector. The event rate in the MicroBooNE detector is then ∝ θ 4 eff . In the final extraction of the limits, we correct for the early decays to obtain the limits as a function of the physical mixing angle θ 2 .
The exponential decay of the HNL flux is negligible for the mixing angles considered here as the HNL lifetime is much longer than the time needed to reach the MicroBooNE detector. The number of HNL decaying before reaching the detector is therefore neglected, and the final event rate is proportional to |U µ4 | 4 .

VIII. SIMULATION AND RECONSTRUCTION
Simulation and reconstruction are performed within the LArSoft [38] framework. The neutrino event generator GENIE [39] simulates the neutrino interactions on argon inside the cyrostat as well as interactions with the surrounding material. The GENIE configuration used in the simulation is found in Ref. [40], which describes a tuning of phenomenological parameters related to chargedcurrent quasi-elastic scattering and scattering on a pair of correlated nucleons, based on a fit to external data.
To obtain the response of the detector to ionization charge and scintillation light, the propagation through the detector material of secondary particles produced in the LLP decays or in the neutrino interactions is simulated by GEANT4 [41].
Cosmic rays crossing the detector in the same readout window are taken from data recorded outside of the beam window and then overlaid on the MC simulation. This overlay also addresses the need for simulating detector noise, as the simulated waveforms are combined with waveforms from data recorded during beam-off times.
The simulated samples and the data are reconstructed using the same algorithms. We use the Pandora pattern recognition framework [42] to combine hits, first clustered independently on each anode plane and then combined across anode planes to build particles reconstructed in 3D. The particles are arranged in a parent-daughter hierarchy based on the topology of the event and classified as track-like (muons, charged pions, and protons) or shower-like (electrons and photons).
Optical hits are constructed from PMT waveforms. Time-coincident optical hits from different PMTs are combined to form "flashes" that are attributed to a single interaction in the detector. The time of the flash, the associated location, and the amount of light are determined for each flash.

IX. DATA SETS
In total, we recorded data corresponding to 2.23 × 10 21 POT with the MicroBooNE detector exposed to the NuMI beam. For this paper, we analyse a subset of the data corresponding to 7.01×10 20 POT, which were taken in two different operating modes, forward horn current (FHC) and reverse horn current (RHC). The FHC data set was recorded during Run 1 in 2015-2016 and the RHC data set during Run 3 in 2017-2018. The CRT was fully installed for the second period, where it is used in the analysis. We thus analyze the two data sets separately to account for differences in neutrino flux, detector configuration, and CRT coverage.
The "beam-on" samples are defined by triggers that are coincident with the neutrino beam; they are used to search for the LLP signal. In addition, we use three samples for each each data period (FHC and RHC) that are designed to describe the background ("beam-off", "incryostat neutrino", and "out-of-cryostat neutrino").
The majority of beam-on events do not contain a neutrino interaction in MicroBooNE and are triggered by a cosmic ray. This source of background is modeled using a sample of events collected under identical trigger conditions but at times where no neutrino beam is present. These samples are referred to as "beam-off". The beamoff sample is scaled so that its normalisation corresponds to the number of beam spills of the beam-on sample. An additional scaling factor of 0.98 is applied to the beamoff sample. This factor takes into account that about 2% of all NuMI spills contain a neutrino interaction in the detector, and the remaining spills contain only cosmic rays [43].
The neutrino-induced background from the NuMI beam is modeled using two samples for each run period. The "in-cryostat ν" sample contains interactions of neutrinos with the argon inside the cryostat, and the "outof-cryostat ν" sample describes interactions with the detector structure and surrounding material. Both samples are generated with the MC simulation and normalized to the number of POT. The number of POT and the number of events passing the online-trigger before scaling factors are applied are summarized in Table I.
After matching the sample size to the number of POT of the beam-on data, the out-of-cryostat ν samples are scaled by an additional factor to ensure that the sum of the beam-off, in-cryostat and out-of-cryostat ν samples matches the data normalization within the NuMI timing window of 5.64-15.44 µs, as shown in Fig. 6. There is a small residual difference of ≈ 2% in the normalization of data relative to the sum of the background samples at this stage, as the scaling factors are derived with only a subset of the data.
Simulated data sets using the procedure described in Sec. VIII are used to evaluate the reconstruction and selection efficiency for the signal LLP decays. In total, we generate twelve samples in the mass range 246 ≤ m HNL ≤ 385 MeV and eight samples in the mass range  The beam-on data are compared to the sum of the beamoff, in-cryostat ν, and out-of-cryostat ν samples. The single bin at the start of the neutrino spill at 5.64 µs, where the expectation significantly exceeds the data, is caused by an effect in the timing structure of the beam that is not included in the simulation. The overall normalisation of the data and prediction agree within ≈ 2%.
212 ≤ m HPS ≤ 279 MeV. The spacing of the mass points takes into account the resolution. Additional points were generated at the edges of phase space, where limits could change more rapidly. The expected time of arrival of HNLs as a function of m HNL is shown in Fig. 7. The MicroBooNE NuMI trigger requires a flash in a timing window of 4.69-16.41 µs that covers the arrival times of the NuMI neutrino beam at the detector. More than 95% of HNLs arrive within the trigger window at low m HNL . Since the time of flight increases with m HNL , the trigger efficiency decreases to ≈ 55% for m HNL = 385 MeV and then quickly goes to zero. As we study a lower mass range for the HPS, this effect is not relevant for the HPS search.

X. CANDIDATE VERTEX IDENTIFICATION
The Pandora reconstruction algorithm groups objects into "slices" after removing cosmic-ray related hits. The slice identification (SliceID) process uses a combination of charge and light information to identify whether a slice corresponds to a neutrino interaction. Slices containing through-going, out-of-time, or stopping muons are removed. The reconstructed charge in the slice is required to be consistent with the location and intensity of the flash that has triggered the event. If this selection leaves ≥ 2 remaining slices in the event, a Support Vector Machine (SVM) is applied to calculate a "topological score" to select the most "neutrino-like" of the slices. This slice is then also used for the LLP search.
The LLP topology is characterized by a single vertex with two decay particles. The Pandora reconstruction frequently misplaces the location of the LLP decay vertex as the signal topology differs from the standard neutrino interactions (see Figs. 1 and 2). The vertex locations are therefore re-calculated for this analysis. All pairs of fitted tracks whose start or end points lie within a 3Ddistance of 5 cm of each other are combined to form signal candidates with a common vertex. If the end point of a track is placed at one of the new candidate vertices, the track direction is reversed. The vertex coordinates are calculated as the mean of the start locations of the two tracks. At this stage, there can be multiple vertices in an event. Approximately (60−70)% of signal events contain at least one vertex and (40 − 50)% contain exactly one vertex for the generated signal samples (Table II).

XI. KINEMATIC AND TOPOLOGICAL VARIABLES
We define several kinematic and topological variables to discriminate between signal and background. These variables are either associated to the slice or to the signal candidate. The slice-related variables use information calculated for the entire slice: • Multiplicity: The total multiplicity of objects in the slice, N tot , and the multiplicities of the objects classified as either tracks or showers (N tr or N sh ).
• Containment: Events are required to be contained inside the TPC's active volume by restricting the maximum and minimum coordinates of the start and end points of the objects within the slice to be 9 < x < 253 cm, −112 < y < 112 cm, 14 < z < 1020 cm. ( We also define the maximum and minimum extent of the slice for each coordinate, denoted as max (i) and min (i), with i = x, y, z. They are defined as the largest or smallest value of the start and end points of tracks in the slice in each dimension.
• Slice energy: The energy of all the objects of the slice, E sl , is reconstructed from the charge read-out by the collection planes. We expect that where E LLP is the energy of the LLP and m 1 and m 2 are the masses of the two decay particles. This assumes the decay particles are stable and their kinetic energy is fully measured in the TPC. The ionization energy contributed by secondary decays will increase E sl .
• Topological score: A Support Vector Machine (SVM) is applied to calculate a topological score that selects the slice that is most consistent with containing a neutrino interaction.
The remaining variables are directly related to the signal candidate.
• Proton likelihood: The log-likelihood particle identification score, S PID , is designed to utilize the calorimetric information to discriminate between protons and minimally ionizing particles (muons) [44]. A score of S PID = −1 indicates that a track is consistent with the proton hypothesis and S PID = +1 that a track is consistent with the muon hypothesis.
• Track length: We calculate the lengths t of the two tracks defining the candidate in the TPC. • Candidate four-momentum: The momenta of the particles associated to the two tracks produced by the LLP decay are determined using the length of the tracks and the continuous-slowing-downapproximation [45]. For the HNL decays the momenta and candidate mass are calculated assuming that the longer track is the muon and the shorter track is the pion. The momenta are summed to calculate the LLP candidate's four-momentum. We define β as the angle between the momentum of the candidate and the vector connecting the centre of the absorber to the candidate vertex. The angle α is the opening angle of the two tracks.

XII. EVENT SELECTION
The first stage of the final selection requires that the flash time coincides with the beam window, and only the flash with the largest number of associated photoelectrons is used. If there is a CRT hit within 1 µs of the flash for the Run 3 (RHC) sample, the event is identified as a cosmic ray and rejected. The 3D distance between the Pandora vertex and the closest tracks reconstructed using the CRT is required to be > 20 cm to remove events where charge from cosmic rays not in time with the beam flash is reconstructed as part of a candidate.
We require 2 ≤ N tot ≤ 4, N tr ≤ 3, and N sh ≤ 2. The event must be contained and the reconstructed energy E sl < 500 MeV. To remove neutrino interactions producing a proton, candidates where at least one track has a score S PID < −0.5 are removed. Candidates are required to have an opening angle cos α > −0.94. This removes events where a single cosmic-ray is split into two tracks with a large opening angle. Finally, the longest track must be shorter than 50 cm. Table III shows the combined efficiency of the selection requirements and vertex reconstruction relative to the SliceID requirements for background and signal. The background rejection is better for Run 3 compared to Run 1 as the CRT improves the rejection of cosmic rays. The CRT also vetoes particles produced in neutrino interactions that either enter or exit the detector, as expected for out-of-cryostat and in-cryostat neutrino interactions. The effect of the CRT veto on the signal is small, as the decay tracks are shorter. The signal efficiency relative to the SliceID requirements are in the range (36-45)%, rising with LLP mass. The total efficiency for an LLP decaying in the detector to be triggered, to be reconstructed and selected by the SliceID, and to pass the final event selection is in the range (13-30)%, again increasing with LLP mass.

XIII. BOOSTED DECISION TREE
We use the XGBoost gradient boosting library [46] to train Boosted Decision Trees (BDTs) that discriminate between the LLP signal and the background candidates passing the initial selection. trained with 30% of the selected in-cryostat ν sample and 50% of each of the generated signal samples. To improve performance of the BDT training, we require that > 90% of the hits in the slice are created by the LLP decay products. This procedure excludes mis-reconstructed signal events.
The background sample contains events where hits from overlaid cosmic events are mis-reconstructed as signal candidates. Therefore, cosmic-ray background is rejected by the BDT, even without training on a beam-off sample. In total, we use 21 BDT input variables: • slice energy E sl , • topological score, • maximum extent of the slice, max (i), and minimum extent, min (i), with i = x, y, z, • multiplicities N tot , N sh , and N tr , • candidate angle β, • candidate opening angle α, • candidate mass m LLP , • angles θ and φ of the longest track, as defined in the MicroBooNE coordinate system, and • length, the PID score S PID and the track score of the two tracks forming the candidate.
To compare data and background simulation, and to demonstrate the sensitivity to an LLP signal, we show eight of the more important variables in the BDT in Fig. 8. In general, we observe good agreement between the background prediction and the data across all variables, both in shape and normalization. The distributions shown are for Run 3 (RHC). The distributions for Run 1 (FHC) are similar with an increased cosmic-ray contribution (see Table III) as the CRT was not yet operational.
In Fig. 8, we overlay the signal distribution for a typical mass point, at m HNL = 304 MeV. As expected, the average reconstructed invariant mass is centered around this value. The slice energy peaks at E sl ≈ 105 MeV, about 20 MeV above the expected value from Eq. 4. This shift is related to the ionization energy deposited by Michel electrons. The slice energy is the variable with the best sensitivity to signal in the BDTs. The candidate opening angle is expected to peak at cos α = −0.34 for this mass. The candidate angle β with respect to the direction from the absorber is peaked at cos β = 1 for a large fraction of the HNL signal candidates.
The lengths of the tracks are measures of the momenta of the two particles, which depend on the decay kinematics. The direction of the longer decay track is included as the angles θ and φ of the track in the detector coordinate system. The direction (in radians) from the absorber cor-responds to θ = 2.20 and φ = 1.15. The angle φ helps to reject cosmic rays, which are aligned with φ = ±π/2. Figure 9 shows the distributions of the BDT scores trained for three different HNL and HPS masses. The BDTs offer strong rejection against background candidates from neutrino interactions. Background candidates from mis-reconstructed cosmic-ray tracks, as found in the beam off sample, are also rejected with high efficiency. The discrimination between signal and background improves slightly with LLP mass due to the higher energy of the decay particles, leading to improved kinematic reconstruction.
In Fig. 9, the full range of BDT scores lies in the range [0,1]. The final score distribution is shown after a transformation using the inverse of the logistic function has been applied, which maps the score to a range [−∞,∞]. We only consider candidates with BDT score > 0 after the transformation, corresponding to > 0.5 in the original distribution, since this region contains > 90% of the signal and hence dominates the sensitivity.

XIV. SYSTEMATIC UNCERTAINTIES
Uncertainty sources are considered for the background samples and signal samples by applying variations that modify the BDT score distributions. For the simulated in-cryostat ν background sample, we consider the impact of the flux simulation, cross-section modeling, hadron interactions with argon, and detector variations: • Flux Simulation. Uncertainties on the neutrino flux arise primarily from the rates and kinematics of hadron production in the beamline. The PPFX package is used to estimate these uncertainties. Each of the parameters in the constrained flux prediction is simultaneously sampled within its estimated uncertainties to produce alternative flux predictions. Flux uncertainties also include variations in the beamline conditions, e.g., changes in horn current, horn position, and beam spot location. However, these beamline condition uncertainties are neglected, as they have been shown to be small [5,6] and were therefore not reassessed.
• Cross-section Modeling. We assess the uncertainties due to the modeling of neutrino crosssections by varying 44 of the parameters used by the GENIE generator. A full discussion can be found in Ref. [40].
• Hadron Interactions. Hadrons produced in neutrino interactions interact strongly, affecting their propagation through argon. The description of hadron interactions therefore impacts the event reconstruction and the description of the neutrino background. These uncertainties are assessed using GEANT4Reweight [47] by considering variations in the GEANT4 cross-section model for charged pions and protons.
• Detector Modeling. Uncertainties arising from detector modeling are estimated by re-simulating the detector configuration using the same generated input event sample. The variations include effects due to the light simulation, reduction in the light yield, increase of the Rayleigh scattering length, and attenuation of the light in the argon.
Uncertainties related to charge reconstruction are assessed using data-driven modifications of the the waveforms on the TPC wires, as discussed in Ref. [48]. Additional variations due to the space charge mapping and ion recombination model are simulated and assessed separately.
We expect only a small number of background events in the signal region of the BDT distribution due to the high purity of the event selection. We therefore extrapolate detector-modeling uncertainties from higher statistics regions of the distribution to the signal region, assuming they are constant.
The beam-off sample is taken from data and therefore has no associated systematic uncertainties other than the statistical fluctuations in the sample. The contribution of the out-of-cryostat sample to the final sample is small. The out-of-cryostat sample normalization (see Section IX) is set as an unconstrained parameter in the final fits and found to be negligible.
The impact of the systematic uncertainties on the BDT score distribution is shown in Fig. 10 for the Run 1 (FHC) and Run 3 (RHC) background samples, with the BDT trained for a signal of m HNL = 304 MeV. The fractional uncertainties are given relative to the total background prediction, which is the sum of beam-off, in-cryostat ν, and out-of-cryostat ν samples.
The dominant uncertainty in the signal region at high BDT scores is due to the statistical uncertainty of the background samples, which is a consequence of the high purity of the signal selection. Detector modeling uncertainties are ≈ (10 − 20)%, neutrino flux and cross-section uncertainties are each ≈ (5 − 10)% for most bins, while all other uncertainties are negligible. The systematic uncertainties are separately evaluated for all signal training points (signal masses and FHC/RHC), with consistent results. As expected, the sensitivity of the final result is dominated by the the RHC data set.
The dominant contribution to the systematic uncertainty on the signal sample arises from uncertainties on the rate of kaon production at rest in the NuMI absorber. This normalisation uncertainty is taken from the evaluation of the MiniBooNE collaboration to be ±30% [34], as discussed in Section VI.
The hadron interaction uncertainty will mainly affect the modeling of pions produced in the HNL decay. The relative uncertainties on the BDT distribution are ≈ 2% for hadron interaction modeling and < 15% for detector modeling. The longest track is assigned to be the muon in HNL decays; the systematic uncertainty due to this choice is also negligble. Finally, a 2% uncertainty on the number of POT delivered is estimated from the uncertainty on the beamline toroidal monitor measurement [35], which is taken as fully correlated uncertainty for all samples.

XV. RESULTS
We apply the BDTs trained for each signal mass point to the data and the simulation. In Figs. 11 and 12 we show a comparison of the resulting BDT distributions for the data to the sum of the background predictions, using a selection of representative signal mass points. In events with more than one candidate, we retain only the candidate with the most signal-like BDT score. The background predictions and the data are in good agreement across both of the run periods studied.
The BDT score distributions are used as input to a modified frequentist CL s calculation [49,50] to set upper limits on the signal strength for each model and mass point. Test statistics are constructed from the log-likelihood ratio (LLR) of the signal-plus-background (S+B) and background-only (B) hypothesis for each LLP mass. The BDT distributions for each run period (FHC and RHC) enter the limit setting as separate channels before their likelihoods are combined.
The systematic uncertainties on the background and signal predictions are taken into account using Gaussian priors. The signal and background predictions are separately fitted to the data distributions under the background-only and signal-plus-background hypotheses. The systematic uncertainties are allowed to vary within their defined priors to maximize the respective likelihood functions. This reduces the effect of the systematic uncertainties on the sensitivity of the search [51].
The background predictions are fitted separately for the two run periods, with uncorrelated systematic uncertainties. We find that this assumption about the correlations has negligible impact on the result. The flux systematic uncertainties on the signal are taken to be correlated between the two data periods. The confidence levels are calculated by integrating the LLR distributions, which are derived using pseudoexperiments, under either the S+B (CL s+b ) or the B-only hypotheses (CL b ). The excluded signal rate is defined as the scaling of the signal strength for which the confidence level for signal reaches CL s = CL s+b /CL b = 1 − 0.9. The observed and median expected 90% CL limits on |U µ4 | 2 are shown for each HNL mass point in Table IV and in Fig. 13, and the corresponding limits on θ 2 eff for the HPS model in Table V. The 1-and 2-standard deviation intervals cover the range of expected limits produced by 68% and 95% of background prediction outcomes around the median expected value. The observed limits are contained in the 1-standard-deviation interval for all mass points with the exception of m HNL = 371. 5   385.0 MeV, and for m HPS = 215 MeV, where the observed limit lies within 2 standard deviations. We use a linear interpolation between the mass points when drawing contours.
We derive the HNL limits assuming that HNLs are Majorana particles. For a Dirac HNL, only decays to the charge conjugated final state µ − π + are allowed in K + decays. The expected number of decays is therefore a factor of two smaller for the same |U µ4 | 2 value. The limits for Dirac HNLs are calculated from the Majorana limit by applying a factor of √ 2 to account for the reduced decay rate, since the difference due to angular distributions of the decay is found to be negligible. Table V shows the values of θ 2 that correspond to the lower and upper bounds of the excluded region for each mass m HPS , where the upper boundary is due to the short lifetime of the HPS. For the highest mass point at m HPS = 279 MeV, the number of HPS reaching the detector before decaying is too small to derive a limit for any of the θ 2 values within the excluded θ 2 eff range.

XVI. COMPARISON WITH EXISTING LIMITS
In Figs. 14 and 15, we compare the observed limits to the existing experimental limits in similar regions of parameter space for both models. The results extend MicroBooNE's sensitivity to |U µ4 | 2 by approximately an order of magnitude compared to the previous MicroBooNE HNL result [8].  BDT Score (m HPS =215 MeV) Like the MicroBooNE detector, the T2K [9] and NuTeV detectors [10] were located in a neutrino beamline. The PS191 experiment [11,12] at CERN was specifically designed to search for massive decaying neutrinos. The NA62 [13] and E949 [14] collaborations performed a peak search for HNLs in kaon decays. The muon spectrum measured in stopped K + → µ + ν decays (K 2µ ) has also been used to set limits on HNLs [15,16].
In the mass range 300 < m HNL < 385 MeV, this search has similar sensitivity as the NA62 experiment [19]. The E949 [14], PS191 [12] and T2K [9] limits are stronger across the range 300 < m HNL < 385 MeV. The T2K collaboration provides no limit point for masses above 380 MeV. Here, the MicroBooNE limit is of equal or greater sensitivity than the NA62 result.
For the HPS model, we constrain a region of parameter space for 212 < m HPS < 275 MeV not previously excluded by any dedicated experimental search. The existing limits in this region are reinterpretations of decades old CHARM [24], LSND [27], and PS191 [26] measurements, performed by authors outside the respective collaborations without access to the original experimental data or MC simulation. Reinterpretations depend on external beamline, flux, and detector simulations. If the signal topology differs from the original selection criteria, the results also depend on estimated detection efficiencies. In the case of the CHARM experiment, the more recent sensitivity estimate in Ref. [24] disagrees by nearly an order of magnitude from the estimate in Ref. [25].  [53], KEK [16], NA62 [19], E949 [14], PS191 [12], T2K [9] and NuTeV [10] collaborations.

XVII. SUMMARY
We present a search for long-lived particles (LLP) using NuMI beam data corresponding to 7.01 × 10 20 POT recorded with the MicroBooNE detector. The results are interpreted within two models, where the LLP is either a heavy neutral lepton (HNL) or a Higgs portal scalar (HPS). The LLPs are assumed to be produced by K + mesons decaying at rest in the NuMI absorber. The signature in the MicroBooNE liquid-argon detector are HNL decays into µ ± π ∓ pairs or HPS decays into µ + µ − pairs.
The main sources of background are neutrino and cosmic-ray interactions, where the majority of the neutrino events are from charged-current muon neutrino interactions. To reject background, we select data recorded in-time with the NuMI beam and consistent with the LLP signatures in the liquid argon. The LLPs originating in the NuMI absorber enter the detector from a different direction than the majority of beam neutrinos. The decay products also have a fixed energy for a given LLP mass.  [24], LSND [27], and PS191 [26] measurements. In other mass ranges, limits are from a MicroBooNE search for the e + e − final state [18] (at the 95% CL), and from searches by the NA62 [19,20] and E949 collaborations [21] for charged kaon decays to pions and an HPS. The LHCb collaboration performed two searches for an HPS with short lifetime, which would be produced and subsequently decay within the detector [22,23]. The joint coverage of the LHCb result is shown at the 95% CL. These kinematic properties are used to discriminate signal from background. The combined reconstruction and selection efficiencies for LLP decays are between 13% and 30%, increasing with LLP mass. To further improve discrimination between signal and background, we train and apply a BDT with 21 input variables for each mass point. No significant excess is observed in the BDT score distributions. In the absence of signal, we employ the modified frequentist CL s method to derive limits on the model mixing parameters |U µ4 | 2 and θ 2 . All limits are presented at the 90% CL.
We also constrain the scalar-Higgs mixing angle θ by searching for HPS decays into µ + µ − final states, excluding a contour in the parameter space with lower bounds of θ 2 < 31.3 × 10 −9 for m HPS = 212 GeV and θ 2 < 1.09 × 10 −9 for m HPS = 275 GeV. These are the first constraints in this region of the θ 2 -m HPS parameter space from a dedicated experimental search. It is also the first search in this mass range using a liquid-argon TPC.
the United Kingdom Research and Innovation; the Royal Society (United Kingdom); and The European Union's Horizon 2020 Marie Sk lodowska-Curie Actions. Additional support for the laser calibration system and cosmic ray tagger was provided by the Albert Einstein Center for Fundamental Physics, Bern, Switzerland. For the purpose of open access, the authors have applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising from this submis-sion. We also acknowledge the contributions of technical and scientific staff to the design, construction, 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.
We particularly acknowledge the many contributions of Salvatore Davide Porzio to devel-oping this research direction on MicroBooNE.