Measurement of neutral current single $\pi^0$ production on argon with the MicroBooNE detector

We report the first measurement of $\pi^0$ production in neutral current (NC) interactions on argon with average neutrino energy of $\lesssim1$~GeV. We use data from the MicroBooNE detector's 85-tonne active volume liquid argon time projection chamber situated in Fermilab's Booster Neutrino Beam and exposed to $5.89\times10^{20}$ protons on target for this measurement. Measurements of NC $\pi^0$ events are reported for two exclusive event topologies without charged pions. Those include a topology with two photons from the decay of the $\pi^0$ and one proton and a topology with two photons and zero protons. Flux-averaged cross-sections for each exclusive topology and for their semi-inclusive combination are extracted (efficiency-correcting for two-plus proton final states), and the results are compared to predictions from the \textsc{genie}, \textsc{neut}, and \textsc{NuWro} neutrino event generators. We measure cross sections of $1.243\pm0.185$ (syst) $\pm0.076$ (stat), $0.444\pm0.098\pm0.047$, and $0.624\pm0.131\pm0.075$ $[10^{-38}\textrm{cm}^2/\textrm{Ar}]$ for the semi-inclusive NC$\pi^0$, exclusive NC$\pi^0$+1p, and exclusive NC$\pi^0$+0p processes, respectively.


I. INTRODUCTION
Neutrino-nucleus cross sections have been the subject of intense study both experimentally and within the theory community in recent years due to their role in interpreting neutrino oscillation measurements and searches for other rare processes in neutrino scattering [1]. While neutrino oscillation experiments primarily rely on measuring the rate of charged current (CC) interactions, it is also important that we build a solid understanding of inclusive and exclusive neutral current (NC) neutrino interactions.
NC neutrino interactions are of particular importance to ν e andν e measurements in the energy range of a few hundred MeV. This is especially true for detectors that cannot perfectly differentiate between photon-and electron-induced electromagnetic showers, and therefore where NC π 0 production followed by the subsequent decay π 0 → γγ can be misidentified as ν e orν e CC scattering. Misidentification of photons as electrons complicates the interpretation of ν e appearance measurements aiming to measure subtle signals. These include sterile neutrino oscillation searches with the upcoming Short Baseline Neutrino (SBN) experimental program [2] and CP violation measurements and mass hierarchy determination with the future Deep Underground Neutrino Experiment (DUNE) [3].
Furthermore, NC π 0 events can contribute as background to searches for rare neutrino scattering processes such as NC Δ resonance production followed by Δ radiative decay, or NC coherent single-photon production at energies below 1 GeV [4]. This is primarily a consequence of the limited geometric acceptance of some detectors, whereby one of the photons from a π 0 decay can escape the active volume of the detector. Depending on a detector's ability to resolve electromagnetic shower substructure, NC π 0 events can further contribute as background to searches for new physics beyond the Standard Model (BSM), such as e þ e − production predicted by a number of BSM models [5][6][7][8][9].
Finally, NC measurements themselves can provide a unique channel for probing new physics. For example, searches for nonunitarity in the three-neutrino paradigm or searches for active to sterile neutrino oscillations are possible via NC rate disappearance measurements [10,11].
Such searches can provide complementary information to nonunitarity or light sterile neutrino oscillation parameters otherwise accessible only through CC measurements.
Using a liquid argon time projection chamber (LArTPC) as its active detector, MicroBooNE [12] shares the same technology and neutrino target nucleus as the upcoming SBN and future DUNE experiments. MicroBooNE's 85 metric ton active volume LArTPC is situated 468.5 m away from the proton beam target in the muon-neutrinodominated Booster Neutrino Beam (BNB) at Fermilab [13] which is also used by SBN. The resulting neutrino beam has a mean energy hE ν i ¼ 0.8 GeV and is composed of 93.7% ν μ , 5.8%ν μ , and 0.5% ν e =ν e . MicroBooNE's cross-section measurements on argon are therefore timely and directly relevant to these future (SBN and DUNE) programs.
We present the first measurement of neutrino-induced NC single-π 0 (1π 0 ) production on argon with a mean neutrino energy in the 1 GeV regime, which is also the higheststatistics measurement of this interaction channel on argon to date. This measurement is relevant to the physics programs of experiments that operate in the few-GeV regime (SBN [2], DUNE [3], NOνA [14,15], T2K [16], and Hyper-K [17]), especially those which share argon as a target material. Additionally, this measurement has been used to provide an indirect constraint to the rate of NC 1π 0 backgrounds in MicroBooNE's recent search for a singlephoton excess [4]. The only previous results for NC 1π 0 scattering on argon are from the ArgoNeuT Collaboration using the NuMI beam which has a much higher mean neutrino beam energy of 9.6 GeV for ν μ and of 3.6 GeV forν μ [18].
The interaction final states that are measured in this analysis are defined as where A represents the struck (argon) nucleus, A 0 represents the residual nucleus, and X represents exactly one or zero protons plus any number of neutrons, but no other hadrons or leptons. The protons are identifiable in the MicroBooNE LArTPC by their distinct ionizing tracks while the π 0 is identifiable through the presence of two distinct electromagnetic showers, one for each photon from the π 0 → γγ decay, with kinematic properties such that they reconstruct to approximately the π 0 invariant mass. These one proton and zero proton samples are used first to perform a rate validation check and subsequently in three distinct cross-section measurements. By leveraging the capability of LArTPCs to detect and identify protons we perform the world's first exclusive NC π 0 þ 0p and NC π 0 þ 1p cross-section extractions and additionally measure the cross section for NC π 0 interactions semi-inclusively using both the one proton and zero proton samples combined. Each of these cross-section extractions utilizes a distinct signal definition. The signal definitions for the two exclusive measurements place a threshold on true proton kinetic energy of greater than 50 MeV, while the semiinclusive measurement allows for any number of protons. The signal definitions for all three measurements also require that there are no other hadrons or leptons in the final state (as noted above). MeV-scale photons, which may arise from nuclear de-excitation processes within the struck nucleus, are allowed in the final state. Finally, the signal definitions allow for interactions of all flavors of neutrinos that are present: ν μ ,ν μ , ν e , andν e .
These definitions are comparable to other historical NC π 0 measurements which typically require one and only one π 0 meson and little hadronic activity in the detector [19][20][21][22][23][24][25][26][27][28]. This differs from the more inclusive approach of the ArgoNeuT experiment motivated both by its higher energy beam as well as the need to mitigate the low statistics of its data sample [18]. Making use of the MicroBooNE LArTPC's power in examining hadronic final state multiplicities and kinematic properties with high resolution, the flux-averaged cross sections extracted in this analysis extend our understanding of this important interaction channel. The exclusive cross sections reported provide new information useful for the tuning of NC 1π 0 production and nuclear final-state interactions (FSIs) in neutrino-argon scattering models, while the semi-inclusive cross section enables a direct comparison to the MiniBooNE measurement of NC π 0 production.

II. ANALYSIS OVERVIEW
This measurement uses data corresponding to a BNB exposure of 5.89 × 10 20 protons on target (POT), collected during the period 2016-2018 and referred to as "Runs 1-3" in many of the subsequent figures. Neutrino-argon interactions are simulated using a custom tune [29] of the GENIE neutrino event generator v3.0.6 [30,31] (based on model set G18_10a_02_11a) adopted by the MicroBooNE Collaboration. This tune specifically targets CC quasielastic (QE) and CC multinucleon interaction models and overall has very little direct effect on this NC-focused analysis. GENIE v3 uses the Berger-Sehgal [32,33] model for resonant production of π 0 and includes improved agreement with an expanded dataset for the A dependence of FSIs, updated form factors [34], updated diagrams for pion production processes [33,35,36], and a new tune to neutrino-proton and neutrino-deuterium cross-section data [31]. The MicroBooNE Monte Carlo (MC) prediction further makes use of GEANT4 v4_10_3_03c [37] for particle propagation and reinteractions within the detector and a custom detector response model all implemented within the LArSoft framework [38].
The MicroBooNE data and MC reconstruction chain begins by reading out and processing the ionization charge signals detected on the 8,192 wires that make up the three anode planes of the MicroBooNE LArTPC. The procedure includes noise removal [39] and signal processing as described in [40] and [41]. Localized regions of interest referred to as "hits" are then identified and fit to Gaussian pulses. The collection of these hits and their characteristics such as readout time, wire channel number, and integrated charge are then used as input to the Pandora pattern recognition framework for further processing [42]. The Pandora framework clusters and matches hits across three 2D projected views of the MicroBooNE active TPC volume to form 3D reconstructed objects. These objects are then classified as tracklike or showerlike based on a multivariate classifier score and aggregated into candidate neutrino interactions. Pandora also reconstructs a candidate neutrino interaction vertex based on the position and orientation of the reconstructed tracks and showers which represents the most likely position of the neutrino interaction.
Being a surface detector, MicroBooNE is subject to a constant stream of high-energy cosmic rays impinging on the detector that substantially outnumber the neutrino interactions and form the largest background to candidate neutrino interactions. To incorporate the effect of cosmicray contamination in the simulation, cosmic-ray data recorded in situ at MicroBooNE, when the beam is not present, are used as overlays (at the wire signal waveform level) to simulated neutrino interactions. During the 2.3 ms that it takes to "drift" ionization charge associated with neutrino interaction final states across the maximum 2.56 m drift distance, Oð10Þ cosmic rays are expected to enter the detector. In order to reduce this cosmic-ray contamination, scintillation light recorded by the MicroBooNE photodetector system is matched to candidate neutrino interactions during reconstruction and is also required to occur in time with the 1.6 μs long BNB neutrino spill.
To select a high-purity sample of BNB neutrino NC 1π 0 interactions, a series of topological, preselection, and boosted decision tree (BDT)-based selections are applied. This results in two mutually exclusive final selection topologies: 2γ1p, which targets two photons and one proton in the final state, and 2γ0p, which targets two photons and zero protons in the final state. The different selection stages are described below, along with the details of the systematic uncertainty evaluation.

A. Topological selection and preselection
The event selection begins with topology-based criteria for candidate neutrino interactions identified by Pandora and targets two mutually exclusive topological definitions: (a) two showers and one track (2γ1p), and (b) two showers and zero tracks (2γ0p). The two showers correspond to the photons expected from π 0 decay. The presence of a track corresponds to a reconstructed proton exiting the nucleus while the zero-track case suggests either a low-energy proton that is not reconstructed or no charged hadrons at all exiting the nucleus.
Once events with the desired signal topologies are identified, a series of loose "preselection" requirements is applied to reduce obvious backgrounds or misreconstructed events. These preselection requirements include shower energy thresholds of 30 MeV for the leading shower and 20 MeV for the subleading shower in both topologies. The preselection also requires that the reconstructed neutrino interaction point be contained in a fiducial volume, defined as at least 5 cm away from any TPC wall, in order to help reduce the number of selected events with tracks that exit the detector. For the 2γ1p topology, the nonzero conversion distance of photons is explicitly used by requiring that each shower has a reconstructed start point of at least 1 cm from the reconstructed neutrino interaction vertex. Typically the reconstructed neutrino interaction vertex is identified as the start of the reconstructed proton candidate track. In order to remove a very small number of poorly reconstructed events in which the candidate track is not consistent with the hypothesis of originating from the candidate neutrino interaction vertex, a requirement is placed to ensure the track start point is always within 10 cm of the reconstructed neutrino interaction vertex. The efficiency of selecting NC 1π 0 þ 0 (1)p events using these preselection requirements is 21.5% (19.9%). Note that the efficiency of the 1p selection is lower because of the additional requirements placed on the track reconstruction.

B. Boosted decision tree-based selection
After applying the preselection requirements, the remaining signal and background are further differentiated and separated using two tailored BDTs trained on simulation. The gradient boosting algorithm XGBoost [43] is used to train each of the BDTs. They take as input various reconstructed kinematic, geometric, and calorimetric variables both for the signal (defined as an NC interaction with identically one π 0 in the final state) and for the background interactions. Because the two tailored BDTs target different topologies, notably including one track in the case of 2γ1p and zero tracks in the case of 2γ0p, the signal definitions used for the two BDTs are slightly different. NC π 0 events with exactly one proton with true kinetic energy above 20 MeV are used as the training signal for the 2γ1p BDT while NC π 0 events with no protons with true kinetic energy above 20 MeV are used as the training signal for the 2γ0p BDT. We note that the 20 MeV threshold used in the BDT training is lower than the 50 MeV proton kinetic energy threshold used later during cross-section extraction, as during training we are aiming to push the threshold as low as possible. Each BDT is trained on ten reconstructed variables. Due to the existence of a proton candidate track in the 2γ1p sample, these ten variables differ for each BDT. They are listed below.
Variables used in both 2γ1p and 2γ0p BDTs: (i) Leading and subleading shower impact parameters: The perpendicular distance between the back projection of the reconstructed shower and the candidate neutrino interaction point which is a metric of how well each shower "points" back to the reconstructed neutrino interaction point. (ii) Leading and subleading shower conversion distances: Defined as the distance between the reconstructed start of the shower and reconstructed neutrino interaction point. (i) Reconstructed energy of the subleading shower.
(ii) Leading and subleading shower geometric length per unit energy: The ratio of each shower's geometric length to its reconstructed energy. The geometric length is an estimate of the 3D extent of the electromagnetic shower. (iii) Pandora "neutrino score": A multivariate classifier in the Pandora reconstruction suite which scores all reconstructed neutrino candidates based on their geometric and kinematic features as to how likely a candidate is due to a neutrino interaction or cosmic in origin. (iv) Reconstructed leading shower vertical angle: Direction in the vertical plane with respect to the beam axis. By construction, BDT scores lie on the interval of [0, 1]. After training, the resulting BDT score distributions, tested on a statistically independent simulation and dataset, are shown in Fig. 1. The simulation and data points agree across the full range of BDT classifier score within systematic and statistical uncertainties (the definition of these systematic uncertainties is described in detail in Sec. II C). The bimodal distribution of the 2γ1p BDT response indicates greater separation power between signal and background compared to that for 2γ0p because the addition of the reconstructed track gives access to an entirely separate handle on background rejection. For this and subsequent MC simulation comparisons to data, the simulation predictions are broken down into the following eight categories, based on GENIE truth-level information: (i) NC 1π 0 : All neutral current interactions that produce one exiting π 0 regardless of incoming neutrino flavor. This is our targeted signal selection, and it is further split into two sub-categories, "NC 1π 0 Coherent" and "NC 1π 0 Noncoherent" contributions, based on their interaction types. Noncoherent scattering occurs when a neutrino interacts with a nucleon inside the argon nucleus, potentially knocking out one or more nucleons. In coherent scattering the neutrino interacts with the nucleus as a whole, leaving it in its ground state. This interactions occurs with low momentum transfer, and as such the resulting π 0 tends to be very forward relative to the incoming neutrino beam. (ii) NC Δ → Nγ: Leading Standard Model source of NC single-photon production below 1 GeV originating from radiative decay of the Δð1232Þ baryon. (iii) CC ν μ 1π 0 : All ν μ CC interactions that have one true exiting π 0 . (iv) CC ν e =ν e Intrinsic: All CC ν e orν e interactions regardless of whether or not a π 0 was emitted. (v) BNB Other: All remaining BNB neutrino interactions that take place in the active TPC volume of MicroBooNE and are not covered by the above five categories, such as multiple π 0 events and η meson decay. See Sec. II F for more details. (vi) Dirt (Outside TPC): All BNB neutrino interactions that take place outside the MicroBooNE active TPC but have final states that enter and interact inside the active TPC detector. This can originate from scattering off liquid argon in the cryostat vessel outside the active TPC volume or from interactions in the concrete and "dirt" surrounding the cryostat itself. (vii) Cosmic Data: Coincident cosmic-ray interactions that take place during a BNB spill but without any true neutrino interaction present. The final NC 1π 0 -enriched samples are selected by placing a requirement on the BDT score distribution that maximizes the product of NC 1π 0 signal efficiency and purity. This corresponds to a threshold on the BDT scores of > 0.854 and > 0.950 for 2γ1p and 2γ0p, respectively. The final distributions are provided and discussed in Sec. II E.

C. Systematic uncertainty evaluation
Systematic uncertainties on the MC simulation prediction include contributions from uncertainties in the neutrino flux, the cross-section modeling, hadron reinteractions, detector effects, and the effect of finite statistics used in the background predictions (both simulations and cosmic-ray data).
The flux systematic uncertainties incorporate hadron production uncertainties where the Booster proton beam hits the beryllium target, uncertainties on pion and nucleon scattering in the target and surrounding aluminum magnetic focusing horn of the BNB, and mismodeling of the horn current. Following [44], these are implemented by reweighting the flux prediction according to neutrino type, parentage, and energy, and studying the propagated effects on the final event distributions.
The cross-section uncertainties incorporate modeling uncertainties on the GENIE prediction [29][30][31], evaluated by GENIE reweighting tools. The default GENIE uncertainties on NC resonant production arising from NC resonant vector and axial mass parameters of m V ¼ 0.840 AE 0.084 GeV and m A ¼ 1.120 AE 0.224 GeV, respectively, were assumed. GENIE uses an effective cascade empirical model for hadronic final-state interactions, called hA2018, which allows for reweighting to estimate the effect on final distributions. For more information on cross-section uncertainties in MicroBooNE, please see Ref. [29].
The hadron-argon reinteraction uncertainties are associated with the propagation of hadrons through the detector, as modeled in GEANT4 [37]. Both charged pions and proton reinteractions during propagation were considered and their impact estimated using the GEANT4REWEIGHT tool [45].
The detector modeling and response uncertainties are evaluated using MicroBooNE's novel data-driven technique for assessing and propagating LArTPC detectorrelated systematic uncertainties [46]. This approach uses in situ measurements of distortions in the TPC wire readout waveform signals-caused by detector effects such as electron diffusion, electron drift lifetime, electric field, and the electronics response-to parametrize these effects at the TPC wire level. This provides a detector modelagnostic way to study and evaluate their effects on the high level variables and, subsequently, the final event distributions. Additional detector systematics corresponding to variations in the charge recombination model, the scintillation light yield, and space charge effects [47,48] are separately evaluated and also included.

D. Shower energy calibration
Electromagnetic shower reconstruction in LArTPCs is known to be a lossy process primarily due to misclustering and thresholding effects. Current reconstruction algorithms often miss small, low-energy hits in an electromagnetic shower when clustering objects, and some of the hits that are reconstructed may fall below the energy threshold. On average, these effects are expected to yield shower energy losses of approximately 20% [49]. This can be seen in Fig. 2, which shows true and reconstructed shower energy for a dedicated high statistics sample of simulated true NC 1π 0 events, where the reconstructed shower energy falls systematically below the true shower energy in simulation. By performing a linear fit to the most probable values of reconstructed shower energy in bins of true shower energy, shown as the pink straight line in Fig. 2, a correction factor is extracted which brings the reconstructed values closer to expectation. This fit results in an energy correction that is applied to all reconstructed showers,  and represents a correction of approximately 20%, as expected.

E. Final selected spectra
After applying the BDT requirements, 1130 selected data events remain with 634 and 496 falling into the 2γ1p and 2γ0p selections, respectively. For the 2γ1p selection, the BDT score requirement efficiency is 85.6% and the purity is 63.5% while for the 2γ0p selection, the efficiency and purity are 58.8% and 52.9%, respectively. The 2γ1p and 2γ0p BDT selection efficiencies and purities are both calculated relative to their signal definition, with proton multiplicity counted with a kinetic energy threshold of 50 MeV. The efficiencies at each stage of the analysis are provided in Table I, and the total efficiency for each selection is shown as a function of (a) true π 0 momentum and (b) true proton kinetic energy in Fig. 3. Overall, the 1p selection is more efficient and of higher signal purity relative to the 0p selection due to the existence of a reconstructed particle track which greatly helps to tag the neutrino interaction point and reject backgrounds. This track information, particularly track calorimetry, provides an additional handle on the neutrino interaction mode; a protonlike track is highly indicative of an NC 1π 0 interaction whereas CC interactions generally have a muon track in the final state.
The resulting distributions as a function of the reconstructed two-photon invariant mass are shown in Fig. 4. The invariant mass is reconstructed from the energy and direction of the two photon candidate showers as where cos θ γγ is the opening angle between the two showers. For the 2γ1p case where a track has been identified as a candidate proton, the directions of the showers and thus the opening angle between them are calculated by constructing the direction between the candidate neutrino interaction point and the shower start point. For the 2γ0p selection, however, no such candidate track exists. Instead, the shower direction and opening angle are entirely estimated from the geometric shape of the showers themselves. A Gaussian plus linear fit is performed to each observed distribution in data to extract the reconstructed π 0 invariant mass while taking into account the non-π 0 background contamination. For the 2γ1p event sample, this fit gives a Gaussian mean of 138.9 AE 2.1 MeV=c 2 with a width of 31.7 AE 2.4 MeV=c 2 . For the 2γ0p event sample, the corresponding fit gives a Gaussian mean of 143.3 AE 3.2 MeV=c 2 with a width of 47.9 AE 4.9 MeV=c 2 . As a goodness-of-fit test, the resulting χ 2 per degree of freedom is 1.20 and 1.45 for the 2γ1p and 2γ0p fits, respectively.   Fig. 6. This is defined as the angle between the lab-frame π 0 momentum direction and the decay axis of the two daughter photons in the CM frame, This quantity should be an isotropic flat distribution for true π 0 → γγ signal events, and any deviation from this can highlight regions of inefficiency in reconstruction or selection. As can be seen in Fig. 6, for both 2γ1p and 2γ0p selections, the distributions taper off at high cos θ CM which corresponds to increasingly asymmetric π 0 decays. When reconstructing asymmetric π 0 decay events, it is more likely that the subleading photon shower is missed due to its low energy. Note, however, that the observed data show the same trend as the simulation within uncertainty. Figure 7 additionally highlights the reconstructed photon conversion distance for all showers in the final 2γ1p selection. Well-reconstructed showers with conversion distances as far as 100 cm from the candidate neutrino interaction are observed. This helps validate the assumption that the reconstructed showers are indeed likely to be true photons as Oð100Þ MeV photons are expected to have a mean free path in argon of ≈20 cm. Note that, as the 2γ0p selection does not have any visible hadronic activity for tagging the interaction point, the corresponding conversion distance is significantly harder to estimate.
Finally, Fig. 8 shows two example event displays of selected events in data for both the 2γ1p and 2γ0p topologies. Each event shows two well-reconstructed showers pointing back to a common interaction point with properties consistent with those being photons from a π 0 decay.

F. Background discussion and validation
In order to validate the background modeling in this analysis we developed a background rich sideband selection by inverting the BDT score cuts shown in Fig. 1. This gives us a high statistics sample of "CC1π 0 " and "BNB Other" background categories with which to compare to data. In addition to inverting the BDT score an additional cosmic rejection cut, where we require the Pandora neutrino score to be > 0.5, is applied to provide a higher purity of the backgrounds we wish to study. These distributions are shown in Figs. 9(a) and 9(b) for 2γ1p and 2γ0p respectively. This is particularly useful for 2γ1p inverted selection as the resulting spectrum is rich in CC1π 0 for higher reconstructed π 0 momenta, and richer in BNB Other at low momenta. For the 2γ0p background rich sample there is still a significant amount of rejected NC π 0 signal events, but the enhanced backgrounds still provide additional validation. The data is observed to be in good agreement with the prediction, within assigned uncertainties, with a χ 2 per degree of freedom (d.o.f.) of 24.61=20 and 17.49=22 for 2γ1p and 2γ0p respectively. This gives us confidence that the backgrounds and uncertainties are sufficiently modeled for a cross-section extraction to proceed.
We can also break down the BNB Other category further in order to improve our understanding of this important background. We find that approximately 75% of BNB Other events contain true photons reconstructed, and in case of 2γ1p 84% have true protons reconstructed. This indicates that despite being a background category the BDTs are indeed selecting events with a very high purity of true photons and protons, as is the target of the selection. The approximate breakdown of BNB Other after applying the BDT cuts is found to be (i) ≈25% are events in which multiple π 0 are exiting the nucleus but only 1π 0 is reconstructed correctly, (ii) ≈25% are events in which no π 0 exits the nucleus but due to baryon or charged pion rescattering in the argon, a π 0 is subsequently created and reconstructed, (iii) ≈25% are events in which there is no π 0 and a NC proton is reconstructed as the track with cosmic contamination resulting in a 2γ event mimicking a π 0 , (iv) ≈20% are events containing a NC η → γγ decay event. These are generally rejected due to them having higher energies, but a small number are selected as they are topologically identical to the π 0 → γγ decay signal, (v) ≈5% Miscellaneous other reconstruction failures, representing less than 1% of total background events. In the case of CC1π 0 we see a similar situation, with 97.3% (98.1%) of 2γ1p (2γ0p) background events having at least 1 shower matched to a true π 0 and with 78% of tracks in the 2γ1p sample being correctly matched to a proton. As such, the vast majority of these events have the correct target particle content, but rather the muon itself is missed. This primarily occurs when the muon is incorrectly clustered into a photon electromagnetic shower due to close proximity, or when the muon is correctly reconstructed but is misidentified as a cosmic muon, with the associated π 0 then being reconstructed as an isolated neutrino event.
III. NC π 0 RATE VALIDATION NC 1π 0 events contribute as a dominant background to NC single-photon production measurements carried out or FIG. 5. The reconstructed π 0 momentum ((a) and (b)) and reconstructed π 0 angle with respect to the neutrino beam ((c) and (d)) for both the 2γ1p ((a) and (c)) and 2γ0p ((b) and (d)) final selected data. The prediction shows agreement with the observed data within assigned uncertainties for the ranges shown, although a systematic deficit is observed in the total event rates as is discussed in Sec. III.
planned by MicroBooNE such as searches for NC Δ radiative decay [4], NC coherent single-photon production, or more rare e þ e − pair production motivated in BSM theories. In addition to using these selected events as a calibration sample for understanding and validating shower reconstruction performance, they are also used to validate the observed overall rate of this process as currently modeled with GENIE. Assuming GENIE provides a sufficient description of the observed data, this sample can and has been used to provide an in situ constraint on NC 1π 0 misidentified backgrounds, e.g. as in [4]. Alternatively, these measurements can be used to increase our understanding of our current NC π 0 modeling and potentially motivate GENIE tuning.
As shown in Fig. 5, both the 2γ1p and the 2γ0p selections see an overall deficit in data relative to the MC prediction. This is more pronounced in the 2γ1p selection where the ratio of the number of selected data events to the number of selected simulated events is 0.79. As it is also known that the GENIE branching fraction of coherent NC 1π 0 production on argon is significantly lower than expectation extrapolated from MiniBooNE's π 0 measurement on mineral oil [23], the possibility of a correction to GENIE predictions on both noncoherent and coherent NC 1π 0 production is explicitly examined. The MC predictions are fitted to data allowing both coherent and noncoherent NC 1π 0 rates to vary. Both normalization-only and normalization plus shape variations to the coherent and noncoherent rates are explored; all yield similar conclusions. This section describes the normalization plus shape variation fit in detail.
The normalization plus shape variation fit is performed as a function of reconstructed π 0 momentum for both 2γ1p and 2γ0p selections, using [0, 0.075, 0.15, 0.225, 0.3, 0.375, 0.45, 0.525, 0.6, 0.675, and 0.9] GeV/c bin limits. In the fit, MC predicted coherent NC 1π 0 events are scaled by a normalization factor N coh , and MC predicted noncoherent NC 1π 0 events are scaled on an event-by-event basis depending on their corresponding true π 0 momentum according to ða þ bj ⃗p true π 0 jÞ, where the true π 0 momentum is given in [GeV/c]. This linear scaling as a function of π 0 momentum was chosen because it was the simplest implementation that was consistent with the observed data-to-MC deficit, as observed in Fig. 5(a).
At each set of fitting parameters (N coh , a, b), a χ 2 is evaluated between the scaled prediction for this parameter set and the observed data using the Combined-Neyman-Pearson χ 2 [51]. The χ 2 calculation makes use of a covariance matrix including statistical and systematic uncertainties and correlations corresponding to the scaled prediction. Flux, cross section, detector and GEANT4 systematic uncertainties are included in the fit including bin-tobin systematic correlations. As the goal of the fit is to extract the normalization and scaling parameters of the coherent and noncoherent NC 1π 0 rates, the cross-section normalization uncertainties of coherent and noncoherent NC 1π 0 are not included. Note that the cross-section normalization uncertainties of coherent and noncoherent NC 1π 0 are only removed for the purposes of this fit and not for the crosssection extraction described in the following section.
The data-extracted best-fit parameters correspond to a ¼ 0.98 and b ¼ −1.0 [c/GeV] for the scaling parameters of the noncoherent NC 1π 0 events, and N coh ¼ 2.6 for the NC coherent π 0 normalization factor with no enhancement, N coh ¼ 1, being allowed within the 1σ error bands. This best-fit gives a χ 2 =d:o:f of 8.46=17. The χ 2 =d:o:f at the GENIE central value (CV) prediction is 13.74=20 yielding a Δχ 2 between the GENIE CV and the best-fit point of 5.28 for 3 d.o.f. Although the goodness-of-fit χ 2 =d:o:f: values for both scenarios are acceptable due to the generally large uncertainties, the momentum-dependent shift is preferred over the GENIE CV at the 1.43σ level. The 1D marginalized Δχ 2 distributions in Fig. 10 also confirm that the GENIE CV FIG. 9. Background rich sidebands for validation created by inverting the BDT selection cuts, for (a) 2γ1p and (b) 2γ0p samples. We observe good agreement between the data and prediction, within the assigned uncertainties giving us confidence in the background modeling.
FIG. 10. The distribution of marginalized Δχ 2 , as a function of flat normalization factor for (a) coherent NC 1 π 0 momentum-independent scaling factor, (b) noncoherent NC 1π 0 momentum-independent scaling factor, and (c) coefficient of momentum-dependent scaling factor for NC noncoherent 1π 0 , marginalized over the other two parameters. The red arrows indicate parameter values expected for the GENIE central value prediction. The 1σ, 90% and 99% confidence level lines are based on the assumption that the distribution follows a χ 2 distribution with 1 degree of freedom. prediction agrees with data within uncertainty. The data and MC comparisons of the reconstructed π 0 momentum distributions scaled to the best-fit parameters are provided in Fig. 11 and, compared to those corresponding to the GENIE CV, show better agreement with data after the fit.
While the data suggest that GENIE may overestimate NC 1π 0 production, the results demonstrate that the GENIE prediction of NC 1π 0 s is accurate within uncertainty. This validates the approach of using the measured NC 1π 0 event rate as a powerful in situ constraint of GENIE-predicted NC 1π 0 backgrounds as in [4]. We stress that while this result motivates a momentum-dependent shift in how we model our NC π 0 events, we do not apply this change to our modeling, relying instead on the fact that the assigned uncertainty covers the observed discrepancy. Future work will investigate further the possibility of tuning GENIE with results such as these to obtain better model predictions. On the other hand, it is natural to extract a data-driven NC 1π 0 cross section on argon using these selections, and compare to a number of neutrino event generators, including GENIE. This is described below.

A. Methodology
The prescription for calculating the cross section is provided in Eq. (5) where the components are defined as follows: N obs NC1π 0 , N cosmic , and N bkg denote the number of selected data events, the number of background events arising from cosmic rays traversing the detector, and the number of expected beam-correlated background events, respectively; ϵ NC1π 0 denotes the efficiency of selecting NC 1π 0 events; Φ denotes the integrated flux; and N targets denotes the number of argon atoms in the fiducial volume of the analysis.
This calculation is performed independently using each of the 2γ1p and 2γ0p selections to measure an exclusive cross section. These measurements are denoted as the NC π 0 þ 1p and NC π 0 þ 0p cross sections, respectively; in each case one or zero protons is explicitly required in the signal definition (described in detail below). Additionally, the calculation is performed using the combined 2γð0 þ 1Þp selection to measure a semi-inclusive cross section, NC π 0 , with no requirement on the number of protons in the signal definition. Note that this semi-inclusive measurement is efficiency corrected to include 2þ proton final states that are not included in the final selected events (11% of the total number of NC π 0 interactions, per GENIE). As noted in Sec. II C, the simulation is run multiple times to encompass the effect of varying underlying sources of systematic uncertainty. The calculation of each cross section is performed separately in each of these systematic "universes" to guarantee that all correlations between components of the cross section are handled correctly. This is done using tools from the MINERvA Analysis Toolkit [52].
Both selections, as well as their combination, correspond to approximately, but not identically, the same POT, provided in Table II (due to differences in the computational processing of the two samples). To extract the semiinclusive cross section from the combined 2γð0 þ 1Þp selection, the relevant 2γ0p distributions are scaled down by the ratio between the POT of the 2γ1p data sample (smaller POT) and the POT of the 2γ0p data sample and then are added to the 2γ1p distributions. This operation is performed for N obs NC1π 0 , N cosmic , N bkg , and the numerator of the efficiency.
N obs NC1π 0 and N cosmic are measured in data and therefore there is no systematic uncertainty attributed to them. These values are reported in Table II. N bkg is extracted from the simulation, and we note that many of the key backgrounds in this analysis are shared with MicroBooNE's search for NC Δ radiative decay [4]. The dominant contributions to the uncertainty on the background event rate for each analysis are from FSI related to inelastic nucleon scattering, pion and nucleon absorption, and pion charge exchange. The axial and vector mass parameters, m A and m V , respectively, in the charged current resonant form factors are also sources of significant uncertainties; this is consistent with expectation because of the large background due to charged-current interactions in which a π 0 is produced.
The efficiency of the selection is constructed using as the numerator the number of signal events passing all reconstruction cuts and analysis BDTs in simulation and as the denominator the total number of signal events preceding the application of any cuts or analysis BDTs. The difference in signal definition between the semiinclusive measurement and each of the two exclusive measurements is contained in the efficiency denominator. The exclusive measurements and the semi-inclusive measurement each use a distinct efficiency denominator, reflecting the total number of simulated events truly satisfying the corresponding signal definition. In each of the exclusive measurements, the signal definition is taken to be NC1π 0 with exactly zero or one final-state proton with a kinetic energy above 50 MeV. In the semi-inclusive measurement, the signal definition is taken to be NC1π 0 , notably allowing for any number of protons in the final state. The efficiency for each analysis is reported in Table II. The integrated flux is calculated separately for all four neutrino species (ν μ ,ν μ , ν e ,ν e ), and the sum of these integrated fluxes is used to normalize each cross-section measurement. This choice was made because of the inability to identify the species of the incident neutrino based on the neutral current final state. The integrated flux is varied within each flux systematic universe, and the correlations between each varied flux and the corresponding variations in the predicted background and efficiency are taken into account when extracting the cross sections.
The number of argon atoms used is calculated as N targets ¼ ρVN A =M Ar , where V ¼ 5.64 × 10 7 cm 3 is the fiducial volume of the analysis, ρ ¼ 1.3954 g=cm 3 is the density of argon at the temperature in the cryostat, and M Ar ¼ 39.948 g=mol is the molar mass of argon. A 1% uncertainty is assigned to the number of targets to reflect variation in the argon density through temperature and pressure fluctuations.

B. Results and interpretation
The calculation of each cross section from its components follows from Eq. (5) and is summarized in Table II. The resulting cross sections are shown in Fig. 12, compared to the simulated cross sections from several neutrino event generators including GENIE, NUWRO [53], and NEUT [54]. NEUT and GENIE both use the Berger-Sehgal model as their foundation for modeling pion production in the Δð1232Þ resonance region, while NUWRO implements a custom model optimized for the Δ resonance peak region. The complete details of the models used in each generator are discussed at length in Ref. [55]. The GENIE curve shown is generated using the MicroBooNE cross-section "tune" [29], which does not modify the GENIE v3.0.6 central value prediction (because the tune did not adjust the NC interaction model), but does define the uncertainty on the prediction. The error bars on the data points include systematic error associated with the modeling of background events that enters into the cross sections via background subtraction as well as error associated with the modeling of the signal events that enters into the cross sections via efficiency correction. The shaded band around the GENIE central value prediction shows the error associated with the prediction of the signal cross section.
We observe a consistent deficit in data compared to GENIE for the combined semi-inclusive measurement and for each of the individual NC π 0 þ 1p and NC π 0 þ 0p exclusive measurements. Overall, the NEUT predictions most closely match the reported measurements across semi-inclusive and exclusive final states. Additionally we note that while NUWRO is generally consistent with the other generators in its semi-inclusive and exclusive 1p predictions, its exclusive 0p prediction is higher compared to NEUT and GENIE predictions. The extracted semi-inclusive NC π 0 cross section is 1.24 AE 0.19ðsystÞ AE 0.08ðstatÞ [10 −38 cm 2 =Ar] which is 26% lower than the GENIE prediction of 1.68 [10 −38 cm 2 =Ar]. We calculate a χ 2 test statistic comparing our data measurement with full uncertainties to the CV of each model prediction for both exclusive NC π 0 þ 0p and NC π 0 þ 1p cross sections simultaneously (i.e. two degrees of freedom). The resulting values are 7.6, 7.7, 2.4, and 5.1 for comparisons against GENIE v3, GENIE v2, NEUT, and NUWRO, respectively.
The corresponding breakdown of uncertainty for each of the measurement channels is shown in Fig. 13. In all cases the flux, GENIE, and statistical uncertainties are dominant. The dominant contributions to the GENIE uncertainties enter into the cross section via the background subtraction and, as noted above, arise from the modeling of final-state interactions and the axial and vector mass parameters governing CC resonant pion production.
To further understand this measurement, it is instructive to compare it to previous experimental measurements of NC π 0 production. We compare our measurement to that performed by MiniBooNE which operated in the same beamline as MicroBooNE but which utilized a different detector material (mineral oil, CH 2 ) as the neutrino scattering target. In MiniBooNE's NC π 0 analysis, they measured NC interactions wherein only one π 0 and no additional mesons exited the target nucleus (no requirement on the number or identity of outgoing nucleons was made). A final flux-averaged cross section of 4.76 AE 0.76 AE 0.05 [10 −40 cm 2 =nucleon] was reported [25]. We can compare this result to our semi-inclusive result by comparing each to the same neutrino generator. This is shown in Fig. 14 where we compare both to the default GENIE v3.0.6 on argon and mineral oil respectively. We observe that while this result on argon lies slightly below the expected central value, both our result and MiniBooNE's agree with GENIE v3.0.6 within assigned uncertainties. FIG. 13. Error budget for the semi-inclusive NC π 0 , exclusive NC π 0 þ 1p, and exclusive NC π 0 þ 0p cross-section measurements.
FIG. 14. Comparison of this semi-inclusive result on argon (left) as well as that from MiniBooNE on mineral oil CH 2 (right), to the same GENIE v3.0.6. While the published MiniBooNE result was originally compared to a prediction made using the NUANCE v3 generator [56], we have instead generated a prediction using GENIE to aid in a comparison between the two experimental results. MiniBooNE's statistical uncertainty is small and only the systematic error bar is visible. Shaded error bands show GENIE uncertainty only.
FIG. 12. Measured semi-inclusive NC π 0 , exclusive NC π 0 þ 1p, and exclusive NC π 0 þ 0p cross sections, each compared to the corresponding GENIE v3 (G18_10a_02_11) cross section and its uncertainty (shaded red bands) as well as other contemporary neutrino generators. Inner error bars on data points are statistical only; outer are statistical and systematic, summed in quadrature.

V. SUMMARY
In summary, we report the highest-statistics measurement to date of neutrino neutral current single pion production on argon, including the first exclusive measurements of this process ever made in argon. These cross sections are measured using the MicroBooNE detector exposed to the Fermilab Booster Neutrino Beamline, which has hE ν i < 1 GeV. As presented within this paper, kinematic distributions of the π 0 momentum and angle relative to the beam direction provide some sensitivity to contributions to this process from coherent and noncoherent pion production and suggest that, given the currently analyzed MicroBooNE data statistics, the nominal GENIE neutrino event generator used for MicroBooNE Monte Carlo modeling describes the observed distributions within uncertainties. This has provided an important validation check justifying the use of this sample as a powerful constraint for backgrounds to single-photon searches in MicroBooNE, e.g. in [4].
Using a total of 1,130 observed NC π 0 events, a fluxaveraged cross section has been extracted for neutrinos with a mean energy of 804 MeV and has been found to correspond to 1.243 AE 0.185ðsystÞ AE 0.076ðstatÞ, 0.444 AE 0.098 AE 0.047, and 0.624 AE 0.131 AE 0.075 ½10 −38 cm 2 =Ar for the semi-inclusive NC π 0 , exclusive NC π 0 þ 1p, and exclusive NC π 0 þ 0p processes compared to 1.678, 0.722, and 0.774½10 −38 cm 2 =Ar in the default GENIE prediction used by MicroBooNE. Comparison to other generators including NEUT and NUWRO show reasonable agreement with the NEUT predictions found to be slightly more consistent with the MicroBooNE data-extracted cross section for all three exclusive and semi-inclusive processes.