Search for an anomalous excess of inclusive charged-current $\nu_e$ interactions in the MicroBooNE experiment using Wire-Cell reconstruction

We report a search for an anomalous excess of inclusive charged-current (CC) $\nu_e$ interactions using the Wire-Cell event reconstruction package in the MicroBooNE experiment, which is motivated by the previous observation of a low-energy excess (LEE) of electromagnetic events from the MiniBooNE experiment. With a single liquid argon time projection chamber detector, the measurements of $\nu_{\mu}$ CC interactions as well as $\pi^0$ interactions are used to constrain signal and background predictions of $\nu_e$ CC interactions. A data set collected from February 2016 to July 2018 corresponding to an exposure of 6.369 $\times$ 10$^{20}$ protons on target from the Booster Neutrino Beam at FNAL is analyzed. With $x$ representing an overall normalization factor and referred to as the LEE strength parameter, we select 56 fully contained $\nu_e$ CC candidates while expecting 69.6 $\pm$ 8.0 (stat.) $\pm$ 5.0 (sys.) and 103.8 $\pm$ 9.0 (stat.) $\pm$ 7.4 (sys.) candidates after constraints for the absence (eLEE$_{x=0}$) of the median signal strength derived from the MiniBooNE observation and the presence (eLEE$_{x=1}$) of that signal strength, respectively. Under a nested hypothesis test using both rate and shape information in all available channels, the best-fit $x$ is determined to be 0 (eLEE$_{x=0}$) with a 95.5% confidence level upper limit of $x$ at 0.502. Under a simple-vs-simple hypotheses test, the eLEE$_{x=1}$ hypothesis is rejected at 3.75$\sigma$, while the eLEE$_{x=0}$ hypothesis is shown to be consistent with the observation at 0.45$\sigma$. In the context of the eLEE model, the estimated 68.3% confidence interval of the $\nu_e$ hypothesis to explain the LEE observed in the MiniBooNE experiment is disfavored at a significance level of more than 2.6$\sigma$ (3.0$\sigma$) considering MiniBooNE's full (statistical) uncertainties.


I. INTRODUCTION
Neutrino flavor oscillation is one of the few observations in particle physics for evidence of physics beyond the standard model. While the majority of neutrino oscillation data can be successfully explained by a threeneutrino framework [1], the exact mechanism for neutrinos to acquire their masses remains a puzzle. In addition, the fact that the mass of the electron neutrino is at least five orders of magnitude smaller than that of the electron [2] is an interesting mystery. The possible existence of additional neutrino flavors may provide a natural explanation of the small neutrino mass [3]. Constrained by precision electroweak measurements [4], these additional neutrinos are expected to be sterile [5] not participating in any fundamental interaction of the standard model. Despite the fact that there is no known mechanism to detect them directly, they may be indirectly observed in neutrino oscillation experiments where sterile neutrinos could mix with the three active neutrinos and affect the way they oscillate.
Because of strong theoretical motivations, there are many dedicated programs searching for sterile neutrinos. While most of the results are consistent with the three-neutrino framework without sterile neutrinos (see Ref. [6,7] among others), there are several experimental anomalies suggesting the existence of an eV-massscale sterile neutrino: i) the observation that calibrated * microboone info@fnal.gov ν e sources ( 51 Cr for GALLEX [8] and BEST [9], 51 Cr and 37 Ar for SAGE [10]) produced lower rates of measured ν e than expected in the three-neutrino framework, which could be explained by ν e disappearance considering light sterile neutrinos; ii) the reactor anti-neutrino anomaly [11], which suggests that the observed deficit in the measuredν e events relative to the expectation based on the recent reactor anti-neutrino flux calculations [12,13] could be explained byν e disappearance considering light sterile neutrinos, although there are recent experimental measurements [14,15] and improved flux calculations [16,17] that disfavor this explanation; iii) the Neutrino-4 [18] anomaly, which suggests reactor ν e oscillation at a few meters; and iv) the LSND [19] and MiniBooNE [20,21] anomalies, which suggest ν e appearance from ν µ to ν e oscillations considering light sterile neutrinos. However, there are significant challenges in explaining all available experimental results with a sterile neutrino oscillation model in a global fit [22].
The MiniBooNE experiment observes an anomalous excess of electromagnetic (electron-or photon-like) events in the data above the prediction for reconstructed neutrino energies below 800 MeV, which is commonly referred to as the MiniBooNE low-energy excess (LEE). The interpretation of the MiniBooNE LEE observation is limited by the capability of its detector technology, which primarily uses the pattern of Cherenkov rings to differentiate between muons and electrons. However, the Cherenkov technology does not perform well in separating an electron from a photon leading to the possibility that the LEE may come from an excess of pho-ton background rather than from the ν e charged-current (CC) interactions. The Liquid Argon Time Projection Chamber (LArTPC) [23][24][25][26] has superb imaging capabilities through the combination of both tracking and calorimetry in a fully active volume. Throughout the years since its inception, many LArTPCs ranging in size from hundreds of liters to hundreds of cubic meters have been constructed and operated for neutrino experiments [27][28][29][30][31][32][33]. Although electrons and photons with energy greater than tens of MeV both induce electromagnetic (EM) showers in LArTPCs, electrons can be differentiated from photons in neutrino interactions through identifying a gap between neutrino and EM shower vertices, topological pattern recognition, and reconstructed ionization energy loss per unit length (dE/dx) measurement [34]. This improves the identification of ν e CC interactions and will enable sensitive measurements such as searches for ( − ) ν µ → ( − ) ν e oscillations. LArTPCs can also reconstruct and identify hadrons at lower energies than Cherenkov detectors, which can only detect charged particles with energies above their respective Cherenkov thresholds. The extended capability of LArTPCs to detect these lower energy hadrons enables the study of various exclusive final states in the detector.
The MicroBooNE experiment was designed primarily to explore the nature of the MiniBooNE LEE observation using an 85-ton active volume LArTPC [32] through its excellent electron versus photon separation. The Micro-BooNE detector is located at the Booster Neutrino Beam (BNB) [35] at the Fermi National Accelerator Laboratory in Batavia, IL, USA. It is in the same beam line with a similar distance to the neutrino source as the MiniBooNE detector. While MicroBooNE is capable of searching for exotic signatures, such as electron-positron pairs [36], the MicroBooNE experiment's main goal is to search for an LEE signal in ν e CC interactions with an electron in the final state (referred to as the eLEE search) as well as in neutral-current (NC) interactions with a single photon in the final state. In these LEE searches, the MicroBooNE experiment has developed four analyses: i) a photon-based-search focusing on the ∆-decay hypothesis [37]; ii) an exclusive search using the ν e chargedcurrent quasi-elastic (CCQE) channel [38]; iii) a semiinclusive search for events with one electron and no pions in the final state [39]; and iv) a fully inclusive search for events with one electron and any final state, which is what this article describes. The results of the three eLEE searches are individually quantified with an empirical model constructed based on MiniBooNE data (see Sec. II) instead of performing a sterile neutrino oscillation model fit and summarized together in Ref. [40]. Looking ahead, a dedicated search of sterile neutrino(s) oscillations through comparison of multiple detectors at different baselines is the goal of the upcoming Short-Baseline Neutrino (SBN) program [41].
This paper is organized as follows. An overview of the search strategy is presented in Sec. II. The MicroBooNE experiment and the Wire-Cell event reconstruction pack-age is introduced in Sec. III. The event selections including ν µ CC and ν e CC selections 1 are described in Sec. IV. Section V describes the systematic uncertainties considered in this work. To search for new physics, validations of the overall model including the predictions of both signal and background are crucial. Various examinations of the model including i) signal prediction, ii) background prediction, and iii) modeling of neutrino energy reconstruction are summarized in Sec. VI. The physics sensitivities and the final result of the eLEE search are presented in Sec. VII followed by the summary in Sec. VIII.

II. OVERVIEW OF THE ANALYSIS
In this paper, we report the results of the eLEE search using the Wire-Cell event reconstruction package [42,43] after analyzing the MicroBooNE data set collected from February 2016 to July 2018. The search for an eLEE signal is equivalent to testing the null hypothesis, which corresponds to the nominal prediction without an eLEE signal. This prediction is formed based on the stateof-the-art understanding of the BNB flux [35,44], the neutrino-argon interaction cross section [45,46], the detector simulation [47][48][49], and data-driven constraints on ν e CC signal and background predictions. An alternative hypothesis is formed based on an eLEE model which represents an anomalous enhancement in the rate of intrinsic ν e CC events at true neutrino energies between 200 MeV and 800 MeV with a fixed spectral shape extracted from the MiniBooNE experiment. This alternative hypothesis has several advantages over a dedicated sterile neutrino oscillation model [50,51]. Firstly, this model is constructed from the MiniBooNE data with a minimum set of assumptions and is agnostic of sterile neutrino oscillation model parameters. Secondly, this model contains many fewer free parameters, which is beneficial for the single detector configuration of MicroBooNE. Thirdly, this model avoids complications introduced by simultaneous modifications in multiple channels including ν µ disappearance, ν µ to ν e appearance, and ν e disappearance in data. Finally, it allows for a quantitative comparison of results between MicroBooNE and MiniBooNE, particularly given the difficulty in explaining the shape of the MiniBooNE LEE in global fits with a sterile neutrino model [22].
This eLEE model used as the alternative hypothesis is obtained by unfolding the observed excess of electronlike events in MiniBooNE to true neutrino energy under a CCQE hypothesis and applying the excess-to-intrinsic ν e ratios directly to the flux of intrinsic ν e expected in the MicroBooNE detector. The MiniBooNE data and simulation results reported in 2018 [52] were used as an input to the D'Agostini's multi-dimensional unfolding 1 In this paper, νx refers to both νx and its counterpartνx. procedure [53]. We note that this eLEE model is constructed using the true neutrino energy and is limited below 800 MeV, beyond which the unfolded LEE signal is negligible. Furthermore, the additional information regarding the lepton kinematics reported by MiniBooNE in 2020 [21] is not taken into account. In addition, a 68.3% confidence interval of the eLEE strength parameter of 1 ± 0.21 (1 ± 0.08) considering both statistical and systematic (only statistical) uncertainties is estimated from the reported 4.8σ (12.2σ) LEE significance from the latest MiniBooNE 2020 result [21] considering both neutrino and anti-neutrino data taking. If only neutrino data taking is considered, the LEE significance would slightly reduce to 4.7σ (11.7σ) considering full (statistical) uncertainty, respectively. In the hypothesis test, we allow the normalization of this eLEE model to float and define a non-negative strength parameter x, such that x = 0 corresponds to the nominal prediction without an eLEE signal (eLEE x=0 hypothesis) and x = 1 corresponds to the prediction with an eLEE signal with magnitude equal to that of the unfolded median of the MiniBooNE LEE result (eLEE x=1 hypothesis). Figure 1 illustrates the prediction of the eLEE x=1 hypothesis with respect to that of the eLEE x=0 hypothesis for ν e CC interactions in the TPC active volume in MicroBooNE assuming 100% detection efficiency.
We fit this eLEE model to our data by minimizing a combined Neyman-Pearson (CNP) χ 2 test statistic [54] that incorporates the MicroBooNE's experimental uncertainties into a covariance matrix thereby obtaining a best fit value of x = x min . We compute the primary nested likelihood ratio test statistic, ∆χ 2 nested = χ 2 eLEEx=x0 − χ 2 eLEEx=xmin , for our data given a varying eLEE strength x = x 0 and obtain frequentist confidence intervals for the eLEE strength x following the Feldman-Cousins pro-cedure [55]. In addition to the ∆χ 2 nested , several other statistical tests are performed to provide additional information. They are i) goodness-of-fit (GoF) tests based on a Pearson χ 2 to examine the compatibility between data measurement and Monte-Carlo (MC) prediction, especially with available constraints using the conditional covariance method [56], and ii) a simple-vs-simple likelihood ratio test, ∆χ 2 simple = χ 2 eLEEx=1 − χ 2 eLEEx=0 with the CNP χ 2 to demonstrate the compatibility between data and the two signature hypotheses (eLEE x=1 and eLEE x=0 ).
Without a near detector to measure the neutrino flux and neutrino-argon interaction cross section, the search for an eLEE in the ν e CC channel would suffer from large systematic uncertainties. In order to compensate for the lack of a near detector and to maximize the physics sensitivity of this search, ν µ CC channel is used to constrain the ν e CCprediction. As to be elaborated in Sec. III A, the prediction of ν e and ν µ flux are strongly correlated, given the parent hadron species are mostly common. In addition, the ν e CC and ν µ CC interaction cross section is also strongly correlated, given the lepton universality. Therefore, the measurement of ν µ CC rates at different energies, which is proportional to the product of ν µ flux and ν µ -Ar cross section, is efficient in constraining the prediction of ν e CC rates at different energies, which is proportional to the production of ν e flux and ν e -Ar cross section. In practice, a seven-channel fit strategy is adopted. The seven channels are: Here, the FC events are defined to be events with the reconstructed TPC activity fully contained within the fiducial volume (3 cm inside the effective TPC boundary [57], which is the corrected boundary that takes a space charge effect [48,58] into account). All non-FC events are defined as PC events. FC events have many advantages such as higher signal-to-background ratio and better energy resolution when compared to the PC events. The primary channel that is sensitive to the eLEE search is FC ν e CC. FC ν µ CC and PC ν µ CC channels are used to provide constraints on the prediction of ν e CC interaction rate and its systematic uncertainties. PC ν e CC channel, which is less sensitive to the eLEE search, also provides constraints to some extent. The other three channels, FC CCπ 0 , PC CCπ 0 , and NCπ 0 , are used to constrain/improve the background prediction since π 0 events are one of the major backgrounds in the ν e CC selection. These seven channels are designed to be orthogonal to each other. The ν µ CC channel excludes ν e CC and CCπ 0 candidates, and the CCπ 0 channel excludes the ν e CC candidates.
To reduce bias, we adopted a blind analysis for this eLEE search. Two sidebands (near and far) are defined in addition to the signal region. The signal region is defined to contain the events with reconstructed neutrino energy E rec ν lower than 600 MeV and passing the ν e CC selection. The far sideband is defined to contain the events that fail a looser ν e CC selection or the events with E rec ν higher than 800 MeV. The far sideband essentially includes ν µ CC events as well as high-energy ν e CC events, which are not sensitive to the eLEE signal. The near sideband covers the remainder of the events between the far sideband and the signal region. A completely open data set with an exposure of 5.33×10 19 proton-on-target (POT) is available for development of the event reconstruction, selection, and statistical analysis. This open data set is less than 10% of the eventual data set of an exposure of 6.369 × 10 20 POT which is used in this analysis. The selection and analysis procedure were frozen after analyzing the open data set, and no changes were made after checking the consistency between the data measurement and MC prediction using the far sideband data set. The overall model including the central value predictions and their associated uncertainties is quantitatively validated with extensive goodness-of-fit tests. The power of these validations are further enhanced through a conditional covariance matrix formalism, which allows for a comparison between data and a constrained model prediction. Constrained by the sideband data, the systematic uncertainties of the model prediction can be significantly reduced, allowing for a in-depth examination of the systematic uncertainties and estimated correlations between data samples. The near sideband was openned subsequently for further examination and, only then, was the signal region unblinded. During the finalization of the analysis procedure, several fake data sets were generated based on different cross section models and injected with or without eLEE signals to test the robustness of the analysis procedure. All details of the fake data sets were initially blind to the analysers, and the analysis results of these data sets were as expected, either extracting or excluding an eLEE signal as appropriate given the truthlevel information released afterwards.

A. BNB Neutrino Flux
The booster neutrino beam line uses 8 GeV kinetic energy protons from the Booster accelerator to bombard a beryllium target. The primary hadrons produced in these interactions and the secondary hadrons further produced through additional interactions of primary hadrons and the surrounding materials are focused by a magnetic horn into a 50-m-long decay pipe. Decays of these hadrons produce a neutrino beam. The neutrino horn current polarity has been set to focus positive charged hadrons, resulting in a beam with a small antineutrino component for all data taking to date. The BNB flux is simulated based on a Geant4 frame-work [59] following the earlier work by the MiniBooNE collaboration [35]. The top panel of Fig. 2 shows the composition of the BNB flux in terms of different neutrino flavors seen by the MicroBooNE detector, which is located at 468.5 m on axis from the target. The dominant neutrino species is the ν µ , which is mostly produced by a π + two-body decay mode as shown in the middle panel in the same figure. For energies higher than ∼2.3 GeV, another two-body decay mode, K + → µ + + ν µ , becomes the main mechanism to generate ν µ . The ν e flux is about 0.5% of the overall neutrino flux. For neutrino energies below about 1.2 GeV, ν e 's are mostly produced by a three-body decay mode, µ + → e + +ν µ +ν e , with µ + originating from the π + two-body decay mode as shown in the bottom panel of Fig. 2. For neutrino energies higher than 1.2 GeV, ν e s are mostly produced by another three-body decay mode, K + → π 0 + e + + ν e . Since the parent hadron species are mostly common, the predictions of ν µ and ν e flux are strongly correlated, which supports the overall strategy of using ν µ CC events to provide constraints on the prediction of ν e CC events.

B. MicroBooNE Detector
The MicroBooNE detector [32] is a 10.4 m long, 2.6 m wide, and 2.3 m high LArTPC, located on-axis along the BNB at a distance of 468.5 m from the beryllium target, 72.5 m upstream of the MiniBooNE detector. It consists of approximately 85 metric tons of liquid argon in the active volume for ionization charge detection along with an array of 32 photomultiplier tubes (PMTs) [61] for the scintillation light detection 2 . Figure 3(a) shows the MicroBooNE TPC, which is housed in a foam-insulated evacuable cryostat vessel. The cathode-plane high voltage is set at -70 kV during normal operation, creating a drift field of 273 V/cm. The ionization electrons drift at a speed of 1.1 mm/µs [62] in the drift field. This corresponds to 2.3 ms drift time for the maximum 2.56 m drift distance.
As shown in Figure 3(b), there are three parallel wire readout planes at the anode side of the TPC. These planes are labeled as the "U", "V", and "W" planes 3 , and each plane contains 2368, 2368, and 3456 wires, respectively. The wire spacing within a plane is 3 mm, and the planes are spaced 3 mm apart. The wires in the W plane are aligned vertically and the wires in the U and V planes are oriented at ±60 • with respect to the vertical direction. The different orientations of the wires allow for determination of the positions of the ionization electrons within the plane that is transverse to the drift direction. Bias voltages for the U, V, and W planes are -110 V, 0 V, and 230 V, respectively, which satisfies the 2 One PMT stopped working after the run 1 data taking. 3 The "W" plane is sometimes also referred to as the "Y" plane. transparency condition that all drifting electrons pass through the U and V (induction) wire planes and are fully collected on the W (collection) plane. The induced current on each wire is amplified, shaped, and digitized through a custom designed front-end application-specific integrated circuit [63] operating at 89 K in the liquid argon. The direct implementation of readout electronics in the cold liquid significantly reduces electronics noise. The equivalent noise charge on each wire is generally below 400 electrons. A minimum ionizing particle usually produces in total 13,000 electrons at a single wire if the particle trajectory is perpendicular to the wire orientation [64]. Figure 3(b) also shows the light-collection system behind the anode wire planes. This light-collection system is used to detect scintillation light from the liquid argon providing the precise timing of particle activity, which is crucial to rejection of the high-rate cosmic-ray background necessary in a surface-operating LArTPC detector like MicroBooNE. Thirty-two 8-inch Hamamatsu R5912-02MOD PMTs [61] provide uniform coverage of the anode plane. An acrylic plate coated with tetraphenyl butadiene is installed in front of each PMT to shift the ultraviolet argon scintillation light to the visible part of the spectrum to which the PMT is sensitive. The magnitude of the detected light on each PMT provides position information for time-isolated particle activities, which is compared with the light pattern predicted from the ionization charge signals in the TPC. A successful match [43] of charge and light signals determines the association between individual TPC activity and light detection and, therefore, the time of the corresponding TPC activity.

C. Event Trigger and Readout
The BNB delivers proton pulses at a rate of 5 Hz and approximately 4×10 12 POT per pulse. Each pulse is called a beam spill and lasts ∼1.6 µs with 82 2-ns wide bunches of protons. The MicroBooNE detector is expected to have one neutrino interaction inside the TPC active volume per about 600 spills. Each beam spill initiates a hardware trigger in the MicroBooNE data acquisition (DAQ), which records 4.8 ms (digitized at 2 MHz) of TPC data of ionization charge signals and 6.4 ms (digitized at 64 MHz) PMT data of scintillation light signals. The PMT data contain two separate trigger streams, a forced trigger covering the beam spill which is referred to as the beam discriminator and a self-trigger which is referred to as the cosmic discriminator. In each event, there are on average 26 cosmic-ray muon induced activities observed in the 4.8 ms TPC readout window.
One such record of TPC and PMT data corresponds to one event which covers not only the beam spill but also the cosmic-ray muon activity in the proximity of the beam spill. The time window of an event accounts for the few-millisecond delay of the TPC readout signal coming  [60]. The maximum drift distance is 2.56 m with a drift electric field of 273 V/cm. The light-collection system, which consists of 32 PMTs, is located behind the three anode wire planes, which detect ionization charge.
from the relatively slow drift of ionization electrons. Following the hardware trigger, a software trigger is applied to reduce the data size. It requires distinct PMT signals within the beam spill window. This results in a reduction of the event rate by a factor of 22 but a negligible efficiency loss for neutrino interactions.

D. Monte Carlo Simulation
The simulated neutrino flux (introduced in Sec. III A) is provided to the event generator Genie v3 [45,46] to generate neutrino-argon interactions inside 4 and outside the cryostat. The latter case is also referred to as the dirt background. Genie v3.0.6, G18 10a 02 11a, was used, which includes improvements on the usage of the Valencia model [65][66][67] for the local Fermi gas nucleon momentum distributions, improvements in the CCQE and CC two-particles-two-holes (CC2p2h) interactions, and improvements in the treatments of final state interaction (FSI) and pion production, when compared to Genie v2 12 2 (used by previous MicroBooNE analyses [68][69][70][71]). In addition to the default configuration, the parameters governing the CCQE and CC2p2h models are adjusted according to the T2K CC0π cross-section results [72]. Given that T2K has a neutrino flux in a similar energy range to that of MicroBooNE, this additional adjustment represents an improved model prediction and uncertainty treatment [73], despite the fact that T2K data is on carbon-hydrogen and MicroBooNE's measurement is on argon. Corrections from carbon to argon targets are applied based on a smooth nuclear dependence (or A-dependence) that are determined in fits to inclusive and semi-inclusive electron scattering data.
The resulting final state particles of each simulated MC 5 event are processed using the LArSoft [74] software framework, which is a toolkit to perform simulation, reconstruction, and analysis of LArTPC events. The final state particles are propagated through the detector using the Geant4 toolkit [59] v4 10 3 03c. The resulting energy depositions are further ported to dedicated detector simulation programs taking into account all detector effects to simulate the ionization charge and scintillation light signals after taking into account the space charge effect. The space charge effect [48,58] is caused by the accumulation of positively charged ions inside the active volume. For on-surface LArTPC detectors such as MicroBooNE, cosmic-ray muons provide a constant source of positively charged ions, which distort the local drift electric field. Consequently, the ionization electrons bend toward the detector center when drifting to the anode plane and the reconstructed positions of ionization electrons appear to be closer to the detector center compared to their true position making the effective detector boundary smaller than the actual active TPC boundary. In addition to the deformation of reconstructed positions, the distortion in the electric field also changes the amount of ionization electrons and scintillation photons through the charge recombination process [75,76]. After comparing simulation (with these effects implemented) and calibration run data, the amount of ionization electrons is further scaled down in the simulation in order to include the effects from conversion of analog-to-digital converter values to the number of ionization electrons and positiondependent energy calibrations [49], which results in an average scaling factor of 0.86. The position and number of ionization electrons modified by recombination and the space charge effect is ported to the TPC detector simulation [47], which takes into account the charge transportation and diffusion [62]. The induced currents on the wires are simulated by convolving the ionization charge distribution at the wire plane with the position-dependent (at 1/10th of the wire pitch resolution) field response function, which are calculated by the dedicated Garfield simulation [77] following Ramo's theorem [78]. The induced current is further convolved with the electronics response before adding the inherent electronics noise from data and is then digitized to generate the final waveform on each wire channel.
The optical detector simulation models the light emitted by charged particles interacting with the detector and produces signals in photomultiplier tubes. A full optical simulation implements the Geant4 simulation of individual optical photons produced along the path of charged particles through both the scintillation and Cherenkov processes. These photons are stepped through several processes with the detector medium including Rayleigh scattering, reflection, and partial absorption in order to produce a realistic detector response to the light source. However, because of the vast number of pho-tons typically produced in a neutrino physics event, the full optical simulation can take hours or days per event. Therefore, a fast optical simulation was developed to overcome this problem for routine simulation tasks. This mode utilizes a library of stored data that represent the PMT acceptance of scintillation light signals 6 to sample an expected detector response given an isotropic emission of light at a certain point in the volume. The PMT response is further convolved with the time distribution of these photons to generate the digitized waveform.
Compared to early MicroBooNE analyses [68][69][70][71], the simulation used in this eLEE search adopts a scheme of overlaying the simulated neutrino interactions with dedicated beam-off data. These data are taken without the neutrino beam and are triggered by a random signal to mimic the neutrino beam gate. The simulated TPC and PMT waveforms from neutrino interactions are overlaid with the data waveform of a beam-off event. Such simulation is also referred to as overlay MC simulation. This scheme eliminates the systematic uncertainties in the simulation of excessive electronics noise and cosmicray muon activity. This scheme limits the statistics of the overlay MC sample because of the finite size of available beam-off data sample to be overlaid. Nevertheless, the statistics of the MC background sample is more than a factor of three larger than that of the data. As elaborated in Sec. V E, a Bayesian approach was developed to obtain an optimal estimation of MC statistical uncertainties, even when the predicted background is zero.

E. Event Reconstruction
An end-to-end automated event reconstruction chain containing various fundamental reconstruction techniques was developed and implemented in this analysis. TPC signal processing, which mitigates the excess noise [64] and deconvolves the detector response from the drift electric field and the electronics readout [47,79], provides the reconstructed ionization charge distributions for each wire to the subsequent calorimetry and topology reconstruction algorithms. A tomographic three-dimensional (3D) image reconstruction algorithm, Wire-Cell [42], is used as the core algorithm of the reconstruction chain. Wire-Cell uses reconstructed ionization charge at different times and readout wire 1D positions to reconstruct the 3D images of ionization electrons without topology heuristic assumptions (e.g. tracks from muons/pions/hadrons or EM showers from electrons/photons) prior to the pattern recognition stage. Other algorithms such as clustering and de-ghosting [43] are implemented to further improve the quality of the 3D images particularly addressing the challenge that gaps occur in the 2D view of the charge signals because of the inefficiency of TPC signal processing or non-functional wires. The space points (representing the 3D voxels with non-zero reconstructed charge) of the reconstructed 3D image are grouped into clusters, each of which represents individual physical activity in the TPC from a cosmic-ray muon or a neutrino interaction. Figure 4(a) shows each 2D projection view of an event's 3D image after imaging and clustering.
As previously discussed, the TPC ionization charge signal data, which provide the topology and calorimetry information, are collected separately from the PMT scintillation light signals, which provide the timing information, because of the longer drift time of ionization electrons relative to the light propagation. This results in a challenge, especially for surface-operating LArTPCs such as MicroBooNE, in identifying neutrino interactions from numerous cosmic-ray muon interactions [68,[80][81][82][83]. A many-to-many charge-light matching algorithm was developed to overcome this challenge by finding the corresponding light signals and providing the interaction time for each charge cluster [43]. About 70% of the cosmic-ray muon events that pass the software trigger are rejected by requiring a charge cluster to have its start time coincide with the beam spill. These clusters are referred to as in-beam clusters. Figure 4(b) shows an example from data, where the in-beam ν e CC candidate cluster is selected out of about 20 cosmic-ray muon clusters after charge-light matching.
The majority of in-beam clusters still originate from cosmic-ray muons after charge-light matching for Micro-BooNE events. Additional algorithms were developed to reject cosmic-ray backgrounds with 5-10% efficiency loss for neutrino interactions [84]. The effective boundary of the TPC active volume considering the space charge effect [48,58] is used to define the fiducial volume in the cosmic-ray rejection as well as the subsequent neutrino selection. The fiducial boundary is defined as 3 cm inside the effective boundary, which leads to a fiducial volume of 94.2% of the active TPC. The through-going muons, which traverse the TPC active volume, are rejected if the two ends of the track exit the fiducial boundary. The rejection of stopped muons, which enter the active volume and stop inside, is based on the identification of an increase of ionization charge loss per unit length (dQ/dx) near the end of the track (i.e., Bragg peak). This is obtained by a newly developed 3D track trajectory and dQ/dx fitting procedure [57], which is also an important ingredient in pattern recognition and particle identification. Note that the external cosmic-ray-tagger [85] system may provide additional rejection of cosmic-ray muon events, but is not included in this work as the system was not installed until late 2017. Using the above mentioned cosmic-ray rejection techniques including chargelight matching, the cosmic-ray background is reduced significantly resulting in less than 15% cosmic-ray contamination in the selected neutrino candidate events while the original neutrino to cosmic-ray muon ratio was about 1:200 for the events passing the software trigger. The efficiency loss for CC neutrino interactions is 10-20% for different flavors of neutrinos up to this stage [57,84]. Pattern recognition is vital for the identification of different flavors of neutrinos, e.g. ν e CC and ν µ CC, for a variety of physics analyses. The details of Wire-Cell pattern recognition can be found in Ref. [86] and highlights are provided in the following. Wire-Cell pattern recognition starts by finding initial end points of track segments by searching for kinks and splits in the selected 3D in-beam cluster. Track segments and their end points are then determined by iterative multi-track trajectory and dQ/dx fitting where linear algebra algorithms and graph theory operations are utilized to achieve a robust performance. Particle identification (PID) is performed based on the dQ/dx, topology information (direction, track or shower, etc.), and allowable particle flow relationships for each track segment. Candidate primary neutrino interaction vertices are concurrently identified as parts of the particle flow tree, which is a series of particles that starts from the neutrino interaction vertex and loops over all identified particles following the particle flow relationship. A deep neural network of SparseConvNet [87] is used to boost the performance of the neutrino vertex identification by predicting the distance from each 3D voxel to the neutrino vertex. It chooses from the neutrino vertex candidates, which are identified based on the above traditional algorithms, and determines the final reconstructed neutrino interaction vertex. The particle flow is then refined if needed. π 0 particles are reconstructed The blue, cyan, green, yellow, and red colors roughly correspond to 1/3, 1, 2, 3, and 4 times the dQ/dx of a minimum ionizing particle (MIP). (e) Reconstructed particle flow starting from the primary neutrino interaction vertex, which is displayed in a rainbow-colored wheel.
relying on the opening angle and topological information of the two decay γ's and other π 0 decay modes like Dalitz decay plays a very minor role in this analysis. A proton (p) can be reconstructed if its length is 1 cm which corresponds to a kinetic energy threshold of 35 MeV. Neutrons (n) are invisible in LArTPCs as they are neutral, but they could be suggested by off-vertex protons that result from the neutron scattering. Figure 5 illustrates the results of pattern recognition at different stages. As reported in Ref. [86], this pattern recognition achieves 60-75% efficiencies of good neutrino vertices (distance between reconstructed and true vertices below 1 cm) and subsequently a 80-90% reconstruction efficiency for primary leptons for charged-current neutrino interactions.

F. Neutrino Energy Reconstruction
Neutrino energy reconstruction for Wire-Cell inclusive neutrino selection adopts a calorimetric approach, essentially adding up the visible energy deposited in the TPC active volume. The reconstructed neutrino energy is given by the following formula: where i represents each identified particle in the reconstructed particle flow, K rec i represents the reconstructed kinetic energy, m i represents the rest mass value from the PDG [1], and B i represents the average binding energy (8.6 MeV) per nucleon after which an argon-40 nucleus is completely disassembled into its constituent nucleons. m i is only added for particles reconstructed as µ ± , π ± , and e ± particles, and B i is added for each reconstructed proton in the particle flow, either a primary proton connected to the neutrino vertex or a secondary proton scattered off from an argon nucleus by a neutron. There are three methods to calculate the kinetic energy, K rec i , for each particle: • Range: the range is used to calculate the kinetic energy of track-like particles µ ± , π ± , and p if they stop inside the active volume. This is based on the NIST PSTAR database [88] with a correction for different particle masses.
• Summation of dE/dx: for short (<4 cm) tracklike particles that stop inside the active volume or any exiting particles, the (visible) kinetic energy is estimated by summing up dE/dx for each piece (∼6 mm) of the track. The energy loss per unit length dE/dx is converted from the ionization charge per unit length dQ/dx considering the recombination effect when the ionization electrons are produced. This method is also used to estimate the energy of long muons with significant delta-rays emitted along the trajectory. This avoids bias in calculating the range. Note that an effective recombination model that also takes into account the overall normalization difference between Micro-BooNE data and simulation charge signals [49] is built by tuning the parameters of the modified box model from ArgoNeuT [76]. This effective recombination model (α=1.0 and β=0.255), as is used in current Wire-Cell energy reconstruction, has an improved agreement with MicroBooNE data; however, it still underestimates the energy by about 10%. Given that the ionization charge in simulation has been scaled to match that in the data, this bias from the effective recombination model ap- pears in both data and simulation in the same way. The residual difference between data and simulation is covered by the systematic uncertainties as discussed in Sec. V.
• Overall charge-energy scaling: for EM showers (non-track topologies), the energy is estimated by scaling the total reconstructed charge of the shower cluster by a factor of 2.50 after multiplying by the 23.6 eV per ionization pair [89,90]. This factor is derived from the nominal simulation and takes into account the bias in the reconstructed charge with respect to the true charge, and the average recombination effect [43] which converts the deposited energy to the true charge. For data events, an additional scaling factor of 0.95, which is calibrated from the reconstructed π 0 invariant mass (Sec. IV C), is applied.
Performance of this neutrino energy reconstruction is evaluated using MC samples by comparing the reconstructed values with the true values. The reconstructed neutrino energy versus the true neutrino energy as well as the energy resolution and bias relative to the true values are shown in Fig. 6. These figures present the results for selected ν µ CC and ν e CC candidates (see Sec. IV), respectively. The neutrino candidates are further divided into FC and PC samples (See Sec. II) based on the containment of the selected in-beam cluster. We should note that since an energetic muon track is much more extended than an electron shower with the same energy, the FC cut actually selects ν µ CC events with higher energy transfer (thus larger missing hadronic energy on average) than the corresponding ν e CC events at the same neutrino energy. For the FC samples, the reconstructed neutrino energy resolution is 15-20% for ν µ CC candidates with ∼10% bias (towards lower energies) and 10-15% for ν e CC candidates with ∼7% bias (towards lower energies). The slightly worse energy resolution of ν µ CC than that of ν e CC reflects the fact that ν µ CC event selection can tolerate more imperfect event reconstruction given its higher signal to background ratio than that of ν e CC. Energy deposited outside the active volume by charged particles or carried away by neutral particles such as neutrons is missing in this calorimetric energy reconstruction. In the energy range of interest, below 800 MeV in true neutrino energy, we have achieved uniform performance in the resolution of reconstructed neutrino energy and its bias, especially in ν e CC FC signal event selection ( Fig. 6(g)). More validation of this neutrino energy reconstruction in terms of the modeling of missing energy using ν µ CC events and data-MC comparisons can be found in Sec. VI B.

IV. EVENT SELECTION
The starting point of the neutrino event selection is the generic neutrino selection [57,84], in which the cosmicray backgrounds are reduced resulting in an overall contamination below 15%. After the generic neutrino selection, the efficiencies 7 for ν µ CC and ν e CC events are approximately 80% and 90% with signal-to-background ratios of about 2:1 (purity of 66%) and 1:250 (purity of 0.4%), respectively.
For the search of the low-energy excess in the ν e CC channel, the event selections are designed to be as general as possible (i.e., inclusive ν e CC), so that more freedom in examining exclusive channels would be available at later stages of the analysis if an excess was to be observed. Since the ν µ CC events are used to constrain the systematic uncertainties in neutrino flux, neutrino-argon interaction cross section, and detector effects, an inclusive ν µ CC selection is also adopted.

A. Charged-Current νe Selection
The development of ν e CC selection involves two stages. The first stage is the categorization of nonν e CC backgrounds, which is informed by hand scans of a small amount of background events. Then, variables with signal-background discrimination capability which represent the characteristic features of each type of background in the first stage are used as input into boosted decision trees (BDTs) trained on large MC simulation samples. Events used in BDT training are removed in making predictions.
The basic selection of inclusive ν e CC events requires an EM shower with a reconstructed energy higher than 60 MeV connected to the (primary) neutrino vertex [86]. The energy threshold is applied in order to exclude Michel electrons. When there are multiple reconstructed EM showers connected to the same neutrino vertex, the EM shower with the highest energy is taken as the primary electron candidate for further examinations.
The backgrounds are categorized into five major types. The first type focuses on primary electron identification, including the examination of the dQ/dx profile at the first few centimeters of the shower (i.e. shower stem) and the identification of a gap between the shower and the neutrino vertex. This gap occurs in photon showers due to the photon conversion length of approximately 18 cm in liquid argon. The second type of background focuses on interactions with multiple EM showers in the final state, most likely from π 0 production. The third type focuses on muon-related misidentification as electrons.
The fourth type focuses on more general background rejection using kinematic information, e.g. comparison of lepton kinematics between an electron candidate and its competing muon candidate in the same event 8 . The last type focuses on interactions with poor pattern recognition, which includes several different failure modes leading to incorrect pattern recognition. Beside these five major categorizations of background featuring in the ν e CC selection, there is a another set of dedicated taggers removing residual cosmic-ray muon induced backgrounds, which will be described in detail in the next section.
The primary electron identification includes: • Gap cut: the beginning of the EM shower in each 2D projection view is examined to search for a gap to the identified neutrino interaction vertex.
• Stem quality cut to remove backgrounds: the beginning of the shower is examined to ensure the quality of the shower stem. The checks include examinations of i) other minimum ionizing particle (MIP) tracks overlapping with the shower stem, and ii) possible track splitting, e.g. the splitting of pair-produced electron and positron instead of traveling in the same direction.
• MIP dQ/dx cut: we examine the dQ/dx profile of the shower stem to ensure a MIP (electron-like) event. We calculate the length of the MIP-like track below a MIP threshold cut (i.e. 1.3 times the expected MIP dQ/dx). The calculation of the length also considers the possibility of a delta ray (i.e. a single sample with high dQ/dx). In addition, a high dQ/dx value at the vertex because of additional vertex activities must be taken into account.
With the hand scanned features selected, we apply BDT techniques to high-statistics MC simulation samples to finalize the ν e CC selection. The usage of machine learning techniques mitigates the limitation of human learning when processing a large number of events. From among the many different machine learning tools, the BDT technique is chosen because it is more robust and approachable for general users compared to other multivariate analysis tools. The BDT package XGBoost [91], which provides fast and robust training through parallel tree boosting, is used. XGBoost also improves the model generalization and overcomes the issues of overfitting in gradient boosting enabling the use of a large pool of variables (over 300 variables from all categorizations of backgrounds in this ν e CC selection) in a single model. To train the BDT, the true ν e CC events in the fiducial volume that pass the generic neutrino selection and have at least one reconstructed electron EM shower are used to define the signal. In order to enhance the performance, the events with bad reconstruction, when the reconstructed neutrino and EM shower vertices are incorrectly reconstructed, are removed from the signal events. The overlay MC simulation after excluding true ν e CC events and the dedicated beam-off data are used as background in training the BDT. Figure 7(a) shows the signal efficiency and purity as a function of the ν e CC BDT score, and the distribution of ν e BDT score is shown in Fig. 7(b). Here, the efficiency is defined with respect to all true ν e CC events with their neutrino interaction vertices inside the fiducial volume. The cut value of 7.0 was chosen in order to maximize the ν e selection efficiency×purity value. A final ν e CC selection with 46% efficiency and 82% purity is achieved. The selected events are categorized into different types which are determined by the truth information. Each category, which will be used throughout this paper, as shown in Fig. 8(a) and Fig. 8(b), is defined as follows: (1) "EXT": cosmic-ray background from the beam-off data set that is external to the BNB data stream and triggered by cosmic-ray activity in coincidence with a fake beam spill, in which case the events have no BNB neutrino interactions; (2) "Cosmic": mistakenly selected cosmic-ray background from BNB overlay MC simulation, in which case each event has a simulated BNB neutrino interaction overlaid with dedicated beam-off data; (3) "Dirt": neutrino interactions with their true neutrino interaction vertices outside the cryostat, as defined in Sec. III D; (4) "out FV": neutrino-argon interactions with vertices outside the fiducial volume but within cryostat; (5) "ν µ CC π 0 in FV": ν µ CC interactions with vertices inside the fiducial volume and with at least one true π 0 in the final state; (6) "ν µ CC in FV": ν µ CC interactions with vertices inside the fiducial volume and with no π 0 in the final state; (7) "NC π 0 in FV": NC interactions with vertices inside the fiducial volume and with at least one true π 0 in the final state; (8) "NC in FV": NC interactions with vertices inside the fiducial volume and with no π 0 in the final state; (9) "ν e CC in FV": beam intrinsic ν e CC interactions with vertices in the fiducial volume; and (10) "LEE": all categories of excessive events originating from eLEE neutrino interactions, which accounts for the difference between the eLEE x=1 hypothesis and eLEE x=0 hypothesis. Except for "EXT", the other categories correspond to MC events which overlay simulated neutrino interactions with randomly triggered beam-off (cosmic) data. As shown in Fig. 8(c) and Fig. 8(d), the other categories for the selected neutrino interactions correspond to different interaction types obtained from the event generator. "CC" or "NC" represent charged-current or neutral-current. "QE", "RES", "MEC", and "DIS" represent quasi-elastic, resonance, meson exchange current, and deep inelastic scattering, respectively. Before applying the ν e CC selection to the BNB data stream, we validated its performance using the off-axis Neutrinos at the Main Injector (NuMI) [92] neutrino beam at FNAL. The NuMI beam is created from collisions of protons accelerated to an energy of 120 GeV with a graphite target. Similar to that of BNB, the charged hadrons are focused by a magnetic field into a 675 m long decay pipe. The distance between the NuMI target and the MicroBooNE detector is about 679 m. At an offaxis location of ∼ 8 • , the ν e 's with a sizable amount of ν e 's above 200 MeV are mostly coming from the 3-body K + and unfocused K 0 L decays. Compared to that of the BNB, the percentage of ν e in the flux is an order of magnitude higher at ∼5 %, which makes it ideal to validate the performance of the ν e CC selection. With 1.917×10 20 POT exposure, we select 269 FC ν e CC and 162 PC ν e CC candidates with reconstructed neutrino energy below 2.5 GeV and with an overall efficiency of 42% and purity of     91%. The overall ratio between data and nominal NuMI MC prediction for FC and PC ν e CC without considering any anomalous enhancement is 0.99±0.06 (stat.) and 1.08 ± 0.08 (stat.) indicating an overall good agreement.

B. Charged-Current νµ Selection
The generic neutrino selection results in a 88.4% selection efficiency and 65.0% purity for ν µ CC events [57,84]. Here, the selection efficiency is defined with respect to all true ν µ CC events with their neutrino interaction vertices inside the fiducial volume. The achieved purity is limited by the residual cosmic-ray muon background, neutrinoinduced background originating outside the fiducial volume, and NC interactions inside the fiducial volume. The precise reconstruction of the ν µ CC neutrino vertex (resolution is less than 1 cm) and particle identification (an integrated efficiency of 90% for primary muons) are leveraged to suppress these backgrounds. The reconstructed neutrino vertex is required to be inside the fiducial volume and the length of primary muons is required to be greater than 5 cm before applying the BDT selection.
In analogy to the ν e CC selection, human scans of the remaining backgrounds were performed to extract the main features of each type of background. The residual cosmic-ray background is typically the result of incorrect charge-light matching where the TPC cluster is placed at an incorrect location along the electric field direction. A through-going muon could have only one track end reconstructed at the detector boundary instead of two track ends reconstructed at the detector boundary, mimicking a single muon starting inside the TPC and exiting the detector. A stopped muon might also appear to be fully contained and, alternatively, be reconstructed with the candidate neutrino vertex at the muon decay vertex connecting the muon and the Michele electron. These topological features are leveraged to do this background rejection. For neutrino-induced background originating outside the fiducial volume, a charged hadron usually enters the detector and undergoes a hadronic interaction. For this kind of events, the neutrino vertex is typically reconstructed at the hadronic interaction point, and the event could then appear to originate inside the fiducial volume with a misidentified muon candidate. Note, with an exiting high-energy charged particle track, one may not achieve a reliable PID for MIP particles such as muons. The kinematics, especially the direction of the muon candidate, can be used to reject such background since most of the hadrons entering the detectors from outside of the detector are not as forward-going as expected.
For NC neutrino interaction background inside the fiducial volume, the main difference from ν µ CC events is the absence of a primary muon at the neutrino vertex. However, the separation of ν µ CC and this NC background, which mainly relies on the discrimination of charged pions and muons, is very difficult if only the dQ/dx information is used. To further reject such NC background, the activities associated with charged pions, e.g. proton scattering, and the relatively large-angle deflection (∼10 degrees) of the trajectory of charged pions can be used to provide additional separation power.
With the identification of the major features of the residual cosmic-ray background, neutrino-induced background originating outside the fiducial volume, and NC events inside the fiducial volume, the BDT was trained and applied to improve the ν µ CC selection. A similar training strategy as discussed in Sec. IV A is used with a signal definition switched to the ν µ CC events. Figure 9 shows the ν µ CC selection efficiency and purity as a function of the ν µ CC BDT score as well as the distribution of the ν µ CC BDT score. The final cut value of 0.9 was chosen for the ν µ CC selection with a 68% efficiency and 92% purity. Figure 10 shows the selected ν µ CC events and selection efficiency as a function of neutrino energy and muon cosθ. The efficiency is generally higher for more   forward-going angles as events with forward-going angles are more likely to have a typical topology of a ν µ CC event to which the BDT input variables are tuned. The "slope" of data-prediction ratios present in the bottom panel of Fig. 10(a) will be discussed in Sec. VI C. For some specific backgrounds, the BDT could select a small number of poorly reconstructed events. This is because the ν µ CC BDT was trained on a data set without an explicit request of a muon being present. Therefore, some ν µ CC candidates may not have a primary muon identified. This results in a non-zero number of entries in the first bin (less than 100 MeV) in Figs. 10(a) and 10(c) in which case the muon rest mass is not considered in the neutrino energy reconstruction and, also, in an absence of some events that have no reconstructed muon cosθ in Figs. 10(b) and 10(d). We can do this because of the high initial signal-to-background ratio of the events.

C. π 0 Selection
The ν µ CC selection described in the previous section is employed to select CCπ 0 events from ν µ CC interactions. Additionally, a NCπ 0 selection is constructed by considering the non-cosmic events that fail the ν µ CC selection (ν µ BDT score smaller than zero). In the reconstruction of π 0 events, the pair of EM showers with highest energies, supposedly from π 0 two-γ decay, is chosen to be the ones pointing to the same vertex. π 0 particles are identified as primary particles by placing a maximum distance cut between the neutrino vertex and the π 0 vertex. More details on the vertexing, clustering, and pattern recognition of π 0 events can be found in Ref. [86]. Further selection cuts use the γ energies, the distances between the neutrino vertex and the γ vertices, the opening angle between the two γ's, and the reconstructed π 0 invariant mass. The distribution of the reconstructed π 0 mass of the selected CC and NC π 0 events can be found in Fig. 11 where a Gaussian fit over the peak region of data returns a mean±width of 131.2±22.1 MeV/c 2 for the CCπ 0 selection and 130.4±19.3 MeV/c 2 for the NCπ 0 selection. The best-fit mass values have <1 MeV differences compared to the MC reconstructed mass values and are consistent with the expected π 0 invariant mass of 135 MeV/c 2 . The best-fit peak width does not include the contribution from the long tails, which is the result of imperfect event reconstruction. Table I lists the efficiency and purity of the π 0 selections. Here, the efficiency is defined with respect to all CC or NC π 0 events with their neutrino interaction vertex inside the fiducial volume. The dominant background component for the CCπ 0 selection is ν µ CC events without a π 0 in the final state, and the dominant background for the NCπ 0 selection is external events originating from cosmic-ray muons or neutrino interactions outside the fiducial volume.

V. SYSTEMATIC UNCERTAINTIES
In this analysis, we consider sources of systematic uncertainties from i) the neutrino flux of the BNB, ii) neutrino-argon cross sections of the Genie event generator, iii) hadron-argon interactions of the Geant4 simulation, iv) detector response resulting from imperfect calibration, and v) the finite statistics of MC samples used for prediction, as well as vi) additional uncertainty for dirt events, which originate from neutrino interactions outside the cryostat. The various sources of systematic uncertainty each have different impacts on the reconstruction or selection efficiency (for both signal and background) as well as the reconstruction of kinematic variables of the predicted events. The uncertainty due to limited statistics of MC samples and beam-off data is particularly important for estimating backgrounds of rare event searches, e.g. the selection of BNB ν e CC, which is ∼0.5% of the total flux. Figure 12 summarizes the relative uncertainties from all systematic uncertainty sources for the seven channels. In the following, we describe each uncertainty in detail.

A. Uncertainties from the Model of Neutrino Beam Flux
The calculation of the BNB neutrino flux and its uncertainties follows earlier work from the MiniBooNE collaboration, which includes a well constrained beamline simulation based on the Geant4 framework [59] as well as techniques to handle systematic uncertainties [35]. Our flux prediction uses the updated flux calculation that takes into account the SciBooNE measurement of p + Be → K + production in the BNB [93,94], which provided a better constraint on kaons produced in the BNB. In addition, our flux prediction evaluates π + and π − production uncertainties directly from HARP pion production data [95] rather than using a fit parameterization. This technique allows the HARP measurement uncertainties to be more properly propagated to the calculated neutrino flux.
As summarized in Table II, the systematic uncertainties in the predicted flux include effects from i) hadron production of π + , π − , K + , K − , and K 0 L and ii) nonhadron production: modeling of the horn current distribution, horn current calibration, and pion and nucleon

Tuning parameter name
Parameter type π + hadron production FLUX π − hadron production FLUX K + hadron production FLUX K − hadron production FLUX K 0 L hadron production FLUX horn current distribution FLUX horn current calibration FLUX nucleon total scattering Xs FLUX nucleon inelastic scattering Xs FLUX nucleon quasi-elastic scattering Xs FLUX pion total scattering Xs FLUX pion inelastic scattering Xs FLUX pion quasi-elastic scattering Xs FLUX MicroBooNE Genie All Genie Xs (µB tune) Strength of the CCQE RPA correction Genie Xs (µB tune) Parameterization of the CCQE nucleon axial form factor Genie Xs Parameterization of the CCQE nucleon vector form factors Genie Xs Changes angular distribution of nucleon cluster in MEC Genie Xs (µB tune) CCMEC Cross-section Shape Genie Xs (µB tune) Angular distribution for RES ∆ → N + π Genie Xs Angular distribution for RES ∆ → N + γ Genie Xs (µB tune) Scaling factor for CC coherent π production Genie Xs (µB tune) Scaling factor for NC coherent π production Genie Xs (µB tune) Second-class vector current Xs Second-class axial current Xs π − interactions Geant4 π + interactions Geant4 proton interactions Geant4 For each tuning parameter, a covariance matrix is generated according to its number of universes. The final covariance matrix is the summation of all individual covariance matrices. For Genie Xs, the label "µB tune" indicates that the tuning parameter and/or its uncertainty is from MicroBooNE Genie tune other than the default treatment from Genie v3.0.6. total, inelastic, and quasielastic scattering cross-sections on beryllium and aluminum. There is also an overall 2% normalization uncertainty associated with the POT counting. Figure 13 shows the fractional contributions to the overall flux systematic uncertainty from each source, for the seven channels of the analysis. In the low-energy region, e.g. E rec ν < 500 MeV, which is relevant to the search for a ν e low-energy excess, the flux systematic is limited by the hadron production of π + , which, with its decay to muons, produces most of the ν e 's and ν µ 's in this energy range. Figure 14 shows the correlations of flux systematics for the seven channels. There are obviously strong correlations between the low energy ranges of ν e events and ν µ events given that they all originate from π + decays. There are strong correlations in the high energy range between ν e and ν µ given that these largely originate from π + and K + decays. There are also strong correlations between the high energy ν µ and the entire energy range of the π 0 events since the neutrino interactions that enter into the π 0 channels are almost all from ν µ interactions (negligible ν e ) which pass the energy threshold of π 0 resonance interactions.

B. Uncertainties from the Neutrino-Argon Interaction Cross Sections
As introduced in Sec. III D, Genie v3.0.6 [46], with parameters governing the CCQE and CC2p2h models adjusted [73] by MicroBooNE, is used to generate exclusive neutrino-argon interactions. Table II summarizes the uncertainties considered in the reweighting procedure for neutrino-argon interaction cross sections (Xs). In particular, the contribution labeled "MicroBooNE Genie All" contains 46 tuning parameters that are simultaneously varied in generating hundreds of universes. The variations in these universes are used to construct the covariance matrix. These tuning parameters cover a wide range of models including CCQE, CC resonance (CCRES), CC non-resonance, CC transition, CC deep-inelastic scattering (CCDIS), NC interactions, and final-state interactions. In addition to "MicroBooNE Genie All", as listed in Table II, there are additional 9 "Genie Xs" parameters taking into account other cross-section uncertainties in Genie, each of which have only two universes corresponding to the nominal value and the 1σ bound of the tuning parameter, respectively. The 1σ uncertainty is taken to be the difference. There are two other tuning parameters labeled "Xs" focusing on second-class currents [96] that may contribute to ν e /ν µ CCQE cross section differences as suggested by Ref. [97]. Figure 15 shows the correlations of cross section systematics for the seven channels. There are generally strong correlations between ν e CC and ν µ CC across the entire energy range, which is a natural result of the lepton universality assumption. In the low-energy region, e.g. E rec ν < 500 MeV which is relevant to the search for a ν e low-energy excess, the neutrino cross section systematic is limited by the level of suppression of the CCQE cross section behavior at low Q 2 (four-momentum transfer from the lepton) because of long-range nucleonnucleon correlations, which is poorly constrained by existing data. We should further note that some inconsistencies were identified in the Genie v3.0.6 reweighing code used to evaluate FSI-related systematics as also discussed in Ref. [73], but the effect of these were found to have a negligible impact on the overall analysis sensitivity and have been ignored.

C. Uncertainties from the Hadron-Argon Interaction
Charged hadrons can scatter, both elastically and inelastically, with external argon nuclei via hadronic interactions. These interactions can lead to the production of additional particles or can cause large angle changes in particle trajectories that may affect the reconstructed neutrino energy. Therefore, it is important to include the uncertainties of these secondary hadron-argon interactions. Studies were performed separately for elastic and inelastic interactions for protons and charged pions.
Geant4 is used to propagate all hadrons through the detector medium based on a semi-classical cascade model [98]. Using the Geant4reweight [99] package, events with inelastic hadronic interactions are reweighted independently for interactions containing protons, positive pions, and negative pions. For each particle type and independent of the particle's energy, the inelastic cross section is varied around its mean by O(20%), based on the uncertainties of world data. Figure 16 shows the correlations of hadron-argon interaction cross section systematics for the seven channels. The contribution (square of the uncertainty) from this source of systematic uncertainty to the total systematic uncertainty is up to 1% as indicated in Fig. 12.

D. Uncertainties from the Model of Detector Effects
There are four major categories of detector systematic uncertainties each targeting one area of the detector calibration effort: i) variations related to the light yield (LY) and propagation simulation, ii) variation in the space charge effect, iii) variation in the recombination model, and iv) variations on the amplitude and width of the deconvolved ionization charge waveforms for other discrepancies between the detector response used in simulation and that in data. We vary the magnitude of the space charge effect based upon the measurements of the spatial distortions at the edge of the TPC, extrapolated to the bulk field [48,58]. A different recombination model, which provides slightly better agreement to the data, is used to estimate the data/MC difference in the dE/dx to dQ/dx conversion. The variations in the wire waveform include effects as a function of x (drift direction), (y, z) (vertical/beam directions), θ xz and θ yz (the angular orientation of the particle's trajectory with respect to the global coordinate system); these variations are constructed by comparing the waveforms in data and simulation. More details on the variations of the wire waveform can be found in Ref. [100].
For each source of detector systematic uncertainty, the same set of MC simulation events are re-simulated with a change to the detector modeling parameter of interest. The systematic uncertainty is estimated by comparing efficiency and reconstructed kinematic variables in the new and old simulation. The change of each detector modeling parameter is treated as 1σ, meaning that there is, in principle, only one degree of freedom in constructing each detector covariance matrix after factoring in the statistical uncertainties. The usage of the same set of events in the old and new detector simulation aims to reduce statistical fluctuation, which is estimated using the bootstrapping method [101].
The central idea of bootstrapping is to re-sample existing MC events in order to estimate the uncertainties on the central value of the prediction. The basic procedure is illustrated in the following: • Choose a common set of events in the central value (CV) and 1σ variation simulation samples.
• Re-sample MC events to form distributions with statistical uncertainty corresponding to that of the expected POT exposure. During this process, one can naturally take into account the event weight in the re-sampling process.
• Apply the event selection requirements to the resampled events in both the CV and the 1σ samples and calculate the difference in the spectra of different channels. Each such vector is considered one "universe." By repeating this procedure many times, many universes are simulated to form the covariance matrix, M R . This covariance matrix essentially represents the uncertainty on the difference between the CV and 1σ sample (nominal difference vector, V nominal D ). In other words, these are uncertainties of the uncertainties.
• In order to generate the detector covariance matrix (M D ), a 2-step procedure is adopted. In the first step, based on M R , generate a random set of vectors (based on decomposition of the symmetric matrix) δ V D to be applied to the nominal difference vector ( V nominal D ) in order to obtain a new vector • In the second step, generate a single random number r using the normal distribution. Then multiply V D by this random number (bin-to-bin fully correlated) to obtain r · V D , which is treated as one universe.

MicroBooNE Simulation
FIG. 17: Correlations of total detector systematics for the seven channels. The correlations between the FC and PC ν µ CC channels are dominated by the components of light yield and recombination.
• By multiple repetitions of the above two steps, many universes are simulated, which can the be used to construct the detector covariance matrix.
In the aforementioned approach, the statistical uncertainties in the detector variation samples are naturally taken into account in the M R matrix. Correlation of the detector variation samples is ensured by using the single random number r. The relative strengths of the correlated and uncorrelated uncertainties are properly dealt with via this procedure. Figure 17 shows correlations of all detector systematics for the seven channels, including statistical uncertainties. The correlations seen between the FC and PC ν µ CC channels and between the FC ν µ CC and ν e CC channels mostly come from the variations in light yield and recombination. Figure 18 shows the fractional contributions to the overall detector systematic uncertainty from each type of variation for the seven channels of the analysis.

E. Uncertainties of Finite Statistics in Making Predictions
For a rare event search with low background predictions such as that in the MicroBooNE eLEE analysis, statistical uncertainties of the MC simulation and dedicated background data can be important. For example, when the predicted background is zero, quoting a zero uncertainty on this prediction following the general error propagation would lead to an underestimation of statistical uncertainties. On the other hand, if we use an upper limit to quote the statistical uncertainty, the size of statistical uncertainties would increase rapidly through the general error propagation procedure when there are multiple components in estimating the background. In order to solve this problem and obtain an optimal estimation of MC statistical uncertainties, we adopt a new approach by combining the Bayesian approach with the covariance matrix. For the Poisson distribution assuming a unit step function as the prior (1 for µ ≥ 0 and 0 for µ < 0), an accurate approximation of the posterior distribution of the expectation µ given the observation of N is given by where the high-order asymptotic expansions (based on Stirling's approximation) of the Gamma function, Γ(µ + 1), are ignored. For example, in the case of N = 0, we have P (µ|0) = e −µ . If we use the covariance matrix formalism to estimate the RMS of µ given a central value of zero, we would have an uncertainty of 1.41 (standard deviation relative to the central value). Similarly, we can estimate the uncertainties for other values of N . In MC simulation, different events may be associated with a different weight value, which can lead to deviations from the simple Poisson distribution. For example, if a MC event is associated with a weight value of two, it would present two events in making a prediction. In the literature, there are many studies of the treatment of a likelihood function given this situation. For a recent review, see Ref. [102]. In this section, we give a prescription for how to deal with this situation. The actual implementation follows Ref. [103]. Given events (labeled by i) with different weights (w i ), we define the mean and variance of the prediction as: which can be written using effective uniform weights as Define so that m Eff = m 2 /σ 2 approximately follows a Poisson distribution with an expectation value of m Eff , and then the likelihood function (or distribution) can be written as where µ is the expected number of events in data with non-uniform weights. Accordingly, Eq. 2 can be modified by replacing "µ" and "N " with m Eff and m Eff . The above estimation can easily be extended to the case where several estimations are added together (e.g., a situation with multiple backgrounds). In practice, we first divide all distributions into two categories: one with non-zero event count predictions and one with predictions of exactly zero event count. These two distributions are then convolved. Finally, we correct for the Bayesian prior. For implementation, we give prior distributions when we convolve multiple (n) distributions together. If the prior information of each distribution is flat, the overall prior of the summation would be max(µ,0) n−1 , or µ n−1 for µ ≥ 0 and 0 for µ < 0. In this case, it is corrected (divided) after the convolution of multiple distributions.

F. Additional Systematics for Dirt Events
The dirt events are neutrino interactions originating outside the cryostat. The biggest uncertainty associated with the dirt events is the modeling of the geometry and materials of the foam insulation covering the cryostat outer surfaces, the nozzle penetrations for cryogenic and electronic services, the supporting structures, and the dirt around the detector facility. In addition to the systematic uncertainties associated with neutrino flux, cross section, and detector, we conservatively assign an additional 50% bin-to-bin uncorrelated uncertainty to dirt events.  Figure 19 shows correlations of all sources of systematic uncertainty across the seven channels. Figure 20 summarizes the fractional contribution of each source of uncertainty to the total systematic uncertainty, for each of the seven channels. In general, without any constraint, the cross section systematic uncertainty dominates in the region where neutrino interaction rate is relatively high, followed by the flux systematic uncertainty. With constraint, as discussed in later sections, cross section and flux uncertainties will be significantly reduced while the data statistical uncertainty becomes predominant in this analysis. For ν e CC channels, the relative contribution of MC statistical uncertainty, which includes statistical uncertainty from both beam-off data and MC simulation in making predictions, is significantly higher than that of the ν µ CC and π 0 channels mainly because of the limited statistics of the samples in estimating the backgrounds.

MicroBooNE Simulation
The equivalent exposure in POT of the beam-off data sample to estimate the beam-off (EXT) background is only 2.5 times of the exposure of the BNB data sample, and the POT of the MC samples to estimate the beam-on background is only 3.7 times of the exposure of the BNB data sample. The situation will be improved with more MC production or data taking. For the ν µ CC and π 0 channels, the detector systematic uncertainties become the dominating uncertainty at high energies where the number of events in the distribution become smaller (See. Fig. 21). This effect is expected with the bootstrapping method and is the result of the limited MC statistics of the detector systematic CV and 1σ detector variation samples in these regions.

VI. MODEL VALIDATIONS
This section summarizes the validations of the overall model, which includes i) the prediction of ν e CC signal through the examination of data and model consistency in ν µ CC events assuming lepton universality with the difference in lepton masses included in the neutrinoargon interaction cross section model, ii) constraints on the ν e CC backgrounds through the examination of data and model consistency in the various π 0 channels, and iii) modeling in the neutrino energy reconstruction.

A. Methodology
The covariance matrix formalism is adopted to construct the χ 2 test statistic: where M and P are vectors of measurement and prediction, respectively. The Cov (M, P ) is the full covariance matrix: (8) Cov full = Cov stat CNP + Cov sys MC stat + Cov sys xs + Cov sys flux + Cov sys det + Cov sys add . The Cov stat CNP is the diagonal covariance matrix constructed based on the combined-Neyman-Pearson (CNP) method [54] with the statistical uncertainty square being 3/(1/M i + 2/P i ) for the ith bin. The Cov sys MC stat is the diagonal covariance matrix containing the statistical uncertainties corresponding to finite statistics in making predictions as described in Sec. V E. The other four covariance matrices Cov sys xs , Cov sys flux , Cov sys det , Cov sys add , are the covariance matrices corresponding to uncertainties from cross section (Sec. V B), neutrino flux (Sec. V A), detector performance (Sec. V D), and dirt (Sec. V F), respectively. The dependence on the LEE strength x is considered in calculating the covariance. It is based on dedicated systematic and statistical covariance matrices for the eLEE x=1 component and correlations between the eLEE x=1 component and the other channels. The goodness-of-fit (GoF) test can be performed to test the compatibility between the data and the overall model. Following the recommendation of Ref. [104], we adopt the Pearson χ 2 construction (instead of the CNP construction) for the data statistical uncertainty: with the statistical uncertainty square being P i for the ith bin. Given the null hypothesis (i.e. eLEEx=0), the χ 2 value can be used to perform the GoF test by comparing to the χ 2 distribution with the associated number of degrees of freedom (ndf ), which is the total number of bins of the measurement. The above GoF test using a Pearson χ 2 construction provides an overall evaluation of the model and the null hypothesis compared to the data. This evaluation can be used to study different parts of the model using the conditional covariance matrix formalism [56]. For example, given the full covariance (stat + sys) containing two channels (X, Y): we can derive the conditional mean and conditional covariance of the prediction of X given the constraints from the measurement of Y : Thus, a goodness-of-fit test can be performed on Y first, and then performed on X after the constraints from Y .
This allows the examination of the model on X and Y individually, which provides more information about the compatibility between the model and data. The current cross-section model we use has conservative uncertainties in general, and along with correlated systematics, the reduced χ 2 values, which is the ratio between χ 2 and number of degrees of freedom, are generally low suggesting that the model describes the data well within the given uncertainty.

B. Validation of Monte Carlo Models with
Goodness-of-fit Tests As previously introduced, we use a seven-channel fit strategy in searching for a low-energy excess. The seven channels are FC ν e CC, PC ν e CC, FC ν µ CC, PC ν µ CC, FC CCπ 0 , PC CCπ 0 , and NCπ 0 . The selections for these seven channels are exclusive from each other. In particular, the candidates in the CCπ 0 excludes the selected ν e CC candidates, and the candidates in the ν µ CC channel excludes ν e CC and CCπ 0 candidates. Figure 21 shows the five non-ν e CC channels in the eLEE search that serve as constraints to reduce the systematic effects in the measurement of ν e CC. The χ 2 /ndf value in each plot is calculated based on the data-MC difference and the full covariance matrix taking into account the bin-tobin correlations from the off-diagonal terms. The error bands in the plot correspond to only diagonal terms in the covariance matrix. Figure 21 shows that the data and MC simulation agree within uncertainties, and the χ 2 /ndf values without any constraint suggest consistency albeit with conservative uncertainties. To further validate the MC mod-   III: χ 2 /ndf without and with constraints, for various selection channels. The central values and uncertainties (covariance) of the MC prediction are both changed with the constraint. Pearson construction of χ 2 , instead of CNP construction of χ 2 , is used for these goodness-of-fit tests. els, i.e. data-MC consistency, GoF tests using the conditional mean prediction and conditional covariance of the prediction, as described in Sec. VI A, are performed on these channels before applying them in the final eLEE search. After constraints, the χ 2 /ndf values based on both of the reduced systematic uncertainty and adjusted central values of the MC prediction are summarized in Table III. The χ 2 /ndf values after constraints are generally increased to a certain degree mostly because of the reduced systematic uncertainty. For the PC ν µ CC, the uncertainty is reduced by about 75% with constraints from the FC ν µ CC. For the three π 0 channels, the uncertainties are reduced by about 30% with constraints from the FC and PC ν µ CC channels. All the χ 2 /ndf values af- ter constraints are still less than 1 indicating consistency between data and MC within the estimated systematic uncertainties of the MC prediction.

C. Validation of energy reconstruction
One may note that the data-to-prediction ratios demonstrate a slope for ν µ CC channels, which motivates further energy reconstruction validation, particularly on inclusive hadronic final states, in addition to the validations on dE/dx and EM shower energy scale using dE/dx profiles from various reconstructed particles and π 0 invariant mass.
As discussed in Appendix, this "slope" can be decomposed using different final state topologies, such as via 0pXπ/NpXπ (N≥1 and X≥0) separation, where p represents reconstructed protons with kinetic energy greater than 35 MeV, and there is no requirement regarding whether charged pions are reconstructed. The enhancement in the data measurement is found to be concen-trated in the 0pXπ channels which essentially correspond to the events with low hadronic energy while NpXπ channels have an overall good agreement between data and MC. To quantitatively assess this issue, a conditional constraint study was performed with the high-statistics ν µ CC events in this analysis following the methodology described in Sec. VI A. The more precise measurements of leptonic kinematics, such as muon energy (E µ ) and muon polar angle relative to the BNB direction (cosθ µ ), are utilized to constrain/adjust the MC prediction on the hadronic kinematics, such as the reconstructed hadronic energy (reconstructed E ν -E µ ), which is expected to contain significant missing energy from invisible neutral particles.
The measurements of muon kinematics and hadronic energy constraint results are shown in Fig. 22. The measurements of E µ and cosθ µ show good data-MC consistency for the MC predictions with or without certain constraints. The measurement of hadronic energy shows an excess for low hadronic energy events. With constraints from the measurements of E µ and cosθ µ , the measure-ment of hadronic energy, as in Fig. 22(c), shows a good agreement (χ 2 /ndf =8.66/32) with the constrained MC prediction even though the common systematic effect between leptonic and hadronic systems (e.g. neutrino flux uncertainties) is significantly reduced resulting in a much smaller residual systematic uncertainty (blue error band). This is an important evidence that the origin of this issue is essentially from the common systematics, most likely cross-section modeling, shared by both leptonic and hadronic final states and can be adequately described/allowed by the current models and their variations in MicroBooNE MC simulation. More discussion can be found in Ref. [105]. This constraint study is sensitive to the potential missing energy in the hadronic system. For example, 5-10% (relative to the true energy transfer) additional invisible energy would result in a noticeable increase of the χ 2 /ndf values in this test. See the supplementary material of Ref. [105] for more details. This stringent test validates the MC model as well as our neutrino energy reconstruction treatment (leptonic and hadronic energy) and gives us confidence using the overall MC model in the seven-channel fit.

VII. RESULTS
In this section, eLEE sensitivity, nue CC selected data and MC comparisons (Sec. VII A), and various statistical analyses (Sec. VII B and Sec. VII C) including the bestfit of the eLEE strength and Feldman-Cousins confidence intervals are reported. The eLEE model used here is generated by unfolding the MiniBooNE's observed excess to true neutrino energy under a CCQE hypothesis and scaling the flux of intrinsic ν e with the excess-to-intrinsic ν e ratios.
As explained in Sec. II, a blind analysis strategy was implemented for the MicroBooNE LEE analysis which sequesters ν e CC candidates in BNB data below a reconstructed neutrino energy of 600 MeV until the analysis had been finalized. This procedure ensured that there was no bias in the reconstruction or event selection while allowing for cross-checks of data selection and simulation performance using sideband events with similar kinematic and/or topological characteristics to those of ν e CC signal.
A. νe CC Selection Results ν e CC events are identified as described in Sec. IV A starting from a generic neutrino selection that is followed by various dedicated selection criteria based on Wire-Cell event reconstruction and BDTs. The ν e CC events are selected with a high purity as shown in Fig. 8. To further maximize the eLEE sensitivity, we divide the ν e CC candidates into two sub-samples: a fully contained (FC) and a partially contained (PC) sample, both of which follow the aforementioned seven-channel fit strategy. The FC and PC ν e CC selection results are shown in Fig. 23 and Fig. 24. The green histogram represents the MC predicted ν e CC signal expectation obtained from the beam simulation and the Genie interaction models tuned to MicroBooNE (eLEE x=0 hypothesis). The figures also contain the data-to-prediction ratios, χ 2 per number of degrees-of-freedom, and error bands calculated based on the simulation of the expected ν e CC signal and background. The MC model of the ν e CC signal with the eLEE x=1 component added to the beam expectation is indicated by the red dashed histogram. The χ 2 /ndf values (assuming eLEE x=0 hypothesis) for FC and PC ν e CC channels without and with the constraints from the ν µ CC and π 0 channels are summarized in Table III.
The neutrino energy spectrum of the data compared to the constrained MC prediction for the FC ν e CC channel, of which the low energy region is most sensitive to the eLEE search is shown in Fig. 25. The applied constraints include the PC ν e CC channel 9 and various ν µ and π 0 channels. The constrained MC prediction with the eLEE component added is indicated by the red dashed histogram in the figure. The number of predicted events for each category of background, beam intrinsic ν e CC, and eLEE ν e CC in the energy region of reconstructed E ν < 600 MeV are summarized in Table IV. Definitions of each category can be found in Sec. IV B. The main background comes from NC π 0 events in which one of the π 0 decay γ's (the other one may exit the active volume) is mis-identified as a primary electron EM shower. The ν e CC selection is also contaminated by CC or NC interactions with no π 0 in the final state in which case the typical failure mode is that the decay electrons from muons or charged pions are mis-identified as primary electron EM showers alongside bad reconstruction of the neutrino vertex. Note that the "Cosmic" category corresponds to the cosmic-ray background estimated from the overlay MC. This background is essentially a mismatch between charge and light signals which accounts for the mis-selected cosmic-ray background from neutrino activity in addition to the cosmic-ray background estimated from beam-off data with no neutrino activity ("EXT" category). Since this "Cosmic" background is associated with neutrino interactions in the overlay MC, the predicted number will change accordingly in the constraint.
In Table. IV, expected number of events for each type of background, beam intrinsic ν e CC, and eLEE x=1 ν e CC in the energy region of reconstructed E ν <600 MeV for FC ν e CC are shown. The uncertainties for each category of the predicted background or ν e CC signal are systematic uncertainties originating from the flux, cross section (including both ν-Argon and hadron-Argon interactions), detector, and the limited statistics of the samples used in making predictions (e.g. MC statistical uncertainty).     "Dirt" background has an additional 50% systematic uncertainty. The statistical uncertainties are presented only for beam-on and beam-off data results. With the constraints, the systematic uncertainty is reduced significantly for the ν e CC signal. However, the uncertainties of the predictions of various backgrounds change only slightly because the uncertainty of each background is dominated by the MC statistical uncertainty. With the constraints, the contributions from flux and cross section systematic uncertainties are largely reduced and the central values of the predictions are adjusted as well. The constraining power is mostly from ν µ CC channels, while π 0 channels help constrain the residual π 0 background in ν e CC candidate events. The PC ν e CC channel has an insignificant effect in the constraint or in the eLEE search, and it serves as a good validation check on the data and simulation. After such constraints, the uncertainty of the result is dominated by statistical uncertainty followed by the systematic uncertainties from the MC statistics, detector response, cross section, and neutrino flux. The data comparisons to the energy spectra of the MC predictions before and after constraint, without and with eLEE x=1 ν e CC, for the FC ν e CC channel are shown in Fig. 26. No significant discrepancy between data and the MC eLEE x=0 hypothesis (obtained from the BNB beam simulation with the MicroBooNE tuned Genie interaction models) was observed after applying all available constraints. The two data points between 500 MeV and 700 MeV show deficits compared to the constrained prediction of the eLEE x=0 hypothesis. A quantitative examination considering full uncertainties yields a local p-value of 0.039 (2.1σ), which is consistent with the hy- pothesis of statistical fluctuation and/or systematic variations. The goodness-of-fit χ 2 /ndf (Pearson format) and p-values for the entire energy region as well as the low energy region (less than 600 MeV) are summarized in Table V. Comparing the χ 2 /ndf values, one can see the impact from the constraint on the sensitivity to the eLEE search as well as the level at which the current data measurement excludes the eLEE x=1 hypothesis. More de-   V: Summary of the goodness-of-fit χ 2 /ndf values of FC ν e CC energy distributions without or with the constraint from the other six channels for eLEE x=0 and eLEE x=1 hypotheses, respectively. In all model goodness-of-fit tests, Pearson construction of χ 2 , instead of the CNP construction of χ 2 , is used. The pvalue for each goodness-of-fit test is also shown.

B. Simple vs. Simple Likelihood Ratio Test
Using the χ 2 (CNP format) as defined in Eq. (7), we construct a simple-vs-simple test statistic: where x represents the expected eLEE strength in the prediction. This test statistic allows one to calculate the p-value using a frequentist approach (i.e. with pseudo experiments) assuming the MC prediction with the eLEE x=0 hypothesis or with the eLEE x=1 hypothesis is true. In generating the pseudo experiments, the events in all channels are randomized according to their associated systematic uncertainties (covariance matrix) and statistical uncertainties (Poisson distribution). In the case the expected counts in a given bin is negative, which can happen for large systematic uncertainties, the corresponding pseudo experiment is discarded and a new psuedo experiment is generated. This test statistic tests the compatibility between the data from all seven channels and the given hypothesis. Note that as opposed to Table V which focuses on the low energy FC ν e CC events, both the FC ν e CC channel and the PC ν e CC channel that have eLEE sensitivity and the other ν µ CC and π 0 channels are simultaneously employed in the hypothesis test here and in the following section, in which case the ν µ CC and π 0 channels effectively serve as constraints.  Figure 27 shows the ∆χ 2 simple values and distributions of data and pseudo experiments for eLEE x=0 and eLEE x=1 hypotheses separately. Assuming the eLEE x=0 hypothesis is true, the p-value of the data measurement is derived to be 0.326 (one-sided), which corresponds to a 0.45σ significance level. This result demonstrates good agreement between the data measurement and the eLEE x=0 hypothesis. Assuming the eLEE x=1 hypothesis is true, the p-value of the data measurement is derived to be 9.0 × 10 −5 (one-sided), which disfavors the eLEE x=1 hypothesis at a 3.75σ significance level.

C. Nested Likelihood Ratio Test and Best-fit of eLEE strength
Besides the simple-vs-simple test statistic, another test statistic is constructed to further test the compatibility between observation in all seven channels and various hypotheses and to find the best fit of an eLEE model where the strength of the eLEE signal, x, is allowed to vary such that x = 1 would be the eLEE model strength based on the extrapolation from MiniBooNE and values less than or greater than that indicate a smaller or larger excess compared to MiniBooNE. The test statistic is defined as follows: The value of x 0 represents the null hypothesis: x 0 = 0 represents the eLEE x=0 hypothesis (expectation from the beam simulation and the tuned Genie interaction models), and x 0 = 1 represents the eLEE x=1 hypothesis. x min is the best-fit value of x in the allowed region (x min ≥ 0) after minimization. This test statistic is a nested likelihood ratio test statistic comparing the null hypothesis (x = x 0 ) with an alternative hypothesis corresponding to x min , which can in principle take any value within the allowed region, and is included to quantify the search for an eLEE along with the Feldman-Cousins approach [55].
The best-fit value of the eLEE strength x min is determined by minimizing χ 2 eLEEx with the covariance matrix as defined in Eq. (8) that accounts for the statistical and systematic uncertainties varying with different eLEE strengths. As shown in Fig. 28, the bestfit eLEE strength is at the boundary x min = 0 which corresponds to the eLEE x=0 hypothesis. The Feldman-Cousins 68.3%, 95.5% and 99.7% confidence level upper limits from the data measurement are calculated to be 0.217, 0.513, and 0.911, with the expected upper limits for the eLEE x=0 hypothesis at 0.243, 0.563, and 0.958, respectively.
In comparison, the 68% confidence interval (CI) of the eLEE strength estimated from the MiniBooNE 2020 result [21] with full and statistical-only uncertainties are shown in Fig. 28 as well. The lower limits of the 68% MiniBooNE full and statistical-only CIs are disfavored at significance levels of more than 2.6σ and 3.0σ, respectively. Note that some of the systematic uncertainties are correlated between the MiniBooNE and MicroBooNE experiments and a direct comparison of the two experiments taking into account correlated systematic uncertainties is unavailable. Therefore two separate significance levels derived based on the published information are used to indicate the bounds of this comparison.
Based on this nested likelihood ratio test statistic, the significance level of rejecting the eLEE x=1 hypothesis is estimated following the Feldman-Cousins procedure. As indicated by the black dashed line in Fig. 29  nested assuming eLEE x=1 is true. The sensitivity of rejecting the eLEE x=1 hypothesis is at 2.93σ following the Feldman-Cousins procedure. Using the same procedure, the data measurement rejects eLEE x=1 hypothesis at the p-value of 0.002, which corresponds to 3.14σ. The results are obtained from 10 million pseudo experiments.
assuming the eLEE x=1 hypothesis is true, the eLEE x=1 hypothesis is disfavored at a p-value of 0.002, which corresponds to a 3.14σ significance level. In comparison with the data result, the sensitivity of rejecting the eLEE x=1 hypothesis is also calculated using the Asimov data set of the eLEE x=0 hypothesis (best-fit of eLEE strength, x min = 0) following the procedure in Ref. [106]. The constraint on FC and PC ν e CC channels from the other five channels are taken into account in the Asimov data set. The resulting sensitivity, as indicated by the red dashed line in Fig. 29, is that the eLEE x=1 hypothesis, which represents the median of the MiniBooNE result, is disfavored at 2.93σ with the ∆χ 2 nested value of 11.26. Various hypothesis test results as well as sensitivity values (using the Asimov data set) are summarized in Table VI. This nested likelihood ratio hypothesis test obtains a similar level of rejection of the eLEE x=1 hypothesis as the aforementioned simple-vs-simple likelihood ratio hypothesis test, and provides a rigorous cross-check of the result.   Fig. 29. The significance level is calculated based on the pseudo experiments corresponding to each null hypothesis. The sensitivity is calculated using an Asimov data set assuming the eLEE x=0 hypothesis is true.

VIII. SUMMARY
In this article, we report a search for a low-energy excess in ν e CC interactions using BNB data in the Micro-BooNE experiment. A data-driven LEE model motivated by the previous observation of an electron-like low-energy excess (eLEE) from the MiniBooNE neutrino experiment is built to quantify the search result. With the single Mi-croBooNE detector, measurements of ν µ CC interactions and CC or NC interactions with a π 0 in the final state are used to constrain the prediction of ν e CC interactions as well as other neutrino backgrounds reducing the systematic uncertainty and thus maximizing the eLEE search sensitivity.
A blind analysis scheme was adopted in analyzing the data with a BNB exposure of 6.369×10 20 POT, which corresponds to the first three years of MicroBooNE data taking. After unblinding the ν e CC events with reconstructed neutrino energy below 600 MeV and using a nested likelihood ratio test statistic using all available information, the best-fit LEE strength is determined to be 0 (eLEE x=0 ) with the Feldman-Cousins 68.3%, 95.5%, and 99.7% C.L. upper limits at 0.217, 0.513, and 0.911, respectively. Using a simple-vs-simple hypothesis test, the eLEE x=1 hypothesis is rejected at 3.75σ, while the eLEE x=0 hypothesis is shown to be consistent with the observation at 0.45σ. Various hypothesis test results can be found in Table VI. Regarding the ν e CC hypothesis to explain the electronlike low-energy excess observed in the MiniBooNE ex-periment, we compare the confidence intervals obtained from the MicroBooNE and MiniBooNE experiments in the context of the eLEE model. The MiniBooNE 1σ statistical uncertainty band is entirely outside the 3σ allowed range of the LEE strength parameter x derived from this analysis in MicroBooNE. Even in the absence of any consideration of cross-experiment correlations, the MiniBooNE 1σ full uncertainty band is well outside the 2σ allowed x parameter range from MicroBooNE. We should note that the current eLEE model only takes into account the excess as a function of true neutrino energy with other kinematics modeled as for the BNB intrinsic ν e CC events. If we consider the lepton kinematic distributions reported by MiniBooNE [21] with an enhanced excess in the forward direction, the achieved ν e CC selection efficiency in this work, which is better at forward angles, would yield a stronger exclusion of the eLEE signal.
While we observe a null result in searching for a ν e LEE signal using the current available data set, further tests of alternative hypotheses explaining the Mini-BooNE low-energy excess using single-photon-like events due to processes either within or beyond the standard model [107][108][109][110][111] will be carried out and reported in future MicroBooNE publications. While this analysis did not perform a fit using a sterile-neutrino oscillation model, future analyses in MicroBooNE with a full data set at 12.25×10 20 POT exposure, as well as the short-baseline neutrino program [41], will explicitly test short-baseline oscillation models in both appearance and disappearance channels. Additional support for the laser calibration system and cosmic ray tagger was provided by the Albert Einstein Center for Fundamental Physics, Bern, Switzerland. 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.
APPENDIX: Separation of 0pXπ and NpXπ samples As shown in Fig. 22 and discussed in Sec. VI C, the enhancement of the ν µ CC events with low hadronic energies can be adequately described by the simulation, in particular by the cross-section model and its allowed variation which is most relevant to the potential missing energy in the hadronic system, employed in the Micro-BooNE MC simulation.
To further validate the reconstruction of interactions with inclusive hadronic final states and cross-check its impact on the analysis results, the inclusive ν µ CC events ( Fig. 21(a) and 21(b)) are further divided into 0pXπ and NpXπ (N≥1) channels where p is defined as reconstructed protons with kinetic energy greater than 35 MeV (corresponding to a length of 1 cm) and X is the number of reconstructed pions (X≥0). The separated 0pXπ and NpXπ channels for the selected ν µ CC events are shown in Fig. 30. Interestingly, the enhancement of ν µ CC events is found to show up in the 0pXπ channel, whereas the NpXπ channel has an overall good agreement between the data and MC prediction. This suggests that the enhancement of the inclusive ν µ CC events in the relatively low energy region is less likely to be caused by potential bias in the neutrino flux prediction. Based on the nominal MC simulation, the ν µ CC 0pXπ FC+PC events, excluding true ν e CC and NC events, have 41% (46%) true 0pXπ events, and the NpXπ FC+PC events have 96% (97%) true NpXπ events in the full energy region (low energy region of E rec ν < 600 MeV). The kinetic energy threshold of the protons in the truth-level study is 40 MeV. To complement the above truth-level studies which indicate contamination of NpXπ events in the 0pXπ ν µ CC selections, we assess from hand-scans of both data and MC simulation that protons which contaminate the selection are generally low energy. It supports the statement that the ν µ CC channels successfully split interactions by the amount of proton hadronic activity, consistent with the goal of this study. This indicates the excess 0pXπ events are not due to the migration of NpXπ events with high energy protons. Finally, the hypothesis of migration of high energy events into the low energy region is disfavored because the deficit in the high energy region for both 0pXπ and NpXπ ν µ CC channels has insufficient events to account for the excess in the low energy region.
The inclusive ν e CC events (Fig. 24) are divided into 0pXπ and NpXπ channels as well in order to investigate the consistency between ν µ CC and ν e CC selections assuming the enhancement of 0pXπ events is largely attributed to the cross section modeling. The separated 0pXπ and NpXπ channels for the selected ν e CC events are shown in Fig. 31, with no constraints applied in the MC predictions. Based on a truth-level study similar to that for ν µ CC events, the ν e CC 0pXπ FC+PC events, excluding true ν µ CC and NC events, have 39% (49%) true 0pXπ events, and the NpXπ FC+PC events have 94% (98%) true NpXπ events in the full energy region (low energy region of E rec ν < 600 MeV). The data-   prediction ratio in the relatively low energy region of the ν e CC 0pXπ channel appears higher than that in the ν e CC NpXπ channel, in analogy to the ν µ CC results. In order to understand the impact of, and account for the observed behavior in the data-prediction ratios in the 0pXπ and NpXπ channels on the eLEE analysis, the eLEE strength fit is repeated separating the inclusive ν µ CC and ν e CC channels by these exclusive final-state topologies. The default 7 channels used in the eLEE analysis are expanded to 11 channels where each ν µ CC or ν e CC channel is separated into two (0pXπ and NpXπ). The results can be found in Fig. 32. The 11-channel result combining separated 0pXπ and NpXπ channels is well consistent with the default 7-channel result. Furthermore, the individual 0pXπ or NpXπ results also have best-fit eLEE strengths at zero albeit the uncertainties and sensitivities are different. Therefore, the consistency between ν µ CC and ν e CC selections with respect to the different hadronic final states is validated, and it supports the inclusive 7-channel fit strategy as well as the overall conclusion of this paper. Additional work studying ν µ CC interaction cross-sections in exclusive final states is ongoing within the collaboration with the goal of further improving our ability to model these processes. FIG. 32: ∆χ 2 nested value as a function of eLEE strength. "7-ch" (black curve) corresponds to the result from the default 7-channel simultaneous fit which is the same as the one in Fig. 28. "11-ch" (green curve) corresponds to the result from a 11-channel simultaneous fit combining separated 0pXπ and NpXπ channels for ν e CC and ν µ CC FC and PC events together with the dedicated π 0 channels. "0pXπ" (red curve) corresponds to the result from the 0pXπ channels for ν e CC and ν µ CC FC and PC events as well as the dedicated π 0 channels. "NpXπ" (blue curve) corresponds to the result from the NpXπ channels for ν e CC and ν µ CC FC and PC events as well as the dedicated π 0 channels.