Search for a Higgs portal scalar decaying to electron-positron pairs in the MicroBooNE detector

We present a search for the decays of a neutral scalar boson produced by kaons decaying at rest, in the context of the Higgs portal model, using the MicroBooNE detector. We analyze data triggered in time with the Fermilab NuMI neutrino beam spill, with an exposure of $1.93\times10^{20}$ protons on target. We look for monoenergetic scalars that come from the direction of the NuMI hadron absorber, at a distance of 100 m from the detector, and decay to electron-positron pairs. We observe one candidate event, with a Standard Model background prediction of $1.9\pm0.8$. We set an upper limit on the scalar-Higgs mixing angle of $\theta<(3.3-4.6)\times10^{-4}$ at the 95% confidence level for scalar boson masses in the range $(100-200)$ MeV$/c^2$. We exclude at the 95% confidence level the remaining model parameters required to explain the central value of a possible excess of $K^0_L\rightarrow\pi^0\nu\bar{\nu}$ decays reported by the KOTO collaboration. We also provide a model-independent limit on a new boson $X$ produced in $K\rightarrow\pi X$ decays and decaying to $e^+e^-$.

We present a search for the decays of a neutral scalar boson produced by kaons decaying at rest, in the context of the Higgs portal model, using the MicroBooNE detector. We analyze data triggered in time with the Fermilab NuMI neutrino beam spill, with an exposure of 1.93 × 10 20 protons on target. We look for monoenergetic scalars that come from the direction of the NuMI hadron absorber, at a distance of 100 m from the detector, and decay to electron-positron pairs. We observe one candidate event, with a standard model background prediction of 1.9 ± 0.8. We set an upper limit on the scalar-Higgs mixing angle of θ < (3.3 − 4.6) × 10 −4 at the 95% confidence level for scalar boson masses in the range (100 − 200) MeV/c 2 . We exclude, at the 95% confidence level, the remaining model parameters required to explain the central value of a possible excess of K 0 L → π 0 νν decays reported by the KOTO collaboration. We also provide a model-independent limit on a new boson X produced in K → πX decays and decaying to e + e − . The Higgs portal model [1] is an extension to the standard model in which an electrically neutral real singlet scalar boson (S) mixes with the Higgs boson with mixing angle θ. Through this mixing, S acquires a coupling to standard model fermions proportional to sin θ and their Yukawa couplings with the Higgs boson. For masses between twice the electron mass and twice the muon mass, and assuming that there are no new dark sector particles lighter than half its mass, S will decay to electronpositron pairs with partial width [2] Γ = θ 2 m 2 e m S 8πv 2 where m S is the scalar boson mass, m e is the electron mass, and v is the Higgs field vacuum expectation value. For these masses, S can be produced from a kaon twobody decay in association with a pion, with the dominant production process being a penguin diagram with a top quark running in the loop. The partial width of the production process is [2] Γ where m K is the kaon mass, m π is the pion mass, m t is the top quark mass, V td and V ts are the elements of the CKM matrix, and λ is the Källen Lambda function.
In 2019, the KOTO collaboration reported [3] the observation of four K 0 L → π 0 + invisible decay candidates, a rate 2 orders of magnitude more frequent than the standard model prediction for K 0 L → π 0 νν decays. One candidate was rejected due to upstream veto activity, but three candidates remained as unexplained. In a recent publication [4], they have reevaluated their background expectation to 1.22 ± 0.26 counts, and the statistical significance of the observed data has reduced to a p value of 0.13. The Higgs portal model could explain [5][6][7] any excess in the KOTO dataset. The value of θ for a central value of 1.78 counts (θ KCV ) is ≈ (4 − 5) × 10 −4 over the m S range of (100−200) MeV/c 2 . The E949 collaboration excludes [8] θ KCV for m S < 120 MeV/c 2 , and the NA62 collaboration excludes it [9] for m S > 160 MeV/c 2 .
This Letter presents the first search for beyond the standard model (BSM) electron-positron pair production in a liquid argon time projection chamber using the Mi-croBooNE detector, and is the second search for BSM physics in MicroBooNE following a search for heavy neutral leptons [10]. We use the results of this search to exclude, at 95% confidence level, the remaining Higgs portal model parameter space where θ KCV has not been excluded.
The MicroBooNE experiment is primarily designed for neutrino scattering measurements in Fermilab's Booster Neutrino Beam (BNB) [11,12]. The detector sits just below the surface, and comprises an 85 ton liquid-argon time projection chamber (TPC) with active dimensions of 2.6 m along the drift direction (horizontal and perpendicular to the beam axis; x coordinate in the detector reference frame), 2.3 m in the vertical direction (y coordinate), and 10.4 m along the direction parallel to the BNB direction (z coordinate). Charged particles traversing the argon produce ionization electrons and scintillation light. Drifted ionization electrons are recorded by three wire planes with 3 mm pitch that are oriented at 60 • rotations relative to each other. An array of 32 eight-inch photomultiplier tubes (PMTs) dis- tributed behind the wire planes provides timing information for scintillation signals produced inside the cryostat. Part way through the detector operations, a cosmic ray tagger (CRT) system [13] was installed, with four walls of plastic scintillator panels situated along the top, bottom, and long sides of the cryostat that provide timing coincidence signals for some cosmic rays that enter the TPC.
In addition to being on the BNB beam line, the Micro-BooNE detector is also situated at a distance of 680 m and 8 • off axis from the target of Fermilab's Neutrinos at the Main Injector (NuMI) neutrino beam [14], which we have previously used to measure electron neutrino interactions [15]. A schematic diagram of the detector position within the beam line is presented in Fig. 1. The Main Injector delivers 120 GeV protons that impact the graphite target, producing secondary hadrons. A system of electromagnetic horns focuses charged particles either toward or away from the beam axis, depending on the horn polarity. In forward (reverse) horn current mode, positively (negatively) charged mesons are bent toward the beam axis to produce a beam mostly of neutrinos (antineutrinos) from the meson decays. A 675 m long helium-filled decay volume is situated downstream of the target and horn system, at the end of which is a 5 m deep hadron absorber. Any surviving hadrons will be stopped in the absorber, and may produce secondary mesons including K + which will decay at rest. The absorber is at a distance of ≈ 100 m from the MicroBooNE detector and at an angle of ≈ 125 • with respect to the BNB direction, such that any particles that enter the detector from the absorber are entering in the opposite direction compared to most neutrino interactions seen by the detector. We exploit the unique decay signature of scalar bosons produced by K + decaying at rest (KDAR) in the NuMI hadron absorber to search for evidence of the Higgs portal scalar model.
We analyze 0.92 × 10 20 protons on target (POT) of exposure during run 1 of MicroBooNE's operations (during 2015-2016), and 1.01 × 10 20 POT of run 3 data (2017-2018). During the run 1 dataset period, the NuMI beam operated in forward horn current mode, and reverse horn current mode was used during the run 3 period. The CRT had been fully installed by the run 3 period and we use its information in the analysis of data from that period. The beam-on data is read out from the detector (an "event") when there is a NuMI beam spill timing signal sent by the Fermilab accelerator complex. An online trigger is employed to record only those events that pass optical trigger criteria based on the total integrated charge summed over all PMTs in a 100 ns window. This trigger requires at least one PMT to produce a signal in time with the beam, and the integrated charge has to be above a configurable photoelectron threshold.
To estimate the cosmic-induced backgrounds, we record a dataset of events produced out of time with both beams that employs the same trigger thresholds, called the beam-off dataset. In addition, there is an unbiased dataset of out-of-time events for which the trigger is not applied. This unbiased dataset forms the basis of the simulated data. The hit pattern of simulated signal decays or background neutrino interactions are overlaid on top of the unbiased data on an event-by-event level, which allows the cosmic contamination of signal or neutrino interactions to be estimated.
To simulate the scalar boson signal, we use the g4numi program [16] which employs a GEANT4 [17] simulation of the NuMI beam line to produce the position and timing distribution of KDAR in the NuMI hadron absorber. For the absolute rate, we use the MiniBooNE estimate [18] of 0.085 muon neutrinos from KDAR in the NuMI hadron absorber per POT. The scalars are emitted isotropically from the kaon decay positions, and the scalar's velocity and lab-frame lifetime is used to distribute the scalar decay position, keeping only those that decay within the detector active volume. The electron-positron pair is simulated isotropically in the rest frame, and boosted by the scalar's momentum.
The g4numi program is also used to simulate the flux of neutrinos that intersect the detector, which produce the other component of the background to this search. We use the same PPFX [16,19] package as used by the MIN-ERvA [20] and NOvA [21] collaborations to correct the central value neutrino flux prediction and provide flux uncertainties. We use GENIE [22] to calculate the neutrino interaction cross sections and final state kinematics, in which the models for charged-current quasielastic scattering and scattering on a pair of correlated nucleons have been tuned based on data from T2K [23]. Neutrino interactions are simulated both inside the cryostat and outside where secondary products enter the detector. For both signal and background simulations, the decay or interaction products are propagated through a GEANT4 simulation of the detector. The response of the detector to both light and charge is simulated.
All three types of data (beam on, beam off, and simulated) are processed through the same chain of reconstruction algorithms. The optical reconstruction uses the PMT waveforms to produce "flashes" of coincident PMT hits. For the TPC information, we apply a two-dimensional deconvolution of the signal waveforms on the wires within each plane [24,25]. Hits are formed from a Gaussian peak finding algorithm applied to the wire waveform. The Pandora framework [26] uses particle flow algorithms to cluster the hits of a single plane, and, then, match clusters across planes into three-dimensional reconstructed objects, which Pandora classifies as "tracks" or "showers" based on a multivariate classifier score. Pandora also "slices" the event into groups of reconstructed objects that it considers to be independent interactions (either cosmogenic or beam induced) and removes well-identified cosmic slices. For any remaining slices, a flash-matching algorithm is applied to produce a PMT hit hypothesis using the reconstructed objects in the slice. The algorithm attempts to match the hypothetical PMT hit distribution with the observed flashes in the beam timing window, and calculates a χ 2 value for the best match. The best matching slice of these remaining slices is labeled as the neutrino slice.
To preselect decay candidates in the event, we use the neutrino slices. The slice has to be matched to a PMT flash with a time of [5.8, 16.8] µs within the 20 µs PMT readout window (where the NuMI prompt neutrino spill produces flashes in the range [6.1, 15.7] µs), and the flashmatching χ 2 has to be less than 10. An additional selection is imposed on the run 3 data, requiring that events cannot have a CRT hit in coincidence with the beam timing. The total number of objects in the slice has to be ≤ 5, and of these, a maximum of 4 can be labeled as tracks. For all possible pairs of objects in the slice, the minimum distance between the object vertices (for reconstructed tracks, both start or end positions, and for showers, only start positions) is calculated. If this distance is less than 5 cm, a decay vertex is produced at the midpoint between the object vertices with the minimum separation. The position of the decay vertex has to be reconstructed within the active volume of the detector. Slices with more than two objects could conceivably form multiple decay candidates, all of which are preselected and passed through the boosted decision tree (BDT) selection.
We apply two different BDTs to the preselected candidates: one trained against cosmic backgrounds and one trained against neutrino interactions simulated inside the cryostat. Each BDT is trained separately over the run 1 events and run 3 events, i.e., there are four BDTs in total. We split the run periods because the use of the CRT in run 3 and the differences between forward and reverse horn current operations can change the topologies and properties of the background distributions that the BDTs are trained against. We use xgboost [27] to train and apply the BDTs. We train the BDTs on ten input variables each. Nine of the ten input variables are the same for the cosmic-focused and neutrino-focused BDTs. These are (1) the opening angle between the two reconstructed objects; (2) the opening angle in the plane transverse to the hadron absorber direction from the detector center; (3,4) the two angles between the two objects and the hadron absorber direction; (5) the Pandora track or shower score of the larger of the two objects (when ordered by number of hits); (6) the number of hits of the larger object; (7) the total number of hits contained in other objects in the slice, not including the two objects that form the decay candidate; (8) the maximum y coordinate, relative to the decay vertex position, of shower start positions or track start or end positions, for any other objects in the slice; and (9) the minimum z coordinate, relative to the decay vertex position, of shower start positions or track start or end positions, for any other objects in the slice. The last two variables are treated as "missing" within xgboost if the slice contains only two objects. The tenth input variable of the cosmic-focused BDT is the length of the larger object. The tenth input variable of the neutrinofocused BDT is the number of tracks in the slice. For all input variables and output BDT score distributions we observe good data-simulation agreement in a control region of data with an early flash time (as the scalar boson signal is delayed by ≈ 600 ns with respect to the prompt neutrino interactions due to time-of-flight differences).
The BDTs are trained on a signal simulation where each decay is of a scalar boson with m S uniformly chosen between 100 and 200 MeV/c 2 , in order to reduce the dependence of the BDTs on m S in the range where θ KCV has not been excluded. The candidates used in the training have to be well reconstructed, with the cosmic contamination of each object below 10%, and the reconstructed vertex and directions close to the generated values. The neutrino-focused BDT is trained against 10% of the simulated statistics of neutrino interactions in the cryostat, with the other 90% along with all the out-ofcryostat simulated interactions used for the sensitivity and limit calculations. Each reconstructed object used in the neutrino background training sample is required to have cosmic contamination below 10%, similar to the signal sample. The cosmic-focused BDT is trained against beam-off candidates that fail the flash-matching χ 2 requirement.
To select candidate decays, we require the two BDT scores to be above a minimum score. We choose the four minimum scores that maximize the sensitivity of the selection to the model parameter θ for m S = 100 MeV/c 2 , as we expect even better sensitivity at higher masses. The 95% confidence level (C.L.) sensitivity and limit are calculated for a single-bin counting experiment with the modified-frequentist CLs method, using the RooStats statistical package [28], including systematic uncertainties as constrained Gaussian nuisance terms.
We consider several classes of systematic uncertainty. We include the simulation statistical uncertainty and beam-off data statistical uncertainty as uncertainties for the model prediction. The flux normalization uncertainty on the signal model is set to 30% as used by Mini- BooNE [18]. The uncertainty on the background neutrino cross section modeling is evaluated by reweighting events using tools included with GENIE [29]. Interaction physics model parameters (both in GENIE and GEANT4) are varied multiple times within their 1 standard deviation uncertainties, and a weight is calculated for each simulated event between the central value and the modified model. The uncertainty on the event count in the selection is calculated from the standard deviation of weighted event counts across the variations. A similar procedure is followed to estimate the uncertainty on the background neutrino flux model due to hadron production uncertainties, using PPFX [19]. The flux uncertainty due to the beam line model (including focusing) is negligible compared to the hadron production uncertainty [15] and is not included. Systematic uncertainties due to the modeling of the detector are evaluated through modified simulations varying parameters of the detector model. They are estimated to be 70% for the neutrino background simulation (dominated by the low statistics of simulated neutrino events in the signal region after selection) and 5% for the signal simulation, taken to be the relative differences of event yields in the signal region between the central value and ten detector model variations summed in quadrature. The first five detector model variations are (1) uncertainties in the space charge mapping [30], (2) the ion-electron recombination model, (3) a decrease in light yield, (4) an increase in the Rayleigh scattering length, and (5) changing the light attenuation between the anode and cathode sides. We also modify the simulated TPC wire waveform amplitudes and widths. The sizes of the modifications are characterized in five dimensions based on hit positions, track angles with respect to the wires, and energy deposited per unit length. The modification sizes are estimated by comparing orthogonal data samples rich in protons and cosmic muons to the central value simulation. These five wire modifications (each dimension independently) are then applied to the signal and background simulation and used to extract the event yield variation. Although the uncertainty on the detector model for the background prediction is large, the final result is statistics-limited, and this uncertainty has minimal impact with respect to repeating the analysis with zero detector uncertainty. The uncertain-  ties in the signal region after the optimal BDT selection are given in Table I.
After applying the BDT selection, the number of events expected for each background contribution and for several signal definitions are shown in Table II. The table also presents the estimated signal selection efficiency. The total expected background-only prediction is 1.9±0.8 candidate events.
In the beam-on dataset, we observe two candidates in the signal region. We reject one candidate because its flash time of 5.84 µs lies in the window between the start of the selection time (5.8 µs) and the start of the neutrino interactions (6.1 µs), making it an obvious cosmic background interaction. This post-selection cut only affects the cosmic background, reducing the cosmic acceptance by 2.7%, with negligible effect on the sensitivity. When we manually inspect the TPC readout of the other candidate event, the two objects have the characteristics of a proton and a photon, and so, it is likely to be a neutrinoinduced background.
With one observed event we set the 95% C.L. upper limit on the Higgs portal model presented in Fig. 2. The observed and expected limits for several scalar boson masses are enumerated in Table III and, for a wider range of masses, in the Supplemental Material [31]. The upper limit is compared with θ KCV , along with other experi-   [4]. The limits for E949 [8] and NA62 [9] are published by the collaborations, whereas the CHARM [6] and LSND [32] limits are reinterpretations of other searches. mental limits, in Fig. 3. In Fig.4, we present our result as a model-independent limit on a new boson X produced in K + → π + X decays, and decaying to e + e − pairs, for X masses in the range [100, 210] MeV/c 2 .
The limit presented in this publication rules out the remaining Higgs portal model parameter space required to explain the central value of a mild excess in KOTO at the 95% confidence level. Our limit is the most constraining for m S ≈ (120 − 160) MeV/c 2 and is directly derived from our own experimental data. The previous most stringent constraints in this range were reinterpretations of decades-old CHARM [6,33] and LSND [32] measurements, performed recently by independent authors without access to the raw experimental data. We have ≈ 2 × 10 21 POT of as-yet unprocessed NuMI data along with ≈ 1×10 21 POT of currently blinded BNB data that we will analyze in the future and expect improved sensitivities [2]. This document was prepared by the MicroBooNE collaboration using the resources of the Fermi National Ac-