When a photon hits the bi-alkali photocathode surface of an R5912-02mod PMT,
an electron is released via the photoelectric effect [162]–this electron is
often referred to as a “photoelectron” (p.e.). Each p.e. is focused toward a
dynode chain–a series of electrodes designed to produce a number $n>1$
electrons for each incident electron, resulting in an avalanche of electrons
by the final anode which can be read out in the form of a current. The
wavelength-dependent probability with which a given p.e. enters the dynode
chain is known as the “quantum efficiency“ of the PMT. At the nominal
operating temperature of 87 K, the MicroBooNE PMTs have an average quantum
efficiency of 15.3% [16]. Each PMT was tested in a liquid nitrogen (77 K)
cryogenic environment, in which the gain and rate of thermal emission (“dark
current”) of the PMT were measured as a function of the supplied high voltage
(HV) across the dynode chain [163]. The HV for each PMT was set to produce a
gain of $3\times 10^{7}$ at 77 K [16], corresponding to a dark current of
$\mathcal{O}$(kHz) [163]. The current output from each PMT passed through
preamp/shaper boards before being digitized via an analog-to-digital converter
(ADC) with a sampling rate of 64 MHz [16]. Thus, light is collected in time
ticks with a length of 15.625 ns [148].

The light collection system was critical in detecting activity in the detector
coincident with a beam spill, indicating the presence of a neutrino
interaction. In order to record a given event, at least 5 photoelectrons must
have been detected across all PMTs [148]. Additionally, a “common optical
filter” was applied to reduce the non-neutrino trigger rate–this filter
required at least one string of six time ticks with greater than 20
photoelectrons detected during the 1.6 $\mu$s beam spill window, and no such
strings of six time ticks in the 2 $\mu$s prior to the beam spill. The
spatiotemporal distribution of light observed by the PMTs can also be used to
augment the signal from the wire planes. For example, the PMT signal from a
cosmic muon which stops in the detector and decays to a Michel electron is
shown in figure 3.5(b).

The scintillation light yield and detection efficiency in liquid argon are
impacted by a number of phenomena. Impurities can dissociate the argon
excimers before they have a chance to emit scintillation light. For example,
O2 molecules in the LAr volume can undergo a two-body collision with an argon
excimer [164],

${\rm Ar}_{2}^{*}+{\rm O}_{2}\to 2{\rm Ar}+{\rm O}_{2}.$ (3.1)

This interaction mainly decreases the decay probability of the longer-lived
triplet state. Because of this, measurements of the delayed scintillation
lifetime in LAr are sensitive to the concentration of impurities in the
detector [160, 164, 16].

Scintillation light can also be absorbed by impurities within the detector. At
concentrations of around 2 ppm, dissolved nitrogen will decrease the
attenuation length of 128 nm scintillation light to around 30 m [165, 158].
TPB emanation from the painted acrylic plates can also lead to a bulk
fluorescence effect within the liquid argon [166]. Additionally, Rayleigh
scattering can deflect scintillation photons, diminishing the detector’s
capability to translate photon detection into spatial information regarding
the path of the original charged particle. The Rayleigh scattering length for
128 nm photons in liquid argon has been measured to be approximately 55 cm
[167]. Given the size of the MicroBooNE detector, scintillation photons will
undergo around five Rayleigh scattering interactions on average before
reaching a PMT.

Figure 3.4: The LAr scintillation light spectrum, TPB ultra-violet absorption
spectrum, TPB emission spectrum, PMT quantum efficiency, and PMT surface
transmission efficiency as a function of photon wavelength. Figure from Ref.
[17].

(a)

(b)

Figure 3.5: Figure 3.5(a) shows a photograph of a single PMT system used in
the MicroBooNE detector. The acrylic window here has not yet been coated with
TPB. Figure 3.5(b) shows the PMT signal from a stopped cosmic muon (top) that
decays to a Michel (bottom). Figures from Ref. [16].

###  3.2 TPC Signal Processing

MicroBooNE has developed a custom cryogenic electronics system to amplify and
digitize the analog signal on each wire of the LArTPC [16]. The wire signals
first pass through front-end ASICs located on cryogenic motherboards attached
directly to the wire carrier boards on the LArTPC. Each ASIC amplifies and
shapes the signal from 16 wires–thus, 516 ASICs are needed to fully instrument
the 8,256 readout channels of the LArTPC. Cryogenic-tested cables carry
signals from the ASICs to intermediate amplifiers directly outside the
cryostat, after which the signals travel through 20 m cables to the readout
system. This system consists of 130 readout modules, each containing 8 custom-
designed 12-bit ADCs and a Front End Module (FEM). Each ADC digitizes the
signal from 8 wires at a rate of 16 MHz. The dynamic range of the ADC is set
to enable the detection both of a minimum-ionizing particle far away from the
anode planes and of a highly-ionizing-particle close to the anode planes
without saturation. The FEM uses a field-programmable gate array (FPGA) to
store and compress the digitized signals. Once a trigger is received, the FPGA
sends 4.8 ms of TPC data to the data acquisition (DAQ) system. The 4.8 ms
window is digitized into 9600 time ticks on each wire, with each time tick
corresponding to 500 ns of integrated charge [148].

####  3.2.1 Noise Filtering

The MicroBooNE LArTPC underwent an initial engineering run from October 2015
to July 2016. Various sources of noise in the TPC were identified using data
from this run, all of which are described in detail in Ref. [18]. The
frequency dependence of the noise was studied to inform the creation of an
offline noise filter.

Some of this noise is inherent to the electronics system, coming from the
cold, front-end ASIC or the two RC circuits in the intermediate amplifier and
ADC board. The gain and peaking time of the front-end ASIC are configurable
from among four pre-defined settings each, which determine the dynamic range
and timing granularity of the ASIC. This choice also fixes the level of
irreducible noise in the ASIC, comprised of multiple sources with different
frequency dependence [18]. The timing response of the RC circuits can distort
large amplitude, long duration signals such as a cosmic muon traveling
parallel to a single wire. This effect can be corrected via a deconvolution
process [18].

The engineering run also revealed a number of excess noise sources above the
inherent electronics noise. The most prominent of these was low frequency
($\lesssim 30$ kHz) noise from the low voltage regulators on the front-end
ASICs. Harmonic noise was also observed corresponding to the ripple frequency
of the HV power supply for the field cage. Offline filters were developed to
subtract the noise induced from both of these sources. As shown in figure 3.6,
these filters are able to remove most of the noise in the wire planes while
having minimal impact on the true signal from ionization electrons. After the
engineering run, hardware upgrades were made to the low voltage regulators and
the HV power supply which greatly reduced the excess noise in the TPC readout
signals. Even after these upgrades, the offline noise filters were still
applied to optimize the TPC noise level.

Figure 3.6: 2D displays of the signal from wires in one of the induction
planes in a single data event, before and after the application of the offline
noise filters. Figure from Ref. [18].

####  3.2.2 Deconvolution

Even after noise filtering, the signal from the TPC is a complex combination
of a number of factors: the ionization electron distribution, the induced
current in each wire due to passing ionization electrons, and the response of
the electronics system. The signal must undergo deconvolution in order to
obtain the quantity of physical interest–the number of ionization electrons
passing each wire as a function of time. The signal processing algorithm
designed to extract this information is described in detail in Ref. [19] and
validated using MicroBooNE data in Ref. [168]. We give a brief overview of the
algorithm here.

The measured signal $M(t)$ on a given wire is a convolution of the original
signal $S(t^{\prime})$ and the detector response function $R(t,t^{\prime})$,
which is intended to capture all the factors impacting signal formation
discussed above. This can be expressed by the integral

$M(t)=\int_{-\infty}^{\infty}R(t,t^{\prime})S(t^{\prime})dt^{\prime}.$ (3.2)

For a time-invariant detector response function
$R(t,t^{\prime})=R(t-t^{\prime})$, one can express this relationship in the
frequency domain via a Fourier transform. In a realistic setup, however, this
treatment tends to amplify the effect of high-frequency components of the
inherent electronics noise. Typically one introduces a filter function
$F(\omega)$ to suppress high frequency response, such that the original signal
can be determined from

$S(\omega)=\frac{M(\omega)}{R(\omega)}F(\omega).$ (3.3)

In this setup, one can consider $F(\omega)$ as a replacement for the detector
response function.

The filter function used in MicroBooNE is inspired by the Weiner filter [169],
which is designed to optimize the signal-to-noise ratio while minimizing
variance and bias and the resulting $S(\omega)$. The exact Weiner filter is
not applicable in the realistic MicroBooNE setup, as it does not conserve the
total number of ionization electrons and can lead to non-local charge smearing
due to suppression of low frequencies [19]. MicroBooNE instead constructs a
Weiner filter $F(\omega)$ from simulation, and then fits the result to a
functional form designed to preserve the normalization of the measured signal
and remove the unwanted low-frequency suppression behavior.

The procedure discussed so far only accounts for the signal generated in a
given wire due to ionization electrons passing by that wire. In reality,
electrons passing by neighboring wires can also induce a current in a given
wire. This is addressed by using a two-dimensional (2D) deconvolution in both
the time and wire dimensions rather than the one-dimensional deconvolution in
the time dimension introduced above. An additional Weiner-inspired filter is
constructed along the wire dimension when performing 2D deconvolution. After
the 2D deconvolution, the algorithm identifies a region of interest (ROI)
around each potential signal. These ROIs act as high-pass filters which
suppress low-frequency noise in the deconvolved signal while preventing any
non-local charge smearing.

The full deconvolution algorithm results in a robust reconstruction of the
amount of charge passing by or deposited on each wire of the LArTPC. Figure
3.7 shows the U plane event display for a neutrino interaction candidate in
MicroBooNE data after the application of the deconvolution algorithm. After
the 2D deconvolution, one can clearly identify all of the charged particle
tracks coming from the neutrino interaction vertex. The shape of the
deconvolved signal for wires in each plane matches closely between data and
simulation, as discussed in Ref. [168]. The signals produced by the
deconvolution algorithm are ready to be used by MicroBooNE physics analyses,
including the search for the MiniBooNE electron-like excess presented in
chapters 4 and 5.

Figure 3.7: A U plane event display of a candidate neutrino interaction in the
MicroBooNE data. The impact of the 1D and 2D deconvolution algorithms on the
post-noise-filtering signal is shown. Figure from Ref. [19].

##  Chapter 4 The MicroBooNE Electron Neutrino Analysis: Overview and
Selection

The MicroBooNE experiment was conceived to investigate the nature of the
excess of low-energy electron-like events observed by the MiniBooNE
experiment, described in detail in chapter 2. The MicroBooNE LArTPC, described
in chapter 3, produces high-resolution event displays of neutrino interactions
which enable the separation of events with electrons in the final state from
those with photons. The collaboration released its first results searching for
a MiniBooNE-like excess of electromagnetic events in October 2021. Four
different analyses were performed, one focusing on photons and three focusing
on electrons. The photon analysis searched for an excess of $\Delta\to
N\gamma$ radiative decays [80]–the only SM photon background not constrained
in situ by MiniBooNE. The three electron analyses searched for an excess of
$\nu_{e}$ charged-current (CC) interactions, each with a different final state
signal definition [116]. One of the electron analyses, hereafter referred to
as the “inclusive analysis”, isolated all $\nu_{e}$ CC interactions using the
WireCell reconstruction framework [170]. Another electron analysis, hereafter
referred to as the “MiniBooNE-like analysis”, used the Pandora reconstruction
framework to isolate all $\nu_{e}$ CC interactions without pions in the final
state [171]. This chapter describes the “two-body CCQE analysis”, which used
deep learning techniques to isolate charged-current quasi-elastic (CCQE)
$\nu_{e}$ interactions consistent with two-body scattering [148]. The results
from this analysis are presented in chapter 5. The two-body CCQE analysis
represents much of my early work as a graduate student.

Publications covered in this chapter for which I either held a leading role or
made major contributions: [148, 116, 172]

###  4.1 Dataset and Simulation

The results covered in this section come from the first $6.88\times 10^{20}$
POT of BNB data taken by the MicroBooNE experiment. These data were separated
into three run periods which spanned from 2016 to 2018: run 1 ($1.75\times
10^{20}$ POT), run 2 ($2.70\times 10^{20}$ POT), and run 3 ($2.43\times
10^{20}$ POT) [148]. Data from two more MicroBooNE run periods exist on tape
and are currently being processed by the collaboration, corresponding to
another roughly $5\times 10^{20}$ POT. The entire dataset was taken during
nominal neutrino mode running conditions for the BNB, in which 8 GeV kinetic
energy protons from the Fermilab Booster strike the beryllium target to
produce a cascade of mesons, and the positively-charged mesons are focused and
decay to produce a beam of neutrinos. The BNB was described in detail in
section 2.1.1. The intrinsic $\nu_{e}$ component of the beam, which makes up
$\sim 0.5$% of the total neutrino flux, can be seen in figure 2.2.

The MicroBooNE simulation relies on the GENIE v3.00.06 [173, 174, 175, 176]
neutrino event generator. MicroBooNE specifically employs the G18_10a_02_11a
model, which uses the Valencia CCQE and meson-exchange-current (MEC) models
[177] as well as a local Fermi gas model to describe the nuclear environment.
GENIE models the final state interactions (FSI) of particles produced in
neutrino interactions with the nuclear environment using the typical $A^{2/3}$
scaling, including corrections from nuclear binding energy [178] as well as
Pauli blocking and the velocity dependence of the nuclear mean field [179]. To
improve the nuclear modeling of CCQE and MEC interactions, MicroBooNE has
developed a custom tune of the G18_10a_02_11a model using a $\nu_{\mu}$
CC$0\pi$ cross section measurement from T2K [180, 181]. This so-called
“MicroBooNE tune” was determined by performing a two-dimensional fit of the
T2K double-differential cross section data in muon momentum $p_{\mu}$ and
angle $\cos\theta_{\mu}$. The fit adjusts four parameters related to the CCQE
and MEC cross section models in GENIE v3.00.06, resulting in a better fit to
both T2K and MiniBooNE double-differential cross section data (though
discrepancies exist at lower muon momenta in the latter case). The MicroBooNE
tune is used for all three MicroBooNE $\nu_{e}$ analyses; it gives better
agreement with MicroBooNE $\nu_{\mu}$ data, which are essential for
constraining the $\nu_{e}$ prediction in each analysis, and updates the
predicted $\nu_{e}$ CCQE and MEC cross sections on argon [180].

Electron scattering data can be used to evaluate the vector part of the
vector-axial neutrino-nucleus cross section models in event generators such as
GENIE v3.00.06 [182, 183]. Earlier versions of GENIE were found to be in
disagreement with electron scattering data in the sub-GeV energy transfer
regime [182]. Updates to the QE, MEC, and $\Delta$ excitation models in GENIE
improved agreement with the QE peak and MEC-dominated dip region, though
disagreements remain at higher momentum transfers [183]. Dedicated electron
scattering measurements using the CLAS detector at Jefferson Lab have also
revealed biases in the GENIE-based reconstruction of the incident electron
energy, especially in the large final state transverse momentum regime [184].

Particle propagation in the liquid argon detector is handled by the Geant4
toolkit V10.3.03c [185], and the MicroBooNE TPC response is modeled using the
general LArSoft [186] package and the internal uboonecode package. The effect
of cosmic particles crossing the TPC during the readout window is modeled by
overlaying off-beam data on top of the simulated events. This is important, as
around 20-30 cosmic rays can cross the detector in a single event due to the
long readout window and near-surface location of the MicroBooNE detector
[148].

###  4.2 The Electron Low Energy Excess Template

The hypothesis studied in each of the three electron analyses is that the
MiniBooNE anomaly comes from an excess of $\nu_{e}$ interactions from the BNB.
While such an excess could come from $\nu_{\mu}\to\nu_{e}$ oscillations in a
$3+1$ model, this is not the only possible source. For example, mismodeling of
$\mu^{+}$ decay-in-flight and $K^{+}$ decay-at-rest in the BNB [14] could lead
to an enhancement of low energy $\nu_{e}$ compared with the expectation.

MicroBooNE has elected to remain agnostic about possible sources of $\nu_{e}$
events for its first results. To do this, the collaboration has developed an
“electron low energy excess (eLEE) model” by unfolding the MiniBooNE excess
under a $\nu_{e}$ hypothesis [148, 20]. This is accomplished using the
D’Agostini iterative unfolding procedure [187], which updates a predicted
distribution in true space using an observed distribution in reconstructed
space and a known detector response matrix connecting the two. The algorithm
relies on Bayes theorem to update the predicted distribution over a specified
number of iterations. A small number of iterations will give an unfolded
distribution biased toward the initial guess, while a larger number of
iterations will give an unfolded distribution that is more sensitive to
statistical fluctuations in the observed data [20]. For the eLEE model, the
MiniBooNE observation as a function of reconstructed neutrino energy, as shown
in figure 2.5, is used to unfold the predicted $\nu_{e}$ interaction rate in
MiniBooNE as a function of the true neutrino energy. This procedure accounts
for the detection efficiency; thus, the unfolded distribution can be
interpreted as the number of $\nu_{e}$ events that interact in MiniBooNE. One
can take the ratio of the unfolded MiniBooNE $\nu_{e}$ prediction to the
nominal $\nu_{e}$ prediction from the MiniBooNE Monte Carlo (MC), as shown in
figure 4.1(a). This gives a set of binned weights as a function of the true
$\nu_{e}$ energy, which can be applied to any BNB MC sample of $\nu_{e}$
interactions to represent a MiniBooNE-like excess. These weights constitute
the MicroBooNE eLEE model of the MiniBooNE anomaly. Figure 4.1(b) shows the
actual eLEE model weights used in all three MicroBooNE $\nu_{e}$ analyses,
calculated using the first $12.84\times 10^{20}$ POT of the MiniBooNE neutrino
mode dataset.

A few comments are necessary regarding the unfolding procedure. First, the
number of iterations is a regularization parameter that must be selected by
the user. MicroBooNE used three iterations of the D’Agostini algorithm to
produce the eLEE model, a choice designed to minimize the variance and bias in
the unfolded spectrum while remaining in good agreement with the MiniBooNE
data, i.e. retaining a $\chi^{2}$ per degree-of-freedom less than 1 in the
reconstructed distribution. While it is true that these criteria select a
unique choice for the number of iterations, the eLEE model shown in figure
4.1(b) is not the unique unfolded distribution consistent with the MiniBooNE
excess. This point is discussed further in Ref. [188]. Additionally, the eLEE
model was developed before MiniBooNE had looked at any $\nu_{e}$ data in the
$E_{\nu}^{\rm QE}<200\;{\rm MeV}$ region [79]. Thus, the unfolded $\nu_{e}$
prediction will lose events that reconstruct below 200 MeV. It is also worth
noting that while the MiniBooNE $e$-like dataset is only sensitive in
principle to $\nu_{e}$ CC interactions without pions in the final state (as
both charged and neutral pions are visible in a Cherenkov detector), the eLEE
model weights are applied to all $\nu_{e}$ CC interactions in MicroBooNE. The
impact of this choice is mitigated by the development of three different
$\nu_{e}$ analyses with an emphasis on different final states. If the
MiniBooNE excess only exists in one class of $\nu_{e}$ interactions, this
should be evident when comparing results from all three analyses. Finally, the
eLEE model does not account for systematic uncertainties on the MiniBooNE
prediction. This is important, as MiniBooNE is a systematics-limited
experiment [10]. The eLEE model weights shown in figure 4.1(b) can be thought
to represent a “median” model of the MiniBooNE excess. One can scale the
weights by an overall factor, hereafter referred to as the “LEE signal
strength” $x_{\rm LEE}$, to capture the 21% statistical $\mathop{\oplus}$
systematic normalization uncertainty on the MiniBooNE prediction. In chapter
5, we will cover results from a one-dimensional signal strength scaling test.

(a) From Ref. [20]

(b)

Figure 4.1: Figure 4.1(a) shows the unfolded $\nu_{e}$ prediction in MiniBooNE
calculated using the first $6.46\times 10^{20}$ POT of neutrino mode data
[20]. Figure 4.1(b) shows the unfolded eLEE model weights derived from the
first $12.84\times 10^{20}$ POT of MiniBooNE neutrino mode data, which
constitute the MicroBooNE eLEE model.

###  4.3 Philosophy Behind the Two-Body CCQE Analysis

The two-body CCQE analysis seeks to do what MiniBooNE could not–isolate a
high-purity sample of $\nu_{e}$ interactions with a low SM photon background,
especially at low energies. To do this, we have selected a very specific final
state definition: $\nu_{e}$ CCQE interactions with reconstructed kinematic
variables that are consistent with two-body scattering.

The $\nu_{\ell}$ CCQE interaction channel is a simple two-body scattering
process,

$\nu_{\ell}n\to\ell p.$ (4.1)

It was chosen because CCQE is the dominant neutrino interaction at lower
neutrino energies, where the eLEE model weights in figure 4.1(b) are the
highest. As shown in figure 4.2, CCQE interactions comprise 77% of the
$\nu_{e}$ CC interactions in MicroBooNE below 500 MeV [148]. Other interaction
channels, such as meson exchange current (MEC) and resonant (RES) scattering,
become more significant at higher neutrino energies. These channels are much
more complicated than CCQE and thus come with larger modeling uncertainties.
MEC events involve pion exchange between two nucleons and are thus sometimes
referred to as “two-particle-two-hole” interactions [189]. RES events involve
the excitation of a nucleon to a $\Delta$ resonance which will subsequently
decay, most often to a pion and nucleon. CCQE events are also easier to
reconstruct; the original neutrino energy is simply the sum of the observed
electron and proton kinetic energies, up to smearing from nuclear effects. In
MEC and RES events, one must account for the kinetic energy of additional
final state protons and pions, and the potential existence of invisible final
state neutrons makes it more difficult to connect the total kinetic energy of
the observed final state to the initial neutrino energy. Thus, the relative
simplicity of CCQE interactions is another advantage of choosing such a signal
definition.

The two-body CCQE analysis has developed a number of techniques to isolate
$\nu_{e}$ CCQE events from other $\nu_{e}$ and $\nu_{\mu}$ interactions in
MicroBooNE. The most basic of these involves leveraging the clean $1e1p$ final
state topology of the signal events. The left panel of figure 4.3 shows an
example candidate $1e1p$ event in MicroBooNE data. One can clearly see the
proton track and electron EM shower emanating from the same vertex, presumably
the location of the neutrino interaction in the LArTPC. Our reconstruction
algorithm takes advantage of this distinct “vee” shape, imposing a strict
requirement that exactly two “prongs” of charge (i.e., tracks or showers) come
from the same 3D-consistent interaction vertex. This process is augmented by a
powerful convolutional neural network to specifically find events with one
track and one shower, indicative of a $1e1p$ final state. The reconstruction
algorithm for this analysis will be discussed in more detail in the next
section.

We also rely on two-body kinematics to select $\nu_{e}$ CCQE interactions. In
theory, the observed $1e1p$ final state from a $\nu_{e}$ CCQE interaction
should obey kinematic constraints imposed in a two-body scattering process. In
a realistic MicroBooNE scenario, the struck neutron lives within the argon
nuclear medium and has non-zero initial momentum. The final state proton must
also travel through the argon nuclear medium, in which it can undergo final
state interactions that alter its momentum and/or create additional hadronic
particles. In order to optimize our reconstruction and minimize uncertainties,
we are interested in events with minimal final state interactions. We also
want to remove backgrounds that will be especially inconsistent with two-body
kinematics, such as $\nu_{e}$ MEC events and $\nu_{\mu}$ events with a
$\pi^{0}$ in the final state in which one of the photons is reconstructed
alongside a proton track. This is achieved by introducing variables to the
analysis that are especially sensitive to deviations from two-body kinematics.
For example, the total transverse momentum of the neutrino interaction
($p_{T}$) will be further from zero for non-CCQE events. Additionally, as
shown in figure 4.4, the proton in a two-body CCQE scattering interaction will
always be forward-going. While nuclear effects can in principle break this
requirement in MicroBooNE, it is much more likely for a non-CCQE background
event to result in a backward-going proton than a true CCQE event. Therefore,
we explicitly require forward-going protons in this analysis. We can also
consider the Bjorken scaling variable

$x_{Bj}=\frac{Q^{2}}{2p_{p}\cdot q},$ (4.2)

where $q^{2}=-Q^{2}$ is the four-momentum transfer and $p_{p}$ is the final
state proton four-momentum. In the case of quasi-elastic scattering, $x_{Bj}$
should be close to unity up to nuclear effects [190], while it can
significantly differ from unity for non-CCQE backgrounds [148].

Finally, we can harness the full power of the LArTPC technology to reconstruct
the $\nu_{e}$ energy in multiple ways. For $\nu_{e}$ traveling along the $z$
axis of the detector, the neutrino energy can be calculated from any of the
following,

$\displaystyle E_{\nu}^{\rm range}$
$\displaystyle=\textnormal{K}_{p}+\textnormal{K}_{\ell}+M_{\ell}+M_{p}-(M_{n}-B),$
(4.3) $\displaystyle E_{\nu}^{QE-p}$
$\displaystyle=\left(\dfrac{1}{2}\right)\/\dfrac{2\cdot(M_{n}-B)\cdot
E_{p}-((M_{n}-B)^{2}+M_{p}^{2}-M_{\ell}^{2})}{(M_{n}-B)-E_{p}+\sqrt{(E_{p}^{2}-M_{p}^{2})}\cdot\cos\theta_{p}},$
(4.4) $\displaystyle E_{\nu}^{QE-\ell}$
$\displaystyle=\left(\dfrac{1}{2}\right)\/\dfrac{2\cdot(M_{n}-B)\cdot
E_{\ell}-((M_{n}-B)^{2}+M_{\ell}^{2}-M_{p}^{2})}{(M_{n}-B)-E_{\ell}+\sqrt{(E_{\ell}^{2}-M_{\ell}^{2})}\cdot\cos\theta_{\ell}},$
(4.5)

where $K_{\ell/p}$ is the kinetic energy of the lepton or proton,
$\theta_{\ell/p}$ is the angle of the shower or track with respect to the $z$
axis, $M_{\ell/p/n}$ is the particle mass, and $B$ is the average binding
energy, which is around 40 MeV in argon [191]. For a true $\nu_{e}$ CCQE
interaction, all three of these formulae should yield a similar reconstructed
neutrino energy. To leverage this effect we define the “QE consistency”
variable $\Delta^{QE}$ as

$(\Delta^{QE})^{2}\equiv(E_{\nu}^{\rm range}-E_{\nu}^{QE-p})^{2}+(E_{\nu}^{\rm
range}-E_{\nu}^{QE-\ell})^{2}+(E_{\nu}^{QE-\ell}-E_{\nu}^{QE-p})^{2},$ (4.6)

which should be close to zero for CCQE events. To mitigate smearing from
nuclear effects, $E_{\nu}^{QE-p}$ and $E_{\nu}^{QE-\ell}$ are calculated in
the rest frame of the struck nucleon. We note here that equation 4.3 is the
official reconstructed neutrino energy used for this analysis; $E_{\nu}$ and
$E_{\nu}^{\rm range}$ are used interchangeably throughout this thesis.

Finally, it is worth mentioning that the two-body CCQE analysis also isolates
a control sample of $\nu_{\mu}$ CCQE interactions with a $1\mu 1p$ final
state. The $1\mu 1p$ sample is constructed analogously to the signal $1e1p$
sample, though with around two orders of magnitude more events (as $\nu_{\mu}$
comprise most of the BNB flux). This sample is used to constrain the intrinsic
$\nu_{e}$ prediction and uncertainties in the $1e1p$ sample, as the $\nu_{e}$
and $\nu_{\mu}$ flux and cross section are highly correlated. We also make use
of dedicated $\pi^{0}$ and Michel electron samples, as will be discussed in
section 4.4.3.

In summary, the goal of the analysis presented here is the isolation of clean
two-body $\nu_{e}$ CCQE interactions in MicroBooNE. The following sections
explain how we leverage the topological and kinematic patterns described above
to identify such events.

Figure 4.2: The expected distribution of $\nu_{e}$ interactions in MicroBooNE
as a function of the true $\nu_{e}$ energy. The dotted line shows the
expectation from the MicroBooNE eLEE model of the MiniBooNE excess discussed
in section 4.2. Figure 4.3: An example candidate $1e1p$ data event in
MicroBooNE, including the raw LArTPC collection plane image (left) and the
pixel labels assigned from SparseSSNet (right). Figure 4.4: The relationship
between the final state lepton (left) and proton (right) energy and scattering
angle, for different neutrino energies, in $\nu_{\mu}$ (top) and $\nu_{e}$
(bottom) CCQE scattering.

###  4.4 Reconstruction

We begin with a broad overview of the reconstruction chain for the two-body
CCQE analysis. An initial set of neutrino event candidates are identified by
requiring sufficient light detection in time with the beam trigger, as
discussed in section 3.1.3. These events are passed through the noise filter
and deconvolution algorithm described in section 3.2 and shown in figure 3.7.
Additionally, a “good runs cut” is employed to remove events during which the
MicroBooNE detector response was especially noisy [148]. These steps are
common to all MicroBooNE analyses. They result in one event display per wire
plane for every event, similar to the rightmost panel of figure 3.7 and the
left panel of figure 4.3. The color axis in the left panel of figure 4.3 is
given in “pixel intensity units” (PIU), which are defined as the integrated
charge in one wire over six 0.5 $\mu$s time ticks [148]. Given the electron
drift speed of 0.11 cm/$\mu$s, this corresponds to a spatial extent of 0.33 cm
along the drift direction. The wire spacing is such that each pixel in the
left panel of figure 4.3 represents a $0.3$ cm $\times$ $0.33$ cm 2D square in
the detector.

In the two-body CCQE analysis, these event displays are first passed through
the sparse semantic segmentation network (SparseSSNet) [21]. SparseSSNet is a
custom convolutional neural network (CNN) developed by the analysis team that
labels the active pixels in each event according to five different charge
deposition hypotheses: highly-ionizing particle (HIP), minimum-ionizing
particle (MIP), EM activity, delta electron, or Michel electron. For the
purposes of this analysis, we group HIP and MIP labels into the parent “track”
label and EM activity, delta electrons, and Michel electrons into the parent
“shower” label. An example SparseSSNet label is shown in the right panel of
figure 4.3.

After going through SparseSSNet, the images are filtered for crossing cosmic
interactions using the WireCell charge-to-light matching algorithm [192, 193].
This algorithm matches charge observed on the TPC to light observed by the
PMTs to identify regions of charge which are not coincident in time with the
beam spill, presumably coming from cosmic rays. We can then mask these pixels
to remove the bulk of cosmic-related charge deposition. The cosmic-filtered
SparseSSNet output is then used as input to a vertex and track identification
algorithm, which looks for a shower-like and track-like prong of charge coming
from a given 3D vertex [24]. Track-like prongs pass through a dedicated step
of the algorithm that performs a step-by-step 3D reconstruction of the extent
and direction of the track, which is then translated into a range-based energy
and angle estimation. The events are then passed through a dedicated EM shower
identification algorithm that identifies cones of shower-like pixels pointing
back to the neutrino interaction vertex and estimates the total visible energy
and angle of the shower. [172]. This algorithm has the capability of finding
up to two showers separated from the vertex, which is necessary for
identifying events with a $\pi^{0}\to\gamma\gamma$ decay.

After finding candidate neutrino vertices, we pass the event displays through
the multiple particle identification (MPID) CNN [22]. This network assigns
five scores to each image indicating the network’s confidence that the event
contains a proton, electron, muon, photon, or charged pion. The MPID scores
are not used until the final signal selection described in section 4.5.

At this point, the reconstruction chain has identified two-prong candidate
events and reconstructed the relevant lowest-level kinematic variables of each
final state particle, i.e. the energy and direction. These kinematic variables
are then combined to calculate more complicated kinematic variables such as
$Q^{2}$, $p_{T}$, $x_{Bj}$, and the reconstructed $\pi^{0}$ mass (in the case
of two-shower events), as well as the most important variable in this
analysis: the reconstructed neutrino energy $E_{\nu}^{\rm range}$, as defined
in equation 4.3. The full suite of reconstructed kinematic variables used to
define the $1e1p$ and $1\mu 1p$ samples is given in table 4.1. These variables
are used as input to the signal selection described in section 4.5. The
following subsections explain in more detail the most important steps of the
reconstruction procedure.

Variable Name | Definition
---|---
Base Variables
$\textnormal{K}_{p}$ | Kinetic energy of proton determined from range [24]
$\textnormal{K}_{\ell}$ | Kinetic energy determined from range for muons [24] and from calorimetry for electrons [172]
$M_{\ell}$, $M_{n}$, $M_{p}$ | Masses of the lepton, neutron and proton
$\cos{\theta_{p}}$, $\cos{\theta_{\ell}}$ |  $p_{p}^{z}/p_{p}$, $p_{\ell}^{z}/p_{\ell}$
$\phi_{p}$, $\phi_{\ell}$ |  $atan2(p_{p}^{y},p_{p}^{x})$, $atan2(p_{\ell}^{y},p_{\ell}^{x})$
$P_{p}=(E_{p},p_{p})$, $P_{\ell}=(E_{\ell},p_{\ell})$ | Reconstructed 4-vector of the proton, lepton
$B$ | Binding Energy for argon; the analysis assumes $B=40$ MeV
Definitions Related to Neutrino Energy
$E_{\nu}^{range}$ (Default value of $E_{\nu}$) | $\textnormal{K}_{p}+\textnormal{K}_{\ell}+M_{\ell}+M_{p}-(M_{n}-B)$
$E_{\nu}^{QE-p}$ | $\big{(}\dfrac{1}{2}\big{)}\/\dfrac{2\cdot(M_{n}-B)\cdot E_{p}-((Mn-B)^{2}+M_{p}^{2}-M_{\ell}^{2})}{(M_{n}-B)-E_{p}+\sqrt{(E_{p}^{2}-M_{p}^{2})}\cdot\cos\theta_{p}}$
$E_{\nu}^{QE-\ell}$ | $\big{(}\dfrac{1}{2}\big{)}\/\dfrac{2\cdot(M_{n}-B)\cdot E_{\ell}-((Mn-B)^{2}+M_{\ell}^{2}-M_{p}^{2})}{(M_{n}-B)-E_{\ell}+\sqrt{(E_{\ell}^{2}-M_{\ell}^{2})}\cdot\cos\theta_{\ell}}$
$\Delta^{QE}$ (2-Body Consistency) | $\sqrt{(E_{\nu}^{range}-E_{\nu}^{QE-p})^{2}+(E_{\nu}^{range}-E_{\nu}^{QE-\ell})^{2}+(E_{\nu}^{QE-\ell}-E_{\nu}^{QE-p})^{2}}$
Event Kinematics
$Q^{2}$ | $2E_{\nu}^{range}(E_{\ell}-P_{\ell}^{z})-M_{\ell}^{2}$
Hadronic Energy (Ehad) | $E_{\nu}^{range}-E_{\ell}$
Björken’s Scaling x (xBj) | $Q^{2}/2M_{n}E_{had}$
Björken’s Scaling y (yBj) | $E_{had}/E_{\nu}^{range}$
Opening angle | $\cos^{-1}(\hat{p}_{\ell}\cdot\hat{p}_{p})$
$p_{T}$ | $\sqrt{(p_{\ell}^{x}+p_{p}^{x})^{2}+(p_{\ell}^{y}+p_{p}^{y})^{2}}$
$p_{L}$ | $p_{p}^{z}+p_{\ell}^{z}$
$\alpha_{T}$ | $\cos^{-1}\left(-\dfrac{\vec{P}_{T}^{l}\cdot\vec{P}_{T}}{|\vec{P}_{T}^{l}||\vec{P}_{T}|}\right)$
$\phi_{T}$ | $\cos^{-1}\left(-\dfrac{\vec{P}_{T}^{l}\cdot\vec{P}_{T}^{p}}{|\vec{P}_{T}^{l}||\vec{P}_{T}^{p}|}\right)$
Boosting Parameters
$p_{fermi}^{T}$ | $p_{p}^{T}+p_{\ell}^{T}$
$p_{fermi}^{z}$ | $p_{p}^{z}+p_{\ell}^{z}-E_{\nu}^{range}$
Table 4.1: The definition of kinematic variables used throughout the two-body
CCQE analysis.

####  4.4.1 Convolutional Neural Networks in LArTPCs

As discussed above, The two-body CCQE analysis relies on two CNNs to assist in
reconstructing the MicroBooNE LArTPC images: SparseSSNet and MPID.

SparseSSNet is described in detail in Ref. [21]. The network uses a
combination of the U-Net [194] and ResNet [195] architectures. The U-Net uses
an encoding step, in which the input image is downsampled to a smaller size,
followed by a decoding step, in which the encoded image is upsampled to the
same dimensionality and size as the original image. This second step is
required to perform pixel-level labeling. The ResNet architecture employs
residual connections that improve the training performance of deeper networks
compared to regular CNNs [195]. SparseSSNet also makes use of sparse
submanifold convolutions [196], which perform better on semantic segmentation
tasks for sparse input data like LArTPC images (where very few pixels in a
given image have nonzero charge). The training sample is generated using the
simulation described in section 4.1; specifically, a random set of visible
final state particles are generated at randomly distributed locations within
the detector. The network performs considerably well on the nominal MicroBooNE
BNB simulation sample, as well as the dedicated intrinsic $\nu_{e}$ simulation
sample, neither of which appeared in the SparseSSNet training sample [21].
While SparseSSNet has not been statistically evaluated against MicroBooNE
data, hand scans of the network output on data event displays do not reveal
any bias in the network performance [21]. The performance of an earlier
iteration of SparseSSNet was evaluated using two data samples: Michel
electrons and $\nu_{\mu}$ CC $\pi^{0}$ events [197]. Good agreement was
observed between pixel labels in data and simulation for both samples, and no
significant differences were found between network-labeled images and
physicist-labeled images. Figure 4.5 shows the SparseSSNet architecture as
well as an example labeled image in the training dataset.

The MPID network is described in detail in Ref. [22]. It uses a standard CNN
architecture which takes as input a 2D LArTPC image cropped around the
neutrino interaction candidate vertex on a single wire plane and outputs five
numbers representing different particle scores for protons, electrons, muons,
photons, and charged pions. The score is meant to represent the likelihood
that a given particle exists in the provided image. Two types of scores are
generated: an “image” score, describing the probability that the particle
exists anywhere in the image, and the “interaction” score, describing the
probability that the particle is connected to the neutrino interaction vertex.
The training sample is generated similarly to that of SparseSSNet; random sets
of particles are generated throughout the detector volume with energies in a
range relevant to MicroBooNE. The MPID score distributions exhibit excellent
agreement between data and simulation in two different samples: $\nu_{\mu}$ CC
$1\mu 1p$ events and $\nu_{\mu}$ CC $\pi^{0}$ events [22]. Figure 4.6 shows
the MPID architecture as well as an example simulated image with MPID image
scores.

It is also worth briefly mentioning energy reconstruction in LArTPCs using
CNNs. My earliest work as a graduate student investigated the ability of the
ResNet [195] and Inception network [198] CNN models to reconstruct EM shower
energies in MicroBooNE. Compared to the standard clustering-based shower
reconstruction algorithm described in section 4.4.3, the CNN method appeared
to better handle EM showers that passed through unresponsive wires in the
MicroBooNE detector. Despite the promise of the CNN method, it was not
developed in time to be implemented in the two-body CCQE analysis.
Nevertheless, the MicroBooNE-specific network was extended by an MIT
undergraduate, Kiara Carloni, to enable a systematic study of CNN-based EM
shower energy reconstruction in a general LArTPC [23]. The corresponding
publication is included in Appendix A. Figure 4.7 shows the main takeaway of
this work: a CNN-based reconstruction method is significantly more robust to
unresponsive wires in the detector compared with a traditional EM shower
reconstruction algorithm, which employs a linear calibration between the total
shower charge and energy.

(a)

(b)

Figure 4.5: Figure 4.5(a) shows a diagram of the SparseSSNet U-ResNet
architecture. Figure 4.5(b) shows the SparseSSNet pixel labels on a simulated
image in the training dataset. Figures from Ref. [21].

(a)

(b)

Figure 4.6: Figure 4.6(a) shows a diagram of the MPID architecture. Figure
4.6(b) shows the MPID image scores on an example simulated event. Figures from
Ref. [22]. Figure 4.7: The fraction of EM showers reconstructed to within 5%
of the true deposited energy as a function of the fraction of unresponsive
wires in the LArTPC collection plane, considering three different methods. The
ResNet and Inception network significantly outperform the traditional linear
calibration between charge and energy. Figure from Ref. [23].

####  4.4.2 Vertex and Track Reconstruction

After running SparseSSNet and filtering cosmic-associated charge, the images
are passed through the vertex reconstruction algorithm described in Ref. [24].
In brief, the algorithm looks for the characteristic “vee” shape of both
$1e1p$ signal events, as shown in the left panel of figure 4.3, and $1\mu 1p$
control sample events. The algorithm begins by isolating a sample of 2D vertex
seeds (i.e., candidate neutrino interaction locations) on each plane. This is
done by clustering regions of charge and then applying two techniques–contour
defect identification using a convex hull and intersection identification
using principal component analysis. The set of vertex seeds is reduced to
those that are consistent in the time dimension across all three planes. From
each triplet of time-consistent vertex seeds on the three wire planes, a final
3D vertex is found by minimizing an angular metric that represents the
likelihood that two clusters of charge emanate radially from that vertex
location. Figure 4.8(a) shows an example minimization of this angular metric,
which involves calculating two different estimates of the opening angle
between the two prongs for a given vertex location. Simulation studies
indicate that this algorithm is able to reconstruct neutrino interaction
vertices to within less than 1 cm with an average efficiency of 56% [24].

After finding a neutrino vertex candidate, we perform a 3D reconstruction of
the prongs of charge emitted from the vertex [24]. An iterative stochastic
algorithm identifies a set of 3D spatial points associated with each prong by
building upon previous points; these points comprise the track. This process
is shown diagrammatically in figure 4.8(b). It is worth noting that the
“tracks” at this stage can correspond to either true tracks (HIPs and MIPs) or
EM showers.

The 3D track can be projected onto each 2D plane. At this point, the
SparseSSNet labels on each plane are used to determine if the event is a
“track-track” or “track-shower” event. The latter case happens when greater
than 20% of the pixels in a given prong are labeled as shower pixels. In the
“track-shower” case, the shower-like prong is considered to be the electron
while the track-like prong is the proton. Thus, these are candidate $1e1p$
events. The “track-track” case corresponds to candidate $1\mu 1p$ event, where
the shorter track is ascribed to the highly-ionizing proton. Electron EM
showers are not sufficiently reconstructed in this step and are left for a
dedicated shower reconstruction algorithm described in section 4.4.3. Muons
and protons can, however, be reconstructed at this stage. The polar and
azimuth angles of each particle can be determined from the direction of the
associated track. The length of the tracks is used in conjunction with the
known stopping power of each particle in liquid argon to determine the kinetic
energy. Simulation studies of $1\mu 1p$ events suggest an angular resolution
of a few degrees and a kinetic energy resolution of a few percent [24].

It is also worth mentioning that two-prong events with an opening angle of
less than 10∘ in any plane are rejected at this point. Thus, signal events in
this analysis consisted of well-separated leptons and protons.

(a)

(b)

Figure 4.8: Figure 4.8(a) shows the angular metric that is minimized to find a
3D neutrino vertex candidate. Figure 4.8(b) shows the iterative track
reconstruction algorithm, which relies on calculating distances (L1 and L2)
and angles ($\theta$ and $\phi$) with respect to the previous point and the
end of the track. Figures from Ref. [24].

####  4.4.3 Publication: Electromagnetic shower reconstruction and energy
validation with Michel electrons and $\pi^{0}$ samples for the deep-learning-
based analyses in MicroBooNE

As mentioned in the previous section, a dedicated algorithm is used to
reconstruct the EM showers coming from the neutrino interaction vertex. The
algorithm is described in detail in Ref. [172], which also covers a data-
driven validation of the shower reconstruction using dedicated $\pi^{0}$ and
Michel electron samples. A brief overview of the shower reconstruction is as
follows. We first mask out pixels in each plane for which the SparseSSNet
shower score is below a value of 0.5. Next, we attach a triangle to the
interaction vertex in each plane and optimize the direction, gap from the
vertex, opening angle, and length of the triangle to encapsulate as many
shower pixels as possible without becoming too large. We calculate the energy
of the shower by performing a linear calibration between the summed PIU in the
collection plane shower triangle and the simulated electron/photon energy. The
collection plane is used because it produces the most robust signals, as
discussed in section 3.2. We reconstruct 3D showers by looking for overlap in
the time dimension between the shower triangles on all three planes. This
algorithm can be run a second time after masking out the optimized triangle
from the first pass, which is essential for isolating a $\pi^{0}$ sample.

This MicroBooNE publication was led by myself and Katie Mason, a graduate
student from Tufts University. Katie focused on the triangle-fitting procedure
of the shower reconstruction and the $\pi^{0}$ mass peak fit, while I focused
on the charge-to-energy scaling of the shower reconstruction and the Michel
energy spectrum fit. The full JINST publication is included below. The most
important result from this study is shown in Figure 16 of the paper, which
demonstrates good agreement between the simulation-derived charge-to-energy
scaling applied to EM showers in this analysis and two standard candles: the
$\pi^{0}$ invariant mass peak and the cutoff of the Michel electron energy
distribution.

See pages - of templates/Publications/Abratenko_2021_J._Inst._16_T12017.pdf

###  4.5 $1e1p$ Event Selection

The signal events for this analysis, $1e1p$ events consistent with two-body
CCQE scattering, are isolated using a series of cuts on kinematic,
topological, and MPID-based variables as well as an ensemble of boosted
decision trees (BDTs). The $1e1p$ signal selection represents one of my main
projects within the MicroBooNE experiment. This section covers in detail the
requirements which define the $1e1p$ sample. The $1\mu 1p$ control sample is
isolated analogously to the $1e1p$ sample and will not be covered in this
thesis, though more information can be found in Ref. [148].

####  4.5.1 Basic Data Selection Criteria

We first apply a series of basic “precuts” meant to remove poorly
reconstructed events and obvious backgrounds. We begin with the collection of
“track-shower” vertices selected by the algorithm described in section 4.4.2,
which represent candidate $1e1p$ events. We remove any vertices that are
within 10 cm of the TPC edges or within a region of unresponsive wires between
$700<z~{}[{\rm cm}]<740.$ This cut defines the fiducial volume. We also
require that the minimal distance of either the track or the shower to the
edge of the TPC is greater than 15 cm. In the case that both prongs are
reconstructed close to an edge, we instead require that the distance of each
is greater than 5 cm, which helps retain neutrino interactions that begin near
the edge of the fiducial volume. Events in the $1e1p$ are further required to
be outside of an inefficient region of the U-plane.

We also employ a cut requiring the reconstructed shower energy of the
candidate electron to be consistent between the three planes. Even though we
rely on the collection plane for the official shower energy estimation, this
cut is helpful in removing events that would not lead to a consistent 3D
picture of the shower. An example of such an event is shown in figure 4.9. The
cut specifically employed in this analysis is

$\frac{\sqrt{(E_{e}^{U}-E_{e}^{V})^{2}+(E_{e}^{U}-E_{e}^{Y})^{2}+(E_{e}^{V}-E_{e}^{Y})^{2}}}{E_{e}^{Y}}<2,$
(4.7)

where $E_{e}^{[U,V,Y]}$ denote the reconstructed electron candidate energy on
the $U$, $V$, and $Y$ plane, respectively. Figure 4.10(a) shows the
distribution of the fractional consistency variable (left-hand side of
equation 4.7) for events above and below a $1e1p$ BDT score cutoff of 0.7,
where the BDT here is an older iteration of the final BDT ensemble used in
this analysis. As can be seen, signal-like events (with a higher BDT score)
tend to have lower fractional consistency values compared to background-like
events (with a lower BDT score). The upper bound of 2 in equation 4.7 was
chosen to retain $\sim 95\%$ of signal-like events while rejecting $\sim 20\%$
of background-like events, as shown in figure 4.10(b).

Next, we make cuts on the reconstructed neutrino energy ($200<E_{\nu}~{}[{\rm
MeV}]<1200$), electron kinetic energy ($K_{e}>35~{}{\rm MeV}$), and proton
kinetic energy ($K_{p}>50~{}{\rm MeV}$). We also require a forward-going
proton and an opening angle between the electron and proton greater than 0.5
rad. Finally, we require that the event can be boosted to the rest frame of
the nucleon using the $p_{fermi}^{T}$ and $p_{fermi}^{z}$ variables defined in
table 4.1.

Figure 4.9: An example of an event that fails the shower energy consistency
cut, because the EM shower passes through an unresponsive region of the
collection plane.

(a)

(b)

Figure 4.10: Figure 4.10(a) shows the distribution of the fractional shower
energy consistency variable for events with an old $1e1p$ BDT score
above/below 0.7. Figure 4.10(b) shows the efficiency with which events
above/below the old $1e1p$ BDT cutoff pass the fractional consistency cut as a
function of the chosen upper bound.

####  4.5.2 Boosted Decision Tree Ensemble

The next stage of the $1e1p$ selection uses an ensemble of BDTs to identify
signal-like events in MicroBooNE. This is the most powerful cut used in the
$1e1p$ analysis. The ensemble consists of 20 BDTs each trained using the
XGBoost gradient-boosting algorithm [199]. Each BDT takes as input a set of 19
kinematic variables, such as the reconstructed energies and angles of the
electron and proton, and 4 topological variables related to the observed
charge, such as the fraction of SparseSSNet-labeled shower pixels in the
event. The full suite of variables used to train the $1e1p$ and $1\mu 1p$ BDTs
is shown in table 4.2. Each BDT itself consists of a collection of regression
trees which each assign a continuous score to an event based on a set of
decision rules developed using the input variables; the assigned BDT score of
the event is given by the sum over the collection [199]. BDTs are trained via
gradient-based optimization of an objective function that includes a
regularization term to penalize overfitting.

A different ensemble is generated for each of the three run periods, as the
detector response changed slightly between them. Thus, there are 60 total BDTs
trained for the $1e1p$ analysis. The BDTs are trained using the MicroBooNE MC
samples, including the nominal ($\nu_{\mu}$-dominated) BNB sample as well as
dedicated intrinsic $\nu_{e}$, $\pi^{0}$, and cosmic samples. The signal for
the training sample is defined as true $\nu_{e}$ CCQE interactions with one
electron and one proton in the final state. We further restrict to well-
reconstructed events for which the vertex is reconstructed with 5 cm of the
true vertex and the neutrino energy within 20% of the true value. Only
non-$\nu_{e}$ CCQE $1e1p$ events are explicitly classified as background. Each
BDT in an ensemble is trained using a randomly-selected half of the MC sample
for the given run period. Figure 4.11 shows the F score distribution for each
variable in one of the $1e1p$ BDTs from the ensemble for each run period. The
F score represents the number of times each variable is used across all trees
in the BDT; thus, it can be thought of as a measure of the importance of each
variable.

The $1e1p$ score for an event is evaluated by calculating the average BDT
score across the ensemble. This reduces the dependence of the ensemble score
on the specific training sample of any single BDT. Ensemble-based methods have
been shown to reduce variance on the output likelihood variable compared to a
single classifier [200]. When computing the average score for events in the
MicroBooNE MC sample, we omit BDTs for which the event appeared in that BDT’s
training sample. It is extremely unlikely ($p=2^{-20}$) for an event to appear
in the training sample of all BDTs in an ensemble; thus, we can use the BDT
ensemble to isolate a $1e1p$ signal sample prediction from the MicroBooNE MC
without throwing away any events. The score is normalized to the range [0,1],
where higher scores indicate a higher signal likelihood. We show the
distribution of the ensemble-averaged $1e1p$ BDT score in figure 4.12, which
shows that true $\nu_{e}$ CCQE events tend to peak toward 1 while all other
events peak toward 0. Signal events are required to have an average score
greater than 0.95, which was chosen to optimize sensitivity to the eLEE model
described in section 4.2.

We have performed a series of tests of the robustness of the BDT ensemble
method. We first examined the signal selection power of a BDT ensemble trained
on simulated events from a given run period when used to infer the signal
likelihood of simulated events from a different run period. This is useful for
two reasons. For one, it provides a method for testing the ensemble on events
that were in neither the training nor the validation sample of the constituent
BDTs. Testing on a sample independent of the validation sample is important
because the BDT training was halted when the classification error on the
validation set didn’t improve after 50 training iterations. It also tests the
impact of removing training events from the MC prediction. When using the BDT
ensemble of a given run period to evaluate the signal likelihood of events
from the same run period, one must remove events that appeared in the training
sample of each constituent BDT. In this study we use the BDT ensemble trained
on run period 2 to evaluate simulated $\nu_{e}$ events from run period 3 and
vice versa, so one does not need to remove training events.

Figure 4.13 shows the results of this run period swap study. As one can see,
the differences between the selection using the correct and incorrect run
period ensemble are small. The correct run period performs slightly better, as
expected due to small differences in the detector status specific to each run
period. However, the relatively small change in $\nu_{e}$ selection efficiency
suggests that the BDT ensemble is able to perform well on simulated events
that do not appear in either the training or validation sample of constituent
BDTs. One can interpret the loss in performance shown in figure 4.13 as a
bound on the potential impact of applying the BDT ensemble to events that are
not in the training or validation samples of the entire ensemble. The drop in
performance is $<10\%$ over most of the energy range–an effect which is
smaller than our systematic error before the $\nu_{\mu}$ constraint described
in chapter 5 ($\sim 15\%$) and much smaller than the statistical error in each
reconstructed neutrino energy bin after the final selection ($\sim$ 50%).

Next, we investigated the impact of using only a subset of the BDTs in the
ensemble to calculate the average BDT score of an event. This is relevant
because we omit BDTs that contained a given MC event in the training sample
when calculating the average BDT score for that event. However, this analysis
uses all 20 BDTs in the ensemble to calculate the average BDT score of events
in the data. Therefore, to verify that comparisons between data and simulation
are robust, one needs to ensure that removing BDTs from the ensemble does not
significantly bias the average score calculation.

To this end, consider $S_{n}$ to be the $1e1p$ BDT ensemble average score
after removing $n$ BDTs from the ensemble. For events in data, the BDT score
is $S_{0}$ while for simulated events, the BDT score is $S_{n}$ for some
$n\in\\{1,...,20\\}$. Figure 4.14 shows the fractional difference
$(S_{n}-S_{0})/S_{0}$ as a function of the number of omitted BDTs $n$ over
signal-like (BDT score $>$ 0.95) events in the simulation from run period 2
and run period 3, respectively. As in the previous study, in order to avoid
bias from BDT training we use the BDT ensemble trained on run period 2 to
evaluate the signal likelihood of simulation events from run period 3 and vice
versa. The data points and error bars in figure 4.14 indicate the average and
standard deviation of the fractional difference over the simulation sample,
respectively. The red histograms show the actual distribution of omitted BDTs
over the simulation from each run period. One can see that the average BDT
score does not exhibit significant bias upon removing BDTs from the
calculation. Also, the standard deviation of the fractional difference only
becomes larger than $1$% when removing $\gtrsim 17$ BDTs from the calculation.
This only happens for $<<1$% of simulated events, as expected.

Variable | Used in $1\mu 1p$ BDT | Used in $1e1p$ BDT
---|---|---
Variables Used in BDTs, Based on Ionization
Charge within 5 cm of vertex | Yes | Yes
Shower charge in event image / |  |
shower charge clustered as electron | No | Yes
Proton shower fraction | No | Yes
Electron shower fraction | No | Yes
Variables Used in BDTs, Related to Energy Measurements
Neutrino Energy | Yes | Yes
Energy of electromagnetic shower | No | Yes
Lepton length | Yes | Yes
Proton length | No | Yes
$p_{z}-E_{\nu}$ | No | Yes
Variables Used in BDTs, Related to 2-Body Scattering Consistency
Bjorken’s $x$ | Yes * | Yes *
Bjorken’s $y$ | Yes * | Yes *
QE Consistency | Yes * | Yes *
$Q_{0}$ | Yes | Yes
$Q_{3}$ | Yes | Yes
Variables Used in BDTs, Related to Transverse Momentum
$\alpha_{T}$ | Yes | Yes
Event $p_{T}$ | Yes | Yes
Event $p_{T}/p$ (“PTrat") | Yes | Yes
$\phi_{T}$ | Yes | No
Variables Used in BDTs, Related to Angles
Proton $\phi$ | Yes | Yes
Proton $\theta$ | Yes | Yes
Lepton $\phi$ | Yes | Yes
Lepton $\theta$ | Yes | Yes
$\phi_{p}-\phi_{\ell}$ | Yes | Yes
$\theta_{p}+\theta_{e}$ | No | Yes
Variables Useful for Comparison, Not Used in Either BDT
$\eta$ (Normalized average ionization difference) | No | No
Opening Angle | No | No
$x$ Vertex | No | No
$y$ Vertex | No | No
$z$ Vertex | No | No
Table 4.2: The suite of variables used to isolate and analyze the $1e1p$ and
$1\mu 1p$ samples. Variables used in the BDT ensemble for each sample are
specified. The “$*$” character indicates that the variable is calculated in
the rest frame of the struck nucleon. The mathematical definitions of many of
these variables appear in Table 4.1.

(a) Run 1

(b) Run 2

(c) Run 3

Figure 4.11: The F score of each variable for one of the $1e1p$ BDTs in the
ensemble from run 1, run2, and run 3.

(a)

(b)

Figure 4.12: Figure 4.12(a) and figure 4.12(b) show the MC distribution of the
$1e1p$ BDT ensemble average score over all three run periods for the full
[0,1] range and zoomed in to the [0.95,1] range, respectively.

(a)

(b)

Figure 4.13: Figure 4.13(a) and figure 4.13(b) show the predicted $\nu_{e}$
event rate in run period 2 and run period 3, respectively, using both the run
period 2 ensemble and the run period 3 ensemble.

(a)

(b)

Figure 4.14: Figure 4.14(a) and figure 4.14(b) show the fractional difference
in average BDT score $(S_{n}-S_{0})/S_{0}$ as a function of the number of
omitted BDTs $n$ over the simulation from run period 2 and run period 3,
respectively. The red histogram shows the actual distribution of the number of
omitted BDTs over the run period 2 and run period 3 simulation samples,
respectively. Scores are calculated using the run period 3 and run period 2
BDT ensemble, respectively.

####  4.5.3 Particle Identification Cuts

Finally, we employ a series of particle-identification requirements to clean
up the remaining backgrounds that survive the $1e1p$ BDT ensemble cut. The
first of these is a cut on the invariant $\pi^{0}$ mass. If the shower
reconstruction algorithm described in section 4.4.3 is able to identify a
second EM shower, the reconstructed $M_{\pi^{0}}$ must be less than 50 MeV.
Another cut requires the ratio of the MPID $\gamma$ and $e^{-}$ image scores
to be less than 2. This helps remove remaining $\pi^{0}$ events that survive
the $M_{\pi^{0}}$ cut. We also require the MPID muon interaction score to be
less than 0.2, which helps remove $\nu_{\mu}$ CC$\pi^{0}$ events where the
muon gets mistaken for a proton. This last cut is only applied for events in
which $E_{e^{-}}>100$ MeV, as MPID tends to assign low electron scores and
high muon scores to true electrons below 100 MeV, as shown in figure 4.15.
This is because electrons stop radiating below this energy, becoming more
track-like and thus difficult to distinguish from muons. The specific values
used for the cuts described here were chosen by optimizing the signal-to-
background ratio in the $1e1p$ sample. Finally, if an event has more than one
reconstructed neutrino vertex candidate passing all selection criteria, we
keep only the vertex with the highest $1e1p$ ensemble BDT score.

(a)

(b)

Figure 4.15: Figure 4.15(a) and figure 4.15(b) show the MPID electron and muon
score, respectively, as a function of the reconstructed electron energy in
intrinsic $\nu_{e}$ MC events.

####  4.5.4 The Final $1e1p$ Sample

The requirements described throughout this section define the $1e1p$ sample;
they are summarized and compared with those in the $1\mu 1p$ sample in table
4.3. This results in a highly-pure $1e1p$ signal sample–75% of events are true
$\nu_{e}$ CCQE interactions. Figure 4.16 shows the predicted $E_{\nu}^{\rm
range}$ distribution for this sample, including the prediction from the eLEE
model. This figure also shows the systematic uncertainty and non-$\nu_{e}$
background prediction, which will be discussed further in chapter 5. The
neutrino energy distribution in figure 4.16 forms the basis of the statistical
results discussed in chapter 5.

The neutrino energy resolution is 16.5% for selected $1e1p$ events, as shown
in figure 4.17. This figure also indicates a slightly higher rate of events
with under-predicted neutrino energy, owing mainly to showers and tracks that
pass through unresponsive regions of the detector. Figure 4.18 reports the
signal efficiency and total number of $\nu_{e}$ CCQE events in this sample
after applying each of the three sets of cuts. The efficiency peaks at lower
neutrino energies $\sim 300$ MeV, which is optimal for testing the eLEE model.
Note that these efficiencies in figure 4.18(a) are reported with respect to
all identified vertices. If we include the effect of the vertex algorithm, the
average efficiency for retaining signal events in the $1e1p$ sample is 6.6%
[148]. This is by design–in choosing a very exclusive signal definition, we
optimize for a high-purity sample at the price of lower efficiency.

Name/Description of Code |  Requirement for $1\mu 1p$ |  Requirement for $1e1p$
---|---|---
Cuts/Selections Applied in External Code
Common Optical Filter |  $>20$ PMT hits w/i 6 ticks in beam window | Same
|  & $\leq 20$ hits w/i 6 ticks in 2 $\mu s$ prior to spill |
WC tagger | Mask cosmic-associated charge | Same
Preselection Cuts for Event Selection
Run Quality Requirements: |  |
Good Run Flag | True | Same
Vertex Requirements: |  |
$1\ell 1p$ consistent | Two prongs of at least 5cm length, no additional prongs | Same
Edge fiducial |  $>10$ cm from active volume edge | Same
No dead region |  outside $z=700$ to 740 cm | Same
Containment Requirements: |  |
Edge containment |  closest approach $>15$ cm from edge | Same
Efficient region | N/A |  outside $y=(1/\sqrt{3})*z-117$
|  |  to $y=1/\sqrt{3}*z-80$ cm.
Particle Energy Requirements: |  |
Neutrino energy (range-based) |  $E_{\nu}<1200$ MeV | Same
Lepton kinetic energy |  $K_{\ell}>35$ MeV | Same
Proton kinetic energy |  $K_{p}>50$ MeV | Same
Analysis Orthogonality Cut |  Max Shower Frac $<0.2$ |  Max Shower Frac $>0.2$
Other Basic Quality Requirements: |  |
Opening Angle |  $>0.5$ radians | Same
Proton $\theta$ | $\cos(\theta_{p})>0$ | Same
Boostable |  $\gamma>1$, $0<|\beta|<1$ | Same
Shower 3-view energy consistency | N/A |  Inconsistency $<$ 200%
Vertex with BDT score that is |  lowest (i.e. best) |  highest (i.e. best)
BDT Requirements for Event Selection
BDT requirement |  $\nu_{\mu}$ BDT $<0.4$ |  $\nu_{e}$ BDT $>0.95$
PID Requirements for Event Selection
$\pi^{0}$ mass | N/A |  $\pi^{0}$ mass $<50$ MeV
muon MPID interaction score | N/A |  $<$ 0.2 if $E_{e}>100$ MeV
$\gamma/e$ image score ratio | N/A | $<2$
proton MPID interaction (max of 3 planes) |  $>0.9$ if $E_{\nu}<400$ MeV | N/A
Table 4.3: The specific used to define the $1\mu 1p$ and $1e1p$ samples. For
definitions of kinematic variables, see Table 4.1.

Figure 4.16: The $E_{\nu}^{\rm range}$ distribution for the $1e1p$ signal
sample, showing only the predicted event rate from the MC. The prediction from
the eLEE model is shown in the dashed blue line.

(a)

(b)

Figure 4.17: Figure 4.17(a) shows the distribution of the fractional error on
the neutrino energy for MC events in the $1e1p$ signal sample, restricted to
$200<E_{\nu}^{\rm Range}\;[{\rm MeV}]<1200$. Figure 4.17(b) shows the 2D
distribution of fractional error as a function of the true neutrino energy.

(a)

(b)

Figure 4.18: Figure 4.18(a) shows the post-vertex-identification efficiency of
true $\nu_{e}$ CCQE selection for subsequent stages of the $1e1p$ cuts. Figure
4.18(b) shows the true $\nu_{e}$ CCQE event rates over the full run 1-3
dataset after subsequent stages of the $1e1p$ cuts.

##  Chapter 5 The MicroBooNE Electron Neutrino Analysis: Results and
Discussion

This chapter presents the first results from the two-body CCQE analysis
described in chapter 4. They come from data taken during the first three
MicroBooNE run periods, corresponding to $6.67\times 10^{20}$ POT after data
quality cuts. The signal data were examined for the first time in Summer 2021
and presented to the community in October 2021 [148]. The results from the
two-body CCQE analysis were accompanied by results from the inclusive [170]
and MiniBooNE-like [171] $\nu_{e}$ analyses. All three MicroBooNE $\nu_{e}$
analyses are summarized in Ref. [116].

This chapter focuses on the statistical results from the $\nu_{e}$ two-body
CCQE analysis. We cover the background prediction, systematic error
evaluation, $\nu_{\mu}$ $1\mu 1p$ constraint procedure, and blinded analysis
approach in detail. We then discuss the series of statistical tests performed
using our $1e1p$ signal sample. We close with a discussion of the MicroBooNE
results, including a presentation of two non-MicroBooNE publications examining
the implications of the MicroBooNE data under two different explanations of
the MiniBooNE LEE: $3+1$ oscillations and excess $\overline{\nu}_{e}$
interactions.

Publications covered in this chapter for which I either held a leading role or
made major contributions: [148, 116, 102, 201]

###  5.1 First Results from the Two-Body CCQE Analysis

The two-body CCQE analysis observed 25 events across data from the first three
run periods. This is in good agreement with the total prediction in the
$200-1200$ MeV range, which is $27.4\pm 3.8_{\rm sys}\pm 5.2_{\rm stat}$
before the $1\mu 1p$ constraint discussed in section 5.1.3 [148]. One of these
events is shown in figure 5.1, including both pixel intensity images from the
TPC readout and SparseSSNet pixel-labeled images on all three planes. Event
displays for the remaining 24 selected events are included in Appendix B. The
$E_{\nu}^{\rm range}$ distribution of these events is shown in figure 5.2,
where the prediction is given both in terms of the underlying interaction type
and the final state event topology. The $1e1p$ BDT ensemble average score
distribution is shown in figure 5.3. Figure 5.4 shows the distributions of
four more interesting variables: $E_{e}$, $E_{p}$, $\theta_{e}$, and
$E_{\nu}^{QE-\ell}$. A full suite of data-prediction comparisons in 36
variables is provided in the Supplemental Material for Ref. [148], including
those described in table 4.1 and those used to train the BDT ensembles.

One can see clearly in figure 5.2 that no excess of $\nu_{e}$ $1e1p$ candidate
events is observed at low energy in the two-body CCQE analysis. Instead, we
observe a slight deficit compared to even the nominal prediction at low
energies, in conflict with the eLEE model represented by the blue dashed line.
This is a trend that is observed across all three $\nu_{e}$ analyses in
MicroBooNE [116]. The extent to which we agree or disagree with the eLEE model
will be discussed further in section 5.2.

One might also notice an apparent excess of data events in the 800-900 MeV bin
in figure 5.2. This appears to be a statistical fluctuation; other kinematic
variables do not exhibit such an excess of events, including the closely
related $E_{\nu}^{QE-\ell}$ distribution shown in figure 5.4(d). If this were
a systematic effect rather than a fluctuation, one would expect similar
behavior in the distributions of other variables. While these statements will
be quantified further in section 5.2, it is useful to include such a
qualitative discussion at this point.

The following subsections discuss the details behind our prediction in the
$E_{\nu}$ distribution. It is important to note that while we have already
shown the observed data in our signal channel, the procedures outlined in the
following subsections were frozen before examining this dataset. Specifically,
we followed the blinded analysis approach outlined in section 5.1.4 to avoid
biasing our selection.

Figure 5.1: Top: pixel intensity (color scale is in PIU as defined in section
4.4); Bottom: SparseSSNet labels; Left to Right: U, V, Y, planes. The white
circle indicates the reconstructed vertex. The horizontal axis corresponds to
the wire plane direction and the vertical axis corresponds to the electron
drift direction, which is measured using the arrival time of charge on the
wires.

(a)

(b)

Figure 5.2: The $1e1p$ sample $E_{\nu}$ distribution, comparing data (black
points) to the unconstrained prediction (stacked histogram) in the
$200<E_{\nu}<1200$ MeV region. The eLEE model prediction is represented by the
dashed blue line. The prediction is presented in terms of both interaction
type (figure 5.2(a)) and final state topology (figure 5.2(b)).  Figure 5.3:
Average $1e1p$ BDT ensemble score distribution comparing data to the
unconstrained prediction.

(a)

(b)

(c)

(d)

Figure 5.4: Comparison between data and unconstrained prediction in the
$E_{e}$ (figure 5.4(a)), $E_{p}$ (figure 5.4(b)), $\theta_{e}$ (figure
5.4(c)), and $E_{\nu}^{QE-\ell}$ (figure 5.4(d)) distributions of the $1e1p$
sample.

####  5.1.1 Background Estimation

In table 5.1 we give a breakdown of the MC events in the “background” category
of figures 5.3 and 5.4 by their interaction channel and event topology. The
“off-vertex” label in this table refers to events for which the reconstructed
neutrino vertex is further than 5 cm away from the true neutrino vertex. These
events can happen when a cosmic muon track happens to cross a photon shower
from a $\pi^{0}$ decay, or a muon decays to a Michel electron; in both cases,
the muon is mistaken for a proton. Table 5.1 indicates that $\nu_{\mu}$
resonant $\pi^{0}$ events are the most common background interaction channel,
followed by off-vertex events. One can also see that events with a $\pi^{0}$
in the final state are the dominant background topology.

There are not many MC events to inform the numbers reported in table 5.1. The
tools developed for this analysis produce a highly pure $\nu_{e}$ CCQE sample,
resulting in a low-statistics simulation sample for assessing the
non-$\nu_{e}$ background to the $1e1p$ signal. This is apparent in the
fluctuations of the background prediction in figures 5.3 and 5.4. In order to
obtain a more robust prediction of this background, we have elected to
leverage information on the energy distribution of non-$\nu_{e}$ background
events at a loose BDT score cutoff of 0.7 and extrapolate this to our signal
cutoff of 0.95. This is accomplished by fitting a parameterized Landau+linear
probability density function (PDF) to the background energy distribution at
the loose cutoff and scaling this prediction to the signal cutoff. As most of
these backgrounds come from $\nu_{\mu}$ interactions, this is referred to as
the $\nu_{\mu}$ background fit.

The Landau+linear shape is motivated empirically by the observation of a rise
in the background rate toward the lowest energies ($\lesssim 500$ MeV) and a
smaller rise toward the highest energies ($\gtrsim 800$ MeV), both for the
loose BDT score cutoff and the signal cutoff, as shown by the blue points in
figure 5.5. The Landau portion of the fit is additionally motivated by the
observation that a majority of $\nu_{\mu}$ background events contain a
$\pi^{0}$ in the final state, for which one of the decay photons is
misinterpreted as an electron shower. The reconstructed neutrino energy of
these events is governed predominately by the energy of this photon, which has
a tail extending to higher energies caused by pions with high momentum in the
lab frame. This tail is a characteristic feature of the Landau function. The
output of the fit and resulting error are used in place of the raw prediction
and statistical error for the $\nu_{\mu}$ backgrounds in this analysis. We
note here that the simulated $\nu_{\mu}$ events with a $\pi^{0}$ in the final
state are weighted according to the observation in our dedicated $\pi^{0}$
sample. The method for calculating these weights is described in Ref. [148].

The predicted background spectrum in each reconstructed neutrino energy bin
$f(E_{i})$ is generated by integrating the Landau+linear PDF $p(E)$ within
that bin. Specifically,

$\begin{split}&p(E)=\exp\bigg{[}\frac{-(E^{\prime}+e^{-E^{\prime}})}{2}\bigg{]}+aE,\\\
&E^{\prime}=(E-\mu)/\sigma,\\\ &f(E_{i})=\int_{E_{i}-\delta
E_{i}}^{E_{i}+\delta E_{i}}p(E)dE,\end{split}$ (5.1)

where $\mu$ and $\sigma$ are the center and width of the Moyal approximation
of the Landau function [202], $a$ is the linear slope parameter, $E_{i}$ is
the center of the i’th energy bin, and $\delta E_{i}$ is half of the bin
width. The Landau+linear fit is carried out using only shape information at a
loose BDT score cutoff of 0.7. In order to get the overall normalization, we
fit the BDT score distribution of the backgrounds to a linear PDF
$\tilde{p}(x)$, which we can integrate to get the total expected number of
background events $\tilde{f}(x)$ for a given BDT score cutoff $x$,

$\tilde{p}(x)=mx+b~{}~{}~{}~{}\tilde{f}(x)=\int_{x}^{1}\tilde{p}(x)dx,$ (5.2)

where $m$ and $b$ are the slope and bias parameters for the BDT score
distribution, respectively. The resulting shape+normalization fit for the
background distribution at the loose BDT cutoff of 0.7 and signal cutoff of
0.95 are shown in Fig. 5.5. One can see that the fit agrees with the raw MC
prediction within statistical error in both cases and that the error on the
fit is generally reduced compared to the simulation statistical error.

The errors on the fits in Fig. 5.5 are obtained by simulating pseudo-
experiments according to the covariance matrix of the fit parameters. This is
accomplished efficiently using Cholesky decomposition [203]. The overall
normalization error on the background rate at a given BDT score cutoff is
calculated using the error on the linear fit to the BDT score distribution of
background events. For the Landau+linear fit, we generate a shape-only
covariance matrix for the $\nu_{\mu}$ backgrounds in the nominal reconstructed
neutrino energy bins. The overall normalization error has been added as a
fully-correlated contribution to each bin of the shape-only $E_{\nu}$
covariance matrix. The uncorrelated uncertainty from the fit, given by the
square root of the diagonal entries of the full (shape $\oplus$ normalization)
covariance matrix, is indicated by the orange band in Fig. 5.5. This
covariance matrix replaces the nominal statistical uncertainty on the
simulated background, given by the quadrature sum of weights in each $E_{\nu}$
bin, which is indicated by the shaded blue band in Fig. 5.5.

The performance of the fit is evaluated on data by examining events slightly
below the signal region, with a BDT score in the range $[0.7,0.95]$. These
events are shown in Fig. 5.6. The raw MC prediction is given by the stacked
histogram, while the prediction incorporating the $\nu_{\mu}$ background fit
on top of the simulated $\nu_{e}$ prediction is given by the red line. The
data agree well with both predictions, indicating the consistency of the fit
method presented here with both the observed and MC-predicted background rate.

Interaction Channel | Predicted Rate
---|---
$\nu_{\mu}$ resonant $\pi^{0}$ | 1.26
$\nu_{\mu}$ resonant $\pi^{\pm}$ | 0.21
$\nu_{\mu}$ CCQE | 0.14
$\nu_{\mu}$ other | 0.19
Off-vertex | 0.93
Event Topology | Predicted Rate
$1\mu N\pi^{0}$ | 0.57
$0\mu N\pi^{0}$ | 1.09
$1\mu 1p$ | 0.14
Off-vertex | 0.93
Table 5.1: Breakdown of MC events in the “background” category of figures 5.3
and 5.4 over the range $200<E_{\nu}<1200$ MeV. The events are partitioned both
by the interaction channel and the event topology.

(a)

(b)

Figure 5.5: The fit to the $\nu_{\mu}$ background distribution to the $1e1p$
analysis. The shape fit is performed at a loose BDT score cutoff of 0.7
(figure 5.5(a)) and scaled to the signal cutoff of 0.95 (figure 5.5(b)). Blue
points represent the prediction from the simulation, with error bars
representing the Gaussian approximation of the statistical error (quadrature
sum of event weights). The orange line and corresponding shaded region
represent prediction and uncertainty, respectively, coming from the
Landau+linear fit. Figure 5.6: The data and MC prediction for events with a
$1e1p$ BDT score inside $[0.7,0.95]$. Good agreement is observed between data
and prediction. The prediction incorporating the Landau+linear background fit
is shown by the red line.

####  5.1.2 Evaluation of Systematic Uncertainties

There are five main sources of systematic uncertainty in this analysis,
related to the beam flux, neutrino-nucleus cross section modeling, hadron re-
interaction modeling, detector simulation, and finite statistics in the MC-
based prediction. These are handled within a covariance matrix formalism. For
each binned distribution, the diagonal entries of the covariance matrix
contain the variance in each bin, while the off-diagonal entries contain the
covariance between different bins. Thus, the covariance matrix is also
symmetric. The use of a covariance matrix is essential for encoding
correlations between binned distributions of selected events in the $1e1p$ and
$1\mu 1p$ samples, which is relevant for the constraint procedure described in
section 5.1.3. The total covariance matrix for a given distribution is simply
the element-wise sum of the covariance matrices for each source of systematic
uncertainty. A few definitions will be helpful for our discussion here; given
a covariance matrix $M_{ij}$ and binned prediction $N_{i}$, the fractional
covariance matrix is given by $F_{ij}\equiv M_{ij}/(N_{i}N_{j})$ and the
correlation matrix is given by $\rho_{ij}\equiv M_{ij}/\sqrt{M_{ii}M_{ij}}$.
The correlation matrix is meant to capture the degree of correlation between
two bins normalized within the range $[-1,1]$, where diagonal entries are 1 by
construction.

The flux, cross section, and re-interaction uncertainties are evaluated using
a reweighting method, in which event weights are modified according to tunable
parameters. Variations in the distributions of reconstructed variables are
used to determine the associated covariance matrix. The detector systematic
covariance matrix is evaluated using a dedicated set of simulation samples
that each alter the detector response in a different way. These alterations to
the detector response are derived using cosmic muon data; the full procedure
is described in detail in Ref. [204].

Flux uncertainties are related to hadron production in the target, secondary
hadron interactions, and focusing efficiency of the BNB magnetic horn.
Neutrino events are reweighted according to the species and momentum of the
parent meson [14, 205]. The largest flux uncertainty comes from $\pi^{+}$
production in the target. The flux uncertainties are also highly correlated
between the $1e1p$ and $1\mu 1p$ samples, as they both come predominately from
the same meson decay chain, $\pi^{+}\to\nu_{\mu}(\mu^{+}\to
e^{+}\nu_{e}\overline{\nu}_{\mu})$.

The cross section covariance matrix is calculated using both interaction-
specific reweightable parameters and FSI-related reweightable parameters in
GENIE v3.00.06. Thus, each MC event is reweighted according to its interaction
channel and any final state particles which underwent FSI. The full suite of
reweightable parameters is described in Ref. [180]. We also evaluate
uncertainties related to the difference between the $\nu_{e}$ CCQE and
$\nu_{\mu}$ CCQE cross sections [206]. The largest cross section uncertainties
come from CCQE-related parameters, which is expected given our signal
definition, and FSI-related parameters. FSI uncertainties are especially
important in this analysis, as FSI can cause $\nu_{e}$ CCQE events to deviate
from two-body scattering kinematics and also cause non-$\nu_{e}$ CCQE events
to look more consistent with two-body scattering. Both CCQE-related and FSI-
related effects impact the $1e1p$ and $1\mu 1p$ samples similarly and thus are
at least partially addressed by the constraint procedure in section 5.1.3.

Hadron re-interaction uncertainties describe the scattering of final state
protons and $\pi^{\pm}$ off of argon nuclei. The covariance matrix is
calculated by reweighting events according to Geant4 truth-level information
about the hadron trajectory [207]. These uncertainties are a small effect in
this analysis.

The detector systematic covariance matrix is calculated using a set of
dedicated simulation samples that vary certain aspects of the detector
response. These include the amplitude and width of TPC signals as a function
of the $x$ and $(y,z)$ positions as well the directions $\theta_{XZ}$ and
$\theta_{YZ}$ of the charged particle; the electron-ion recombination rate;
the electric field map in the TPC; and the light yield, attenuation, and
Rayleigh scattering length. As these samples are computationally expensive to
generate, we are not left with many MC events with which to construct a
covariance matrix. Thus, we rely on kernel density estimation [208, 209] to
smooth the predictions in each detector variation sample. This procedure does
not work for the $\nu_{\mu}$ backgrounds to the $1e1p$ analysis since we are
especially starved for MC statistics here, as discussed in section 5.1.1. We
instead assign a conservative uncorrelated 20% fractional detector systematic
uncertainty on these events, based on the extent of differences in the total
rate of selected background events in the detector variation samples.

The variance due to finite MC statistics is simply given by the quadrature sum
of event weights in each bin. These uncertainties are uncorrelated and give a
diagonal covariance matrix. The exception is finite MC statistics
uncertainties on the $\nu_{\mu}$ backgrounds in the $1e1p$ analysis, which are
calculated using the procedure described in section 5.1.1.

Finally, we note that to evaluate the reweightable systematic uncertainties on
the $\nu_{\mu}$ backgrounds to the $1e1p$ analysis, we use the selected MC
events to calculate a fractional covariance matrix, which is then scaled by
the parameterized background prediction derived in section 5.1.1. When
performing the constraint in section 5.1.3, we do not consider off-diagonal
systematic uncertainties between $\nu_{\mu}$ backgrounds in the $1e1p$ sample
and bins in the $1\mu 1p$ sample, as it is not straightforward to compute this
correlation.

The fractional uncertainty in each bin of the $E_{\nu}^{\rm range}$
distribution in the $1e1p$ and $1\mu 1p$ samples is shown in figure 5.7. These
correspond to the square root of the diagonal entries of the total systematic
covariance matrix calculated from all sources described above, divided by the
predicted event rate in each bin. In the $1e1p$ sample, flux, cross section,
and detector uncertainties dominate at low energies while detector
uncertainties dominate at high energies. The joint $1e1p$ and $1\mu 1p$
covariance and correlation matrices are shown in figure 5.8. The large off-
block-diagonal correlations between the $1e1p$ and $1\mu 1p$ samples will be
exploited in section 5.1.3.

(a)

(b)

Figure 5.7: The uncertainty in each bin of the $E_{\nu}^{\rm range}$
distribution of the $1e1p$ (figure 5.7(a)) and $1\mu 1p$ (figure 5.7(b))
samples.

(a)

(b)

Figure 5.8: The joint covariance (figure 5.8(a)) and correlation (figure
5.8(b)) matrices for the $E_{\nu}^{\rm range}$ distribution of the $1e1p$ and
$1\mu 1p$ samples.

####  5.1.3 Constraint from the $1\mu 1p$ Sample

As alluded to thus far, we use the observation in the $1\mu 1p$ control sample
to constrain the prediction and uncertainties in the $1e1p$ signal sample. We
rely on the conditional covariance method, which has also been used by the
MiniBooNE collaboration [76]. We perform the constraint in the $E_{\nu}^{\rm
range}$ distribution of each sample. Let us represent the joint covariance
matrix between both samples shown in figure 5.8(a) as

$M=\begin{pmatrix}M^{ee}&M^{e\mu}\\\ M^{\mu e}&M^{\mu\mu}\\\ \end{pmatrix}.$
(5.3)

Given a measurement $d_{\mu}$ and prediction $\mu_{\mu}$ in the $1\mu 1p$
channel, the constrained prediction $\mu_{e,{\rm const}}$ and covariance
matrix $M^{ee,{\rm const}}$ are given by [170]

$\begin{split}&\mu_{e,{\rm
const}}=\mu_{e}+M^{e\mu}\cdot(M^{\mu\mu})^{-1}\cdot(d_{\mu}-\mu_{\mu}),\\\
&M^{ee,{\rm const}}=M^{ee}-M^{e\mu}\cdot(M^{\mu\mu})^{-1}\cdot M^{\mu
e},\end{split}$ (5.4)

where $\mu_{e}$ is the nominal prediction in the $1e1p$ channel before the
constraint.

Figure 5.9 shows the comparison between data and prediction in the $1\mu 1p$
$E_{\nu}^{\rm range}$ distribution. This observation is used as input to the
constraint procedure in equation 5.4. The diagonal uncertainties in the $1e1p$
$E_{\nu}^{\rm range}$ distribution are reduced by up to a factor of $\sim 2$
across most bins, as shown in figure 5.10. In figure 5.11 we compare observed
data to the final constrained prediction and uncertainties in the $1e1p$
$E_{\nu}^{\rm range}$ distribution. As can be seen, the constraint procedure
slightly increases the $1e1p$ prediction at lower $E_{\nu}^{\rm range}$
values, though this is not a large effect. The result shown in figure 5.11
forms the basis of the statistical tests discussed in section 5.2.

Figure 5.9: The $E_{\nu}^{\rm range}$ distribution in the $1\mu 1p$ channel,
comparing data to the MC prediction. Figure 5.10: Fractional systematic
uncertainty in the $1e1p$ $E_{\nu}^{\rm range}$ distribution before and after
the $1\mu 1p$ constraint. Figure 5.11: Comparison between data and prediction
in the $1e1p$ $E_{\nu}^{\rm range}$ distribution after applying the $1\mu 1p$
constraint procedure.

####  5.1.4 Blinded Analysis Approach

Here we briefly cover the blind analysis procedure and related cross-checks
performed before examining data in the $1e1p$ signal sample. The main strategy
here was to use a sequential unblinding approach. Data were first examined for
events with a low BDT score ($<0.7$) or a high reconstructed neutrino energy
($>700$ MeV). Next, we examined data with a medium BDT score (within
$[0.7,0.95]$) or medium energy ($500-700$ MeV). Only after establishing good
agreement between data and prediction in all of these samples were potential
signal events first examined, where signal events are defined to have a $1e1p$
BDT score greater than 0.95 and reconstructed neutrino energy below 500 MeV.
This sequence was chosen because sensitivity to the eLEE model is minimized at
lower BDT scores and higher neutrino energies. Moving on to each subsequent
stage required signoff from the entire MicroBooNE collaboration. Figure 5.12
shows the comparison between data and prediction in the $1e1p$ BDT ensemble
average score distribution for this “medium BDT score” sample. Figure 5.6 is
another example of one of the “medium BDT score” plots.

We also explored a sample that made a number of cuts on kinematic variables in
order to mimic the $1e1p$ BDT ensemble cut. As humans are poorer optimizes
than BDTs, this sample contained a much larger rate of $\nu_{\mu}$ background
events compared to the official signal selection. Good agreement between data
and prediction was observed across all variables explored within this sample,
thus giving confidence in the MC description of backgrounds to the $1e1p$
signal sample.

Finally, before any $\nu_{e}$ analysis was able to look at data in their
signal sample, we each had to pass a series of fake data tests. In these
tests, five different data-like MC samples were generated with various unknown
alterations to the nominal MC. These included, for example, the presence of a
large $\nu_{e}$ eLEE signal and alterations to underlying cross section
models. All three $\nu_{e}$ analyses ran through the statistical results
presented in section 5.2 for each fake data sample. In all cases, the
statistical conclusions determined from each fake dataset were consistent with
the truth-level injected events. At this point, the $\nu_{e}$ analyses
(including the two-body CCQE analysis) were ready to move on to the real
signal dataset.

Figure 5.12: Comparison between data and prediction in the $1e1p$ BDT ensemble
average score distribution within the range $[0.7,0.95]$.

###  5.2 Statistical Interpretation

We perform three different statistical tests of the eLEE model: a $\chi^{2}$
goodness-of-fit test, a $\Delta\chi^{2}$ two-hypothesis test, and a signal
strength scaling test. We use the combined Neyman-Pearson (CNP) chi-square
test statistic [210],

$\chi^{2}_{\rm CNP}=(\vec{d}-\vec{\mu})^{T}(M^{\rm stat}_{\rm CNP}+M^{\rm
sys})^{-1}(\vec{d}-\vec{\mu}),$ (5.5)

where $\vec{d}$ and $\vec{\mu}$ are the observation and prediction for a
particular distribution, $M^{\rm sys}$ is the systematic covariance matrix
described in section 5.1.2, and

$(M^{\rm stat}_{\rm
CNP})_{ij}\equiv\frac{3\delta_{ij}}{\frac{1}{d_{i}}+\frac{2}{\mu_{i}}}.$ (5.6)

This formulation of the statistical error is meant to more closely approximate
the Poisson likelihood, which is important given the low-statistics nature of
the $1e1p$ analysis. Unlike the Poisson likelihood, the CNP chi-square test
statistic is compatible with the covariance matrix formalism developed for our
systematic uncertainty treatment.

All statistical tests use the constrained $1e1p$ prediction and uncertainties
shown in figure 5.11 when comparing to the 25 observed data events. We define
two distinct hypotheses which we will refer to in the following discussion:
$H_{0}$, the prediction without the eLEE model, considering only intrinsic
$\nu_{e}$ events and non-$\nu_{e}$ backgrounds (i.e. the stacked histogram in
figure 5.11); and $H_{1}$, the prediction including the eLEE model (i.e. the
dashed blue line in figure 5.11). Throughout this section, we report
probabilities, or $p$-values, calculated using the following frequentist
method. To test the distribution of a specific variable under a given
hypothesis, $10^{5}$ pseudo-experiments are generated by sampling the
multivariate distribution defined by the prediction and covariance matrix of
the hypothesis. This gives a predicted rate in each bin of the distribution;
we use this to sample from a Poisson distribution to get an integer number of
“observed” events in each bin. We calculate the desired test statistic for
each pseudo-experiment to build up the probability distribution of that test
statistic, which incorporates both systematic and statistical uncertainties.
Then, we can compare the observed value of the test statistic in the real
dataset to the calculated probability distribution to determine the
probability of observing a more extreme value–this defines the $p$-value.

####  5.2.1 Goodness of Fit

For the first test, we evaluate the goodness of fit in the $1e1p$ sample for
$H_{0}$ and $H_{1}$ using the $\chi^{2}_{\rm CNP}$ test statistic. This is
done both before and after the $1\mu 1p$ constraint. We also consider the
goodness-of-fit in both the full analysis range, $200<E_{\nu}^{\rm
range}\;[{\rm MeV}]<1200$, and the low energy region $200<E_{\nu}^{\rm
range}\;[{\rm MeV}]<500$, where the eLEE model prediction is highest. Table
5.2 shows the $\chi^{2}_{\rm CNP}$ per degree-of-freedom (DOF) for each of
these scenarios, as well as the corresponding $p$-values calculated via the
frequentist method described above.

The first thing to note there is significant tension with $H_{1}$, especially
after the $1\mu 1p$ constraint, as indicated by the low $p$-values. This is
true both over the full analysis energy range and in the low energy range.
After the constraint, there is greater than $3\sigma$ tension between $H_{1}$
and the observed data across the full energy range.

There is also disagreement between $H_{0}$ and the observed data, though not
to the same level as the $H_{1}$ case. Over the full analysis energy range,
the $p$-values in table 5.2 indicate $2.1\sigma$ ($2.4\sigma$) tension between
$H_{0}$ and the data before (after) the $1\mu 1p$ constraint. This seems to
come from two sources: the deficit in data compared to prediction below $500$
MeV in figure 5.11, and the excess in the 800-900 MeV bin in figure 5.11. The
first of these is possibly a real effect–similar low-energy deficits are
observed by all three $\nu_{e}$ analyses [116]. It is also observed in the
$E_{\nu}^{QE-\ell}$ distribution in figure 5.4(d) and is therefore not likely
to be related to proton misreconstruction. Our Michel electron sample covered
in section 4.4.3 indicates robust modeling of low-energy electron showers, so
this is also unlikely to be the culprit. The excess in the 800-900 MeV bin,
however, appears to be a statistical fluctuation. The other variables examined
in this analysis all show good agreement between data and MC and do not have
similar narrow excesses. For instance, the $E_{e}$, $E_{p}$, $\theta_{e}$, and
$E_{\nu}^{QE-\ell}$ distributions in figure 5.4 have a $\chi^{2}_{\rm
CNP}$/DOF ($p$-value) of 10.68/10 (0.42), 5.77/10 (0.80), 10.19/10 (0.48), and
11.46/10 (0.52), respectively. The $\chi^{2}_{\rm CNP}$/DOF and $p$-value of
all 36 variables examined in this analysis are reported in the Supplemental
Material of Ref. [148].

Nominal Predictions
---
Range | $H_{0}$ | $H_{1}$
$\chi^{2}_{\text{CNP}}/\text{DOF}$ |  $p$-value | $\chi^{2}_{\text{CNP}}/\text{DOF}$ |  $p$-value
200–500 MeV | 6.06/3 | 0.138 | 8.30/3 | 0.053
200–1200 MeV | 23.02/10 | 0.024 | 25.37/10 | 0.014
Constrained Predictions
Range | $H_{0}$ | $H_{1}$
$\chi^{2}_{\text{CNP}}/\text{DOF}$ |  $p$-value | $\chi^{2}_{\text{CNP}}/\text{DOF}$ |  $p$-value
200–500 MeV | 7.91/3 | 0.075 | 17.3/3 | 0.002
200–1200 MeV | 25.28/10 | 0.014 | 36.35/10 | $5.0\times 10^{-4}$
Table 5.2: Results from goodness-of-fit tests comparing observed $1e1p$ data
to the $H_{0}$ and $H_{1}$ predictions, reported via the
$\chi^{2}_{\text{CNP}}$ test statistic and the frequentist $p$-value. The top
half of the table considers the nominal prediction and uncertainties before
the $1\mu 1p$ constraint described in section 5.1.3, while the bottom half
considers the post-constraint prediction and uncertainties.

####  5.2.2 Two Hypothesis Test

The next test we perform is a two-hypothesis $\Delta\chi^{2}$ test between
$H_{0}$ and $H_{1}$. Specifically, we use the test statistic

$\Delta\chi^{2}=\chi^{2}_{{\rm CNP},H_{0}}-\chi^{2}_{{\rm CNP},H_{1}}.$ (5.7)

Thus, a positive (negative) value indicates better agreement with $H_{1}$
($H_{0}$). For this test, we consider the post-$1\mu 1p$-constraint prediction
and errors for the $1e1p$ sample. Figure 5.13 shows the distribution of this
test statistic under both hypotheses, generated using the frequentist method
described above. As expected, the $H_{0}$ distribution peaks below zero while
the $H_{1}$ distribution peaks above zero. These distributions can be used to
determine the one-sided exclusion sensitivity of each hypothesis in the two-
body CCQE analysis. If we denote the median $\Delta\chi^{2}$ value under
$H_{1}$ as $(\Delta\chi^{2})_{H_{1}}^{\rm med}$, the probability of $H_{0}$
generating a larger $\Delta\chi^{2}$ than $(\Delta\chi^{2})_{H_{1}}^{\rm med}$
is 0.003. This corresponds to an $H_{0}$ exclusion sensitivity of $2.7\sigma$.
Conversely, if we define $(\Delta\chi^{2})_{H_{0}}^{\rm med}$ in an analogous
way, the probability of $H_{1}$ generating a smaller $\Delta\chi^{2}$ than
$(\Delta\chi^{2})_{H_{0}}^{\rm med}$ is 0.017, corresponding to an $H_{1}$
exclusion sensitivity of $2.1\sigma$.

As shown in figure 5.13, the observed data result in a test statistic of
$\Delta\chi^{2}=-11.08$. The probability of observing a smaller
$\Delta\chi^{2}$ under $H_{0}$ and $H_{1}$ is 0.020 and $1.6\times 10^{-4}$,
respectively, corresponding to a Gaussian significance of $2.1\sigma$ and
$3.6\sigma$. Thus, the data exhibits mild tension with $H_{0}$ and significant
tension with $H_{1}$. A similar conclusion was reached for the goodness-of-fit
test; the tension with $H_{0}$ likely comes from the same sources discussed in
section 5.2.1. Due to the fact that non-negligible tension is observed with
respect to both hypotheses, we report our exclusion significance for $H_{1}$
using the CLs method [211, 212]. The CLs test statistic is defined as

${\rm CL}_{s}=\frac{p_{H_{1}}}{p_{H_{0}}},$ (5.8)

where ${\rm CL}_{s}=0.008$ for our observed data. Using this method, we reject
$H_{1}$ in favor of $H_{0}$ with a Gaussian significance of $2.4\sigma$.

Figure 5.13: Distributions of the $\Delta\chi^{2}$ test statistic defined in
equation 5.7 for $H_{0}$ (red) and $H_{1}$ (blue), calculated by generating
$10^{5}$ pseudo-experiments under each hypothesis. The $\Delta\chi^{2}$ value
of the data is also shown.

####  5.2.3 Signal Strength Scaling Test

The final statistical test we perform on our $1e1p$ dataset is a one-parameter
signal strength scaling test. Specifically, we multiply by the eLEE weights in
figure 4.1(b) by an overall normalization factor $x_{\rm LEE}$ which can vary
within $[0,\infty]$. In this construction, $x_{\rm LEE}=0$ corresponds to
$H_{0}$ and $x_{\rm LEE}=1$ corresponds to $H_{1}$. We can use the Feldman-
Cousins procedure [25] along with the test statistic

$\Delta\chi^{2}(x_{\rm LEE})=\chi^{2}_{\rm CNP}(x_{\rm LEE})-\min_{x_{\rm
LEE}\in[0,\infty]}\\{\chi^{2}_{\rm CNP}(x_{\rm LEE})\\}$ (5.9)

to determine confidence intervals on $x_{\rm LEE}$. We can use the Asimov
dataset, in which the observation is exactly the predicted value in each bin,
to determine the sensitivity to each hypothesis in this test. Under $H_{0}$,
the expected upper bound on $x_{\rm LEE}$ at the 90% ($2\sigma$) confidence
level is 0.75 (0.98). Under $H_{1}$, the expected confidence interval at the
$1\sigma$ ($2\sigma$) confidence level is $[0.53,1.66]$ ($[0.28,2.67]$), and
the expected sensitivity to rule out $x_{\rm LEE}=0$ is $2.8\sigma$.

Figure 5.14 shows the confidence intervals on $x_{\rm LEE}$ derived from the
observed data. The best-fit value is $x_{\rm LEE}=0$, which is expected given
the deficit in data at the lowest neutrino energies where the eLEE model
weights are the largest. The 90% ($2\sigma$) upper bound on $x_{\rm LEE}$ is
0.25 (0.38). These upper bounds are stronger than the expected confidence
intervals under $H_{0}$, which is mainly due to the aforementioned deficit at
low neutrino energies.

Figure 5.14: Confidence intervals on $x_{\rm LEE}$ calculating using the
Feldman-Cousins procedure. The solid and dotted lines indicate the confidence
level with which a given $x_{\rm LEE}$ is disfavored, calculated using the
Feldman-Cousins method [25] and Wilks theorem [26], respectively. The
MiniBooNE statistical and systematic errors are shown as a band around $x_{\rm
LEE}=1$.

###  5.3 Discussion and Outlook

The observation in the $1e1p$ signal channel of the two-body CCQE analysis is
inconsistent with a MiniBooNE-like excess of $\nu_{e}$ CCQE interactions.
Under a two-hypothesis test, we rule out the nominal eLEE model ($H_{1}$) at
the $2.4\sigma$ using the CLs method. We can set an upper bound on the signal
strength scaling parameter $x_{\rm LEE}$ of 0.38 at the $2\sigma$ level.

This is consistent with the results from all three MicroBooNE $\nu_{e}$
analyses [116]. As shown in figure 5.15(a), three of the four signal channels
from the three $\nu_{e}$ analyses observe a deficit of selected events at the
lowest neutrino energies, where the eLEE model is supposed to peak. The
exception is the $1e0p0\pi$ channel from the MiniBooNE-like sample; however,
this is the least sensitive of the four channels, as it is dominated by
non-$\nu_{e}$ backgrounds. Figure 5.15(b) shows the results from the signal
strength scaling test in all four channels. With the exception of the
$1e0p0\pi$ channel, all analyses rule out $x_{\rm LEE}=1$ at greater than the
$2\sigma$ confidence level. The most sensitive analyses observe an upper bound
$x_{\rm LEE}\lesssim 0.5$ at the $2\sigma$ confidence level.

The first MicroBooNE results also included a search for $\Delta\to N\gamma$
decays as an explanation for the MiniBooNE excess. The analysis did not
observe an excess of $\Delta\to N\gamma$ candidate events and ruled out the
required enhancement to explain the MiniBooNE anomaly at the 95% confidence
level. The full details of this analysis can be found in Ref. [80].

Thus, the MicroBooNE results disfavor both an enhancement of $\nu_{e}$ CC
interactions and $\Delta\to N\gamma$ decays as the sole explanation of the
MiniBooNE excess. Potential explanations of the MiniBooNE anomaly must now
remain consistent with these MicroBooNE results. One promising avenue involves
dark sector models which produce a photon or $e^{+}e^{-}$ pair in the
MiniBooNE detector. Such event signatures have not yet been explored by
MicroBooNE, though this work is in progress. Crefch:neutrissimos explores a
photon-based explanation of MiniBooNE, involving the decay of a heavy neutral
lepton to a single photon via a transition magnetic moment. The theoretical
and experimental neutrino physics communities are actively working to develop
innovative models to explain the MiniBooNE LEE and equally innovative
measurements to test those models.

The following subsections close with two projects I co-led with various
collaborators to follow up on the state of the MiniBooNE excess after the
MicroBooNE results. The first examines the $3+1$ model simultaneously in
MiniBooNE and MicroBooNE [102]. The second explores the implications of a
$\overline{\nu}_{e}$-based MiniBooNE LEE for the MicroBooNE results [201].

(a)

(b)

Figure 5.15: Figure 5.15(a) shows the observation compared to the nominal
($H_{0}$) prediction in all four signal channels from the three MicroBooNE
$\nu_{e}$ analyses, including statistical errors on the data points and
systematic errors on the prediction. The eLEE prediction ($H_{1}$) is also
indicated by the red line. Figure 5.15(b) shows the observed $1\sigma$ and
$2\sigma$ confidence intervals on $x_{\rm LEE}$ from all four signal channels.
The $2\sigma$ expected sensitivity of each channel is shown in red.

####  5.3.1 Publication: MiniBooNE and MicroBooNE Combined Fit to a $3+1$
Sterile Neutrino Scenario

The initial results from the MicroBooNE $\nu_{e}$ analyses did not explicitly
test short-baseline $\nu_{\mu}\to\nu_{e}$ oscillations from the $3+1$ model.
Shortly after the MicroBooNE results came out, two studies appeared assessing
the sensitivity of the MicroBooNE $\nu_{e}$ analyses to the $3+1$ model [188,
213]. Following up on this, the MiniBooNE collaboration released the first
combined fit to a $3+1$ model using MicroBooNE and MiniBooNE data [102]. This
study was led by myself and Austin Schneider. The full Physical Review Letters
publication is included below. We considered the two-body CCQE and inclusive
MicroBooNE $\nu_{e}$ analyses, as these placed the most stringent constraints
on the eLEE model and provided data releases with the requisite detail to
perform a $3+1$ fit [116]. The most important result from this study is Figure
2, which indicates that there are still allowed regions in $3+1$ parameter
space at the $3\sigma$ confidence level in the MiniBooNE-MicroBooNE combined
fit. This was also the first MiniBooNE analysis to account for all three
oscillation probabilities in equation 1.36: $\nu_{e}$ appearance, $\nu_{\mu}$
disappearance, and $\nu_{e}$ disappearance. Previous $3+1$ fits to MiniBooNE
data have only considered $\nu_{e}$ appearance, as this is the dominant effect
in the electron-like channel. However, $\nu_{e}$ and $\nu_{\mu}$ disappearance
can have a non-negligible impact on the intrinsic $\nu_{e}$ and misidentified
$\nu_{\mu}$ background rate in the MiniBooNE electron-like prediction. The
$\nu_{e}$ disappearance probability is also important when considering
MicroBooNE data, as the prediction in each MicroBooNE $\nu_{e}$ analysis is
dominated by intrinsic $\nu_{e}$ interactions [148, 170]. MicroBooNE has since
released official constraints on the $3+1$ model using data from the inclusive
$\nu_{e}$ analysis [214]–these constraints are consistent with our
implementation.

See pages - of templates/Publications/PhysRevLett.129.201801.pdf See pages -
of templates/Publications/MBuB_supplemental.pdf

####  5.3.2 Publication: Implications of MicroBooNE’s low sensitivity to
electron antineutrino interactions in the search for the MiniBooNE excess

The first MicroBooNE results indeed set strong limits on $\nu_{e}$ CC
interactions as an explanation for the entire MiniBooNE LEE. However, as
pointed out in Ref. [201], this conclusion changes if the MiniBooNE LEE comes
instead from $\overline{\nu}_{e}$ CC interactions. This is driven mainly by
differences in the low-energy suppression of $\overline{\nu}_{e}$-nucleus and
$\nu_{e}$-nucleus cross sections, coming from interference between the vector
and axial form factors as well as nucleon binding energy effects. Binding
energy effects are especially important in MicroBooNE–within a non-isoscalar
nucleus like 40Ar, it takes more energy for a $\overline{\nu}_{e}$ to turn a
proton into a neutron than it does for a $\nu_{e}$ to turn a neutron into a
proton. This is not necessarily the case in MiniBooNE, as 12C is an isoscalar
nucleus. This study was led by myself and Matheus Hostert. The full Physical
Review D publication can be found below. The most important results are shown
in Figure 5 and Figure 6, which indicate that the MicroBooNE sensitivity to
the MiniBooNE LEE decreases as more of the excess is attributed to
$\overline{\nu}_{e}$ rather than $\nu_{e}$ CC interactions. Specifically, the
MicroBooNE data are consistent at the $2\sigma$ confidence level with a
scenario in which the entire MiniBooNE LEE comes from $\overline{\nu}_{e}$
interactions.

See pages - of templates/Publications/PhysRevD.107.092002.pdf

##  Chapter 6 Neutrissimos: Heavy Neutral Leptons with a Dipole Moment

In this chapter, we discuss a phenomenological explanation of the MiniBooNE
excess involving both an eV-scale sterile neutrino and an MeV-scale heavy
neutral lepton (HNL) with an effective transition magnetic moment coupling to
active neutrinos. This “mixed model” was introduced in Ref. [27] and refined
in Ref. [31]. The HNL in this model is hereafter referred to as a
“neutrissimo” [31], a term which has previously been used to describe
interacting HNLs in the literature [215, 216].

We begin with a brief overview of the theoretical framework behind
neutrissimos. We then motivate the mixed model and discuss its implications
for global neutrino experiments. A robust evaluation of the neutrissimo signal
in MiniBooNE is performed. For a neutrissimo mass $m_{\mathcal{N}}\sim
500~{}{\rm MeV}$, single photons from visible neutrissimo decay are shown to
provide a reasonable explanation of the energy and angular distributions of
the MiniBooNE excess. World-leading constraints on neutrissimos are derived
using elastic scattering data from the MINER$\nu$A experiment [217, 218, 219].
While the MINER$\nu$A data exclude a large region of parameter space for
$m_{\mathcal{N}}=\mathcal{O}(100~{}{\rm MeV})$, they do not rule out the
MiniBooNE solution at the 95% confidence level. The details behind this
analysis are given in Ref. [31]; the full published manuscript is included in
section 6.4.

Publications covered in this chapter for which I held a leading role: [27, 31]

###  6.1 Dipole-Portal Neutrissimos

The neutrissimo $\mathcal{N}$ is a right-handed neutrino that couples to the
active neutrinos through a transition magnetic moment. This coupling is
described by the dimension-five dipole operator

$\mathcal{L}_{D}=d_{\alpha\mathcal{N}}\,\overline{\nu_{\alpha}}\sigma_{\mu\nu}F^{\mu\nu}\mathcal{N}_{R}+\text{
h.c.},$ (6.1)

where $d_{\alpha\mathcal{N}}$ is the dipole coupling between the right-handed
neutrissimo $\mathcal{N}_{R}$ and the left-handed weak flavor eigenstate
$\nu_{\alpha}$, and $F^{\mu\nu}=\partial^{\mu}A^{\nu}-\partial^{\nu}A^{\mu}$
is the EM field strength tensor. The dipole coupling has units of inverse
energy, as equation 6.1 is an effective operator. This means it is agnostic to
the ultraviolet (UV) completion: the short-distance physics generating the
dipole coupling at a higher energy scale. Possible UV completions for the
dipole operator include (but are not limited to) $SU(2)_{L}\times SU(2)_{R}$
theories [220], leptoquark models [81], and models introducing additional
Higgs multiplets or other charged scalars [221, 222, 223, 224]. We will treat
neutrissimos within the effective field theory framework throughout this
chapter. Figure 6.1 shows a Feynman diagram of the effective operator in
equation 6.1. It facilitates a number of new interactions involving the
neutrissimo, the most relevant of which are shown in figure 6.2.

Equation 6.1 is only valid below the electroweak scale, as it does not respect
$SU(2)_{L}\times U(1)_{Y}$ gauge invariance. Before electroweak symmetry
breaking (EWSB), the dipole coupling is described by the gauge invariant
dimension six operator

$\mathcal{L}\supset\frac{1}{\Lambda^{2}}\overline{L_{\alpha}}\widetilde{H}\sigma^{\mu\nu}\mathcal{N}_{R}\left(C_{B}^{\alpha}\,B_{\mu\nu}+C_{W}^{\alpha}\,W^{a}_{\mu\nu}\sigma_{a}\right)+\text{
h.c.},$ (6.2)

where $\Lambda$ is the UV completion scale, $L_{\alpha}$ is the left-handed
lepton $SU(2)_{L}$ doublet, $\widetilde{H}=i\sigma_{2}H^{*}$ is the conjugate
Higgs field, and $C_{B}^{\alpha}$ and $C_{W}^{\alpha}$ are the Wilson
coefficients for couplings to the $B_{\mu}$ and $W_{\mu}$ vector gauge bosons,
respectively. The dipole coupling in equation 6.1 after EWSB is then given by

$d_{\alpha\mathcal{N}}=\frac{v_{h}}{\sqrt{2}\Lambda^{2}}(C_{B}^{\alpha}\cos\theta_{W}+C_{W}^{\alpha}\sin\theta_{W}),$
(6.3)

where $v_{h}$ is the Higgs vacuum expectation value and $\theta_{W}$ is the
Weinberg angle. Note that couplings are also generated with the $W$ and $Z$
bosons after EWSB, but interactions involving these couplings are suppressed
by $G_{F}$ and are thus not as important as the photon coupling.

Transition magnetic moments between different neutrino species have been
discussed extensively in the literature [225, 226, 227, 228, 229, 81, 130,
223, 230]. The imprint of heavy neutral leptons with transition magnetic
moments in experimental data has received particular attention [130, 81, 231,
123, 124, 126, 127, 232, 233, 234, 235, 236, 135, 237, 238, 239, 240]. One
potential complication of this model is that large transition magnetic moments
often lead to large Dirac masses $m_{\nu\mathcal{N}}$ for neutrinos [81, 130,
31]. This is because the photon line can typically be removed from figure 6.1,
resulting in a contribution to the Dirac mass term $\mathcal{L}\supset
m_{\nu\mathcal{N}}\overline{\nu}_{L}\mathcal{N}_{R}$. Dipole coupling
strengths of interest to MiniBooNE ($d\sim 10^{-6}$ GeV) would generate Dirac
mass contributions $m_{\nu\mathcal{N}}=\mathcal{O({\rm MeV}})$ [81].
Considering a neutrissimo mass $m_{\mathcal{N}}\sim 500$ MeV, the typical
seesaw relation of equation 1.23 would predict active neutrino masses
$m_{\nu}>>1~{}{\rm eV}$, in conflict with existing limits. This can be
remedied by the inverse seesaw mechanism [241, 242], in which the smallness of
active neutrino masses comes from approximate lepton number conservation
without relying on large right-handed neutrino masses. [81, 31].

Another complication is the potential existence of mass mixing between the
neutrissimo and the active neutrinos. This mixing is heuristically given by
$U_{\alpha\mathcal{N}}\sim m_{\nu\mathcal{N}}/m_{\mathcal{N}}$ and is large
even in the inverse seesaw mechanism [81]. If this operator is not suppressed,
the interactions of the neutrissimo are dominated by those involving the $W$
and $Z$ bosons. Any potential mass mixing would be subject to constraints from
accelerator neutrino experiments, which can be strong for
$m_{\mathcal{N}}=\mathcal{O}(100~{}{\rm MeV})$ [243]. Specifically, searches
for HNLs produced in kaon decay at T2K (using the ND280 near detector) [244],
E949 [245], and PS191 [246] set upper bounds on the muon-flavor mixing of
$|U_{\mu\mathcal{N}}|^{2}\leq 10^{-8}-10^{-9}$ up to the kaon mass threshold
of $m_{\mathcal{N}}\sim 400~{}{\rm MeV}$. Within the range $400\leq
m_{\mathcal{N}}~{}[{\rm MeV}]\leq 2000$, NuTeV is able to set upper bounds of
$|U_{\mu\mathcal{N}}|^{2}\leq 10^{-6}-10^{-7}$ by searching for HNLs produced
in $D$ meson decays.

In order to suppress mass mixing, the Dirac mass itself must be small compared
to the neutrissimo mass. One possible method to accomplish this is the
approximate horizontal $SU(2)_{H}$ symmetry proposed by Voloshin [247]. The
idea of an approximate $SU(2)_{H}$ symmetry initially gained popularity as a
means to enable the $\nu_{e}$ diagonal/transition magnetic moment solution to
the solar neutrino problem [248, 249, 250]. For the heavy neutral lepton case,
we consider a scenario in which $(\nu_{L}^{C},\mathcal{N}_{R})$ transforms as
a doublet under $SU(2)_{H}$ [251, 81, 224]. The antisymmetric nature of the
dipole operator
$\overline{\nu}_{\alpha}\sigma_{\mu\nu}F^{\mu\nu}\mathcal{N}_{R}$ in flavor
space means that it transforms like a singlet under $SU(2)_{H}$, while the
symmetric nature of the Dirac mass term makes it a triplet under $SU(2)_{H}$.
Thus, approximate $SU(2)_{H}$ conservation will suppress the Dirac mass
contribution from the photon-less version of figure 6.1 while allowing for
large contributions to the transition magnetic moment. Though it is not clear
the level of fine-tuning required within this model to accommodate the
MiniBooNE solution [224], we consider this sufficient motivation to consider a
scenario in which the mass mixing operator of the heavy neutrino is negligible
compared to the dipole operator. It is also worth pointing out the leptoquark
UV completion of the dipole operator proposed in Ref. [81], which can
naturally suppress the mass mixing of neutrissimos via an $SU(2)_{H}$ symmetry
without the need for significant fine-tuning.

$\nu_{L}$$\mathcal{N}_{R}$$\gamma$ Figure 6.1:  Feynman diagram depicting the
effective dipole operator of equation 6.1.

(a)

(b)

(c)

Figure 6.2: New interactions involving the neutrissimo that are enabled by the
dipole operator in equation 6.1, including three-body $\pi^{0}$ decay (figure
6.2(a)), Primakoff-like upscattering (figure 6.2(b)), and neutrissimo decay
(figure 6.2(c)).

###  6.2 Overview of the Mixed Model

The mixed model consists of the eV-scale sterile neutrino described in section
1.4 and the MeV-scale neutrissimo described above in section 6.1. This model
is motivated by the significant tension observed between appearance and
disappearance experiments in $3+1$ global fits [93, 12, 11, 252], which comes
from strong signals indicating $\nu_{e}$ appearance and $\nu_{e}$
disappearance coupled with the lack of such a signal in the $\nu_{\mu}$
disappearance channel. This tension is measured using the parameter goodness-
of-fit (PG) test [253], which uses the test statistic

$\begin{split}&\chi^{2}_{PG}=\chi^{2}_{\rm global}-(\chi^{2}_{\rm
app}+\chi^{2}_{\rm dis}),\\\ &N_{PG}=(N_{\rm app}+N_{\rm dis})-N_{\rm
global},\end{split}$ (6.4)

where $\chi^{2}_{\rm global}$, $\chi^{2}_{\rm app}$, and $\chi^{2}_{\rm dis}$
are the total $\chi^{2}$ at the $3+1$ best-fit parameters for global,
appearance-only, and disappearance-only experiments, respectively. The number
of degrees-of-freedom $N$ is defined similarly. The probability that the
appearance and disappearance experiments come from the same underlying $3+1$
model can be calculated using the $\chi^{2}_{PG}$ distribution for $N_{PG}$
degrees-of-freedom.

This tension appears to be driven significantly by MiniBooNE data [27, 12].
Within the Harvard-Columbia-MIT global fit that existed at the time of Ref.
[27], the PG test tension decreased from $4.8\sigma$ to $2.5\sigma$ upon
removing the MiniBooNE dataset. The allowed regions in $\Delta
m_{41}^{2}$-$\sin^{2}2\theta_{\mu e}$ parameter space for the global,
appearance-only, and disappearance-only experiments both before and after
removing MiniBooNE are shown in figure 6.3 and figure 6.4, respectively. One
can see visually that there is more overlap in the $\Delta m^{2}$ dimension
between the three different allowed regions after removing MiniBooNE. This is
likely because the MiniBooNE excess peaks strongly at low energies, an effect
which is difficult to produce in the $3+1$ model, as discussed in section 2.2.
The latest MiniBooNE $3+1$ fits discussed in section 5.3.1 specifically show
that the oscillation hypothesis under-predicts the excess in the lowest bins
of the MiniBooNE electron-like dataset. This behavior drives the global fit to
the $\Delta m_{41}^{2}<1~{}{\rm eV}^{2}$ region, in contrast to the higher
values ($\Delta m_{41}^{2}\gtrsim 1~{}{\rm eV}^{2}$) preferred by the reactor
and gallium anomalies.

The global fit without MiniBooNE, shown in the leftmost panel of figure 6.4,
strongly prefers a solution at $\Delta m_{41}^{2}\sim 1~{}{\rm eV}^{2}$ and
$10^{-4}\lesssim\sin^{2}2\theta_{\mu e}\lesssim 10^{-3}$. The best fit point
is specifically at $\Delta m_{41}^{2}=1.3~{}{\rm eV}^{2}$ and
$\sin^{2}2\theta_{\mu e}=6.9\times 10^{-4}$. We take these as our benchmark
$3+1$ parameters, which fix the $\nu_{\mu}\to\nu_{e}$ oscillation contribution
to the MiniBooNE electron-like channel. The oscillation prediction alone is
not sufficient to explain the MiniBooNE LEE; it can only accommodate $\sim
10\%$ of the total excess events and peaks at $E_{\nu}\sim 600~{}{\rm MeV}$,
inconsistent with the low-energy nature of the MiniBooNE excess. However, the
benefit of these specific oscillation parameters is that they can explain
LSND, reactor, and gallium data while remaining consistent with null results
from $\nu_{\mu}$ disappearance searches. The impact of different mixing angles
within the allowed region of the global fit without MiniBooNE is discussed in
Ref. [27], and the impact of the (larger) oscillation contribution predicted
by the MiniBooNE-MicroBooNE combined best fit [102] is discussed in Ref. [31].

The remaining majority of the MiniBooNE excess must then come from another
source. In the mixed model, we ascribe the bulk of the LEE to the visible
decays of neutrissimos. This model works well for MiniBooNE, as neutrissimos
can be produced abundantly in the BNB through the first two diagrams in figure
6.2: Primakoff upscattering $\nu A\to\mathcal{N}A$ and the three-body decay
$\pi^{0}\to\gamma\nu\mathcal{N}$. The third diagram, neutrissimo decay
$\mathcal{N}\to\nu\gamma$, could then provide a source of extra photons in
MiniBooNE. As discussed in chapter 2, these photons would show up in the
electron-like sample due to MiniBooNE’s reliance on Cherenkov rings for
particle identification. Models of heavy neutrino decay to single photons have
been previously considered as a solution to the MiniBooNE anomaly [123, 254,
124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135]. We extend upon
these studies by performing a robust analysis of the neutrissimo model in
MiniBooNE within the context of the mixed model. This analysis is described in
detail in the next section. It turns out that compared to the $3+1$ model,
neutrissimo decays are more consistent with the energy and angular
distributions of the MiniBooNE excess. Thus there are two distinct advantages
of the mixed model: (1) it can provide an explanation for other anomalies in
the neutrino sector while relieving tension in global $3+1$ fits, and (2) it
can provide a better explanation of the MiniBooNE excess itself.

Figure 6.3: $3+1$ global fits including MiniBooNE, considering global (left),
appearance-only (middle), and disappearance-only (right) experiments. The
allowed regions in $3+1$ parameter space at the 90%, 95%, and 99% confidence
levels are shown by the red, green, and blue points, respectively. The best-
fit point is indicated by the star.

Figure 6.4: $3+1$ global fits without MiniBooNE, considering global (left),
appearance-only (middle), and disappearance-only (right) experiments. The
allowed regions in $3+1$ parameter space at the 90%, 95%, and 99% confidence
levels are shown by the red, green, and blue points, respectively. The best-
fit point is indicated by the star.

###  6.3 Neutrissimos in MiniBooNE

The analysis presented in Ref. [27] used a custom simulation to estimate the
single photon event rate from neutrissimo decays in the MiniBooNE detector.
The updated MiniBooNE analysis in Ref. [31] relied instead on the public
LeptonInjector 111github.com/Harvard-Neutrino/LeptonInjector software package,
which will be described in more detail in section 6.3.1. In Ref. [27] we
considered neutrissimo production via both $\pi^{0}\to\gamma\nu\mathcal{N}$
decays and $\nu A\to\mathcal{N}A$ Primakoff upscattering. The former process
was estimated using the Sanford-Wang parameterization of the $\pi^{0}$
production rate at the BNB ] [110, 14] and was found to be negligible for the
dipole parameters of interest to MiniBooNE. We therefore explicitly ignored
this process, considering Primakoff upscattering as the sole neutrissimo
production mechanism.

The $\nu_{\alpha}A\to\mathcal{N}A$ upscattering rate was calculated using the
differential cross section with respect to the Mandelstam variable
$t=Q^{2}=(p_{\mathcal{N}}-p_{\nu})^{2}$ [27, 130, 255, 232],

$\begin{split}\frac{d\sigma_{\nu_{\alpha}A\to\mathcal{N}A}}{dt}=\frac{2\alpha
d_{\alpha}^{2}}{M}\bigg{[}&F_{1}^{2}(t)\Big{(}\frac{1}{E_{r}}-\frac{1}{E_{\nu}}+m_{\mathcal{N}}^{2}\frac{E_{r}-2E_{\nu}-M}{4E_{\nu}^{2}E_{r}M}+m_{\mathcal{N}}^{4}\frac{E_{r}-M}{8E_{\nu}^{2}E_{r}^{2}M^{2}}\Big{)}\\\
+&\frac{F_{2}^{2}(t)}{4M^{2}}\Big{(}\frac{2M}{E_{\nu}^{2}}((2E_{\nu}-E_{r})^{2}-2E_{r}M)+m_{\mathcal{N}}^{2}\frac{E_{r}-4E_{\nu}}{E_{\nu}^{2}}+\frac{m_{\mathcal{N}}^{4}}{E_{\nu}^{2}E_{r}}\Big{)}\bigg{]},\end{split}$
(6.5)

where $\alpha$ is the fine structure constant, $d$ is the dipole coupling,
$E_{\nu}$ is the SM neutrino energy, $m_{\mathcal{N}}$ is the mass of the
heavy neutrino, $M$ is the target mass, $t=-(p_{\mathcal{N}}-p_{\nu})^{2}$ is
the momentum transfer, $E_{r}=-t/2M$ is the target recoil energy, and
$F_{1/2}(t)$ are the electric/magnetic target form factors, respectively. Note
that the term proportional to $E_{r}m_{N}^{4}$ in the $F_{1}$ line only exists
for spinless nuclei, and must be replaced for nonzero spin nuclei [127]. In
the case of coherent scattering off of a nucleus, $F_{1}$ receives a $Z^{2}$
enhancement and is therefore dominant over $F_{2}$. Thus, we explicitly
neglect the $F_{2}$ term for this analysis. In Ref. [27], $F_{1}$ was
calculated using a dipole parameterization. However, since the dipole
approximation tends to underestimate the cross section at high squared
momentum transfer $Q^{2}$, Ref. [31] updated to the Fourier-Bessel
parameterization,

$F^{\rm
FB}(Q)=Z^{2}N\frac{\sin(QR)}{QR}\sum_{n}\frac{(-1)^{n}a_{n}}{n^{2}\pi^{2}-Q^{2}},$
(6.6)

where $Z$ is the proton number of the struck nucleus, $N$ is a normalization
factor ensuring $F^{\rm FB}(0)=1$, and $a_{n}$ are coefficients obtained from
nuclear scattering data [256, 257, 258, 259]. For the case of inelastic
upscattering off of individual nuclei, we use the nucleon electric and
magnetic form factors discussed in the appendix of Ref. [27]. Inelastic
upscattering exhibits only a linear enhancement in the number of protons and
neutrons, which suppresses its effect compared to the $Z^{2}$ enhancement of
the coherent upscattering cross section. Ref. [31] also includes contributions