Neutron measurements from anti-neutrino hydrocarbon reactions

Charged-current anti-neutrino interactions on hydrocarbon scintillator in the MINERvA detector are used to study activity from their final-state neutrons. To ensure that most of the neutrons are from the primary interaction, rather than hadronic reinteractions in the detector, the sample is limited to momentum transfers below 0.8 GeV/c. From 16,129 interactions, 15,246 neutral particle candidates are observed. The reference simulation predicts 64\% of these candidates are due to neutrons from the anti-neutrino interaction directly, but also overpredicts the number of candidates by 15\% overall, which is beyond the standard uncertainty estimates for models of neutrino interactions and neutron propagation in the detector. Using the measured distributions for energy deposition, time of flight, position, and speed, we explore the sensitivity to the details those two aspects of the models. We also use multiplicity distributions to evaluate the presence of a two-nucleon knockout process. These results provide critical new information toward a complete description of the hadronic final state of neutrino interactions, which is vital to neutrino oscillation experiments.


I. INTRODUCTION
Neutrons are the last essential (anti-)neutrino interaction final-state particle to have their number and energy distribution studied in detail. Neutrons play a special role in those oscillation measurements which depend on comparing distortions of the antineutrino energy spectrum to the neutrino energy spectrum [1][2][3][4][5][6]. The antineutrino reactions' calorimetric response is heavily suppressed relative to the neutrino case because of the prevalence of neutrons, consistent with their generic role in hadron calorimetry [7]. The neutron energy deposits are not proportional to their kinetic energy and are not always observed at a location consistent with their initial trajectory from the interaction point. Sometimes they are not observed at all because they escape the detector volume unseen, their interaction products are below Cerenkov detection threshold, or they thermalize without producing enough ionization to be reconstructed. In addition to energy determination, neutrons' presence in the final-state also impacts event selection and background rejection for oscillation, interaction, and rare process analyses.
In this paper, we present the first direct measurements of the neutron content from neutron-rich antineutrino charged-current reactions. Neutrons with ten to hundreds of MeV kinetic energy are observed as they rescatter off hydrogen and carbon in the MINERvA detector. This study is done with a low three-momentum transfer sample [8,9] but is otherwise inclusive with no selection on the number or type of hadrons in the finalstate. Low overall hadronic activity means that neutron activity is easily separated and likely to be due to neutrons from the original neutrino interaction rather than secondary hadronic reactions. The results are compared to the predictions of a full simulation that consists of a modified version of the GENIE neutrino event generator [10], a GEANT4 simulation of particle propagation in the detector material [11,12], and a simulation of the calibrated response of our scintillator and electronics [13,14].
The energy deposit, time, position, speed, and multiplicity distributions are sensitive to the details of neutron production in the initial reaction and to models for neutron scattering in the detector. When the neutrino reaction occurs in a nucleus, modeling the hadronic finalstate is complicated. Different charged-current weakinteraction processes produce different numbers of neutrons after the final-state lepton gains its charge. The charge-changing antineutrino quasielastic (QE) and twoparticle two-hole (2p2h) process must turn at least one proton into a neutron. Often the neutron has all of the hadronic final-state energy. In resonance production and deeply inelastic scattering (DIS) the charge can be exchanged with the resulting meson or the struck nucleon or quark, producing a higher number of neutrons per event than the equivalent neutrino case. These outcomes for carbon are summarized in Table I. For antineutrino reactions on hydrogen, the 2p2h and FSI processes do not occur; also charged-current neutrino-hydrogen QE reactions are not possible. neutron content process antineutrino neutrino QE 1 0 resonance 1 or 0 0 or 1 2p2h 2 or 1 1 or 0 after FSI 0 to 7 0 to 5 Sensitivity to the mix of reaction types is reduced and indirect because the resulting hadrons will frequently reinteract on their way out of the nucleus. In event generators, these rescattering processes are referred to as final-state interactions (FSI). Such reinteractions can be soft scatters that do not change the outgoing charge state, full charge exchange reactions where the energetic neutron becomes a proton or vice versa, and knockout reactions where multiple nucleons and mesons exit the nucleus. Calorimetric measurement of the interaction is affected and the hadron topology also changes.
Neutrons from neutrino and antineutrino reactions have been measured before. The earliest technique was tagging the capture of thermal and "fast" (up to 10 MeV) neutrons, used from the very first neutrino observations [15] to the present and near future [16][17][18][19]. Higher energy neutrons from neutrino reactions caused the most important np → np background to the discovery [20,21] of weak neutral current reactions in the freonfilled Gargamelle detector at CERN. The collaboration made measurements of "associated neutrons" and used a cascade simulation [22] to translate that measurement into an estimate for neutron production from neutrinos interacting in the material upstream, and so obtained the crucial constraint on the background. Similar studies were needed for followup neutral current measurements including those with two liquid scintillator detectors [23][24][25] in the Brookhaven National Lab neutrino beam, and the "dirt backgrounds" from MiniBooNE's measurement [26,27]. These neutrons' time and spatial distributions in the detector were measured and simulated in order to constrain and subtract backgrounds, but were not correlated with simultaneous measurement of their original interaction in the material upstream of the detector. MiniBooNE's paper presents a comparison to a modern Monte Carlo simulation using the Nuance neutrino event generator [28], and GEANT3 [29] using the GCalor [30] option, and found that a 30% reduction of the neutron component was needed to describe the data.
Recently, a measurement [31] was presented by Ar-goNeuT for neutrino-argon reactions which is similar to the one in this paper. The ArgoNeuT analysis studied low-energy photons produced by the deexcitation of the argon nucleus struck by the neutrino and by the deexcitation of argon nuclei struck by neutrons generated in the neutrino reaction. Unlike the MINERvA antineutrino data described in this paper, the neutron component is only half the sample, and the photons from deexcitation of the argon nucleus struck by the neutrino are evident in the spatial and multiplicity distributions. In the ArgoNeuT study, both components were simulated by FLUKA [32,33]. Another recent paper [34] breaks down the simulation of neutrino and antineutrino reactions with argon to provide details of the pathways to missing energy, especially via neutrons.
MINERvA has several advantages relevant to detecting neutrons from (anti-)neutrino interactions. The 5.3-ton, fully-active tracking volume is much larger than the neutron interaction length of approximately 10 cm for neutron kinetic energies near 20 MeV. Larger volumes have been used (Super-K, MiniBooNE) but their Cerenkov technique has too high a threshold for detecting the protons scattered by neutrons. MINERvA's active volume is polystyrene (C 8 H 8 ) n , shortened to CH when used in nuclear and particle physics. For low-energy neutron detection, the hydrogen presents as significant a target as the carbon, despite a carbon nucleus having twelve nucleons. And MINERvA is sensitive down to a low threshold of 1 MeV because it is an underground detector with low noise overall and well-constrained contributions from other sources of neutron production in the beam.

II. THE MINERVA EXPERIMENT
MINERvA is a dedicated neutrino-nucleus cross section experiment. Its goals are to make cross section measurements needed for neutrino oscillation experiments and to probe the environment of the nucleus, complementary to what the electron scattering community has accomplished. The experiment is located in the highintensity NuMI neutrino beam at Fermilab.
The centerpiece of the detector [13] is a 5.3-ton fully-active scintillator tracker with excellent calorimetric properties of its own and surrounded by additional electromagnetic and hadronic calorimeters (ECAL and HCAL respectively). The experiment also has passive targets made of iron, lead, water, graphite, and helium, which enable the study of the A-dependence of neutrino and antineutrino reactions.
The tracker is built from planes of polystyrene scintillator strips. With Lexan sheets, epoxy, tape, and reflective titanium dioxide, the target consists of 8.2%, 88.5%, and 2.5% hydrogen, carbon, and oxygen respectively (by mass), plus small amounts of heavier nuclei. The strips are triangular in shape with 3.3-cm base and 1.7-cm height and up to 245-cm length. The strips are nested with alternating orientation to make 1.7-cm thick planes. With this arrangement, ionization activity in the plane is split between strips and tracking resolution is improved. An entire plane is a hexagon containing 127 strips and one module consists of two planes. The second plane in one module is oriented with the strips vertical, producing an X-coordinate of the observed energy deposits. The plane in front of it is rotated 60 degrees one way to form a U coordinate, or the other way to form a V coordinate. The resulting modules are themselves alternated to produce repeating UX,VX sets of planes, with the detector Z axis running normal to these planes. The MINOS Near Detector [35] is located 2 m downstream of MINERvA and measures the charge-sign and momentum of muons selected in this analysis.
These data were obtained from the NuMI beam [36] operating in antineutrino mode. The primary 120-GeV proton beam interacts in a graphite target producing pions and kaons. Two magnetic horns focus the the negativelycharged mesons toward a decay pipe, leading to an antineutrino spectrum that peaks near 3.0 GeV. In total, these data are from an exposure of 1.04 × 10 20 protons on target between November 2010 and February 2011. This beam configuration also produces neutrino chargedcurrent interactions in the detector which are over 10% of the events in the antineutrino peak energy range used in this analysis. These are removed with high efficiency because their measured curvature in the MINOS Near Detector is the wrong direction.
The flux prediction is GEANT4-based [11,12] with central values and uncertainties adjusted [37] using thintarget hadron production data [38][39][40][41] and an in situ neutrino-electron scattering constraint [42]. The design of this analysis is relatively insensitive to the resulting 8 to 10% uncertainties in the energy spectrum or the absolute flux.

III. SIMULATION
The reference simulation combines the GENIE 2.8.4 neutrino interaction model [10] with modifications, the GEANT4 9.4.p2 particle transport model [11,12] with modifications, a calibrated detector time and energy response model [13,14], and an "overlay" of data events to reproduce the unrelated activity that might overlap in time with the simulated event.
A. GENIE event generator GENIE's simulation of the charged-current quasielastic process is from Llewellyn Smith [43] with vector form factors parameterized as in [44], and a dipole form factor with axial mass of 0.99 GeV. A relativistic Fermi gas model [45] is implemented for interactions on carbon and other nuclei. The ∆ and higher resonances are from Rein and Sehgal [46], with a nonresonant component added from the DIS model as the resonances are phased out from invariant mass 1.4 < W < 2.0 GeV. The DIS cross section comes from the Bodek-Yang model [47], where the hadronic system [48] is produced using KNO scaling [49] transitioning to PYTHIA [50] between 2.4 and 3.0 GeV.
Modifications are made to the above default GENIE 2.8.4 model. We refer to this set as MINERvA tune version 1.1 (MnvGENIE-v1.1). These modifications can also be applied to the default 2.12.6 version of GENIE. The QE process is modified to include a screening effect based on the random phase approximation (RPA) technique. The suppression is based on the calculations of Nieves and collaborators [51,52] for a Fermi-gas nucleus and implemented by reweighting GENIE QE events [53], including an uncertainty on the RPA screening derived from comparison to neutron capture data [54,55]. The uncertainty on the QE axial form factor is set to 9% following the analysis of [56]; additional uncertainty on the highest Q 2 component is not needed for this sample. Nonresonance pion production is reduced based on the reanalysis of bubble chamber neutrino data [57,58]. Coherent pion production with pion kinetic energy below 450 MeV is also reduced based on analysis of MINERvA data [59,60] and consistent with the Berger-Sehgal [61] modifications of the original Rein-Sehgal model [62].
Of special interest is the interaction of the antineutrino with two nucleons, knocking them both out and leaving two holes (2p2h) in the nucleus. In this analysis, the base model for the 2p2h component with no pion is from the IFIC Valencia group [52,63] implemented in GENIE [64]. This process is further enhanced in the region between the QE and ∆ resonance components based on a fit [65] to the reconstructed neutrino data presented in [9]. In all but the last figure in this paper, the error band includes an uncertainty on this fit that varies the fraction of reactions on pn and pp initial states. This enhancement has been applied to this analysis and successfully describes a wide variety of data for other MINERvA observables [8,[66][67][68][69][70].
Reinteractions of nucleons and mesons (final-state interactions or FSI) as they exit the target nucleus are the other important mechanism for neutron production. The GENIE simulation of final-state reinteractions of hadrons leaving the nucleus is a parameterized, effective cascade model which is called "hA", short for hadron-nucleus interaction. The model steps hadrons through a nucleus with its associated radius and nuclear density function. The hadron's mean free path is then determined from tabulated hadron-proton and hadron-neutron cross sections from SAID [71]. The resulting probability to interact within the nucleus is high, 73% for a neutron from a 3 GeV quasielastic event in carbon, and 88% in iron. Because the code originated with the MINOS experiment [35], when an interaction is specified, the fates (absorption, pion production, knockout, charge exchange, elastic scatter) are chosen according to their proportions for iron. In-medium effects that might relatively favor particular fates in carbon or lead are not included in this original version, nor is multiple scattering. When one interaction occurs, the number of outgoing nucleons is drawn from a distribution that favors single nucleon, deuteron, and alpha states, but allows a chance for complete breakup. Compared to the cross sections for the initial reactions, the simulation of FSI differs the most among different event generators. The GENIE FSI model produces more low-energy nucleons than other commonly used neutrino event generators. Figure 1 shows the energy and multiplicity spectra for a newer (2.12.10) but equivalent version of GENIE with its models and modifications configured for the MnvGENIE-v1.1 tune. It is compared to two other common generators NEUT [72,73] and NuWRO [74], produced using the NUISANCE framework [75]. The inputs are monoenergetic 3-GeV antineutrinos on a CH target, and the q 3 < 0.8 GeV selection is made. Neutrons below 2 MeV are not included in the kinetic energy distribution. Appropriate to this analysis, neutrons below 10 MeV are not included in the multiplicity distribution. Below 25 MeV, the NEUT generator stops any reinteraction process (cross section is set to zero), so those neutrons are placed outside the nucleus rather than producing yet more and lower energy neutrons. In contrast, above 50 MeV, each generator is using similar models for the primary production of neutrons before FSI, and their predictions for the neutrons' kinetic energy spectra converge.

B. GEANT4 neutron propagation model
Particles in the final state are passed to GEANT4 9.4.p2 with the Bertini cascade option (QGSP BERT) which propagates them through a model of the MIN-ERvA detector geometry and materials. Uncertainties in the interaction cross section in scintillator are implemented using a reweighting scheme that considers proton and charged pions (±10%), and neutron (±15-25%) scattering separately. Confidence in the proton and pion response at these energies also comes from analysis of calorimetric data taken with MINERvA detector elements in a test beam at Fermilab [14].
The energy dependence of the total cross section for neutrons in GEANT4 has been improved since this older version. We have compared this version with a late 2016 release (v10.2.p2) to understand changes in the neutron interaction code. The new versions now match the Abfalterer et al. high-precision neutron scattering total cross sections [76] on carbon and hydrogen from 5.2 to 560 MeV. Using the same reweighting tools as to evaluate hadron systematic uncertainties, the older GEANT4 cross section is changed to represent these same data. After making this correction we assign uncertainties of 25% below 10 MeV, 20% from 10 MeV to 25 MeV, and 15% above 25 MeV. These are applied to the total cross section, but are larger than the remaining discrepancy with the data. In this analysis they play the effective role of uncertainties on the elastic vs. inelastic components.
In principle, another uncertainty arises from the outcomes of the neutron-hydrogen and neutron-carbon interactions. The Abfalterer et al. total cross section includes elastic scatters that deflect the neutron as little as 0.12 degrees, which involve sub-keV scale energy transfers. Thus, even with an accurate total cross section, the GEANT4 model also makes a prediction for the fraction of scatters that are above and below our experimental threshold.
A second, prominent feature of the GEANT4 neutronscattering model is the production of deexcitation photons following neutron-carbon reactions. Such photons Compton scatter and account for 50% of the neutron candidates that originate from a GENIE neutron. The simulation predicts the rest are from protons and up to 5% from nuclear fragments. For comparison, two recent accountings of neutron induced activity are presented in [31,34] for argon; the latter has especially detailed discussion of simulated activity.
GEANT4 provides an alternate, high precision neutron simulation (HP in the GEANT4 option names), originally designed for studies of fission reactions up to about 20 MeV. After accounting for the Abfalterer cross sections, the differences between fully simulated HP configuration and the default configuration appear where the lowest-energy neutrons are substantial parts of the sample, especially near the interaction point. But these differences are modest, similar to the other uncertainties, and difficult to disentangle from the rest of the GEANT4 predicted response.

C. Detector response
The simulation also must produce energy deposits consistent with our calibrated scintillator and electronics response. The photoelectron yield and absolute energy scale is tuned using a comparison of data and simulated muons at near-normal incidence to the planes. Nonlinearities for energy deposits above minimum ionizing muons are accounted for via individual calibration of the digitizing electronics [13] and by the test beam calibration [14] of Birks' quenching in scintillator.
The simulation reproduces the detailed response of light propagation in the scintillator bars and the photomultiplier tube, so that the simulated activity is reconstructed using identical steps as with the data. When a particle is fully tracked, every energy deposit's location is known in all three dimensions. In this case, the reconstruction estimates the effect of light reflection and attenuation in the scintillator strip and optical fiber and produces a more accurate estimate of the actual energy deposit. This is rarely possible for the neutron energy deposits in this analysis. When not a part of the track, the reconstruction approximates each energy deposit to happen at a position halfway along the scintillator bar. Geometrical fluctuations will therefore be present in the energy and timing distributions in both data and simulation.
The width of the time distribution of digitized activity after light propagation depends on the photo-electron yield in the PMT. The distribution used in the simulation is based on the observed distribution in data from fully-tracked muons. Because averaging over the entire track yields a precise time and location for the passage of the muon through any one plane, the correlation of the fluctuations in time and the light yield in a single scintillator strip can be obtained directly, without resorting to individual photon modeling of the optical elements and photomultiplier tube response. The width of the time distribution is 10 ns for 1 to 2 photoelectrons, 3 ns for 6 to 12 (about 1 to 1.5 MeV in a single strip), then approaches the electronics limit of 2.2 ns. More detail on the nanosecond timing response and its use can be found in [77]. With geometry-based fluctuations, the final time resolution for neutron candidates is 4.5 ns.
A subdominant contribution to neutron candidates in this analysis comes from other neutrino interactions inside and outside the detector which produce their own neutrons and photons. They randomly happen at the same time as the charged-current reactions selected for this analysis. These "accidental" backgrounds are added directly to the simulation; 16 µs of activity from one pulse of the beam from data are selected randomly from the same months. This data activity is added on top of the activity from the simulated event. The resulting set of reconstructed times, locations, and energies are given to the same reconstruction algorithm. Thus this background and its dependence on the intensity of the beam are reproduced. Using alternate selections to isolate four regions that are predicted to be high (> 70%) in this particular background leads to the assignment of a ±10% uncertainty applied to the simulation.

D. Further modifications of the simulation
The configurations within the simulation packages do not contain enough explicit uncertainties or tunable parameters to describe the data presented in this paper. Described in detail with the data in later sections, we will make two heuristic reductions in the number of neutron candidates in the simulation. One mimics a mismodeling of either or both the energy deposit spectrum and the number of deexcitation photons as simulated by GEANT4. The other reflects the wide range of predictions for low-energy neutrons from the neutrino interaction models suggested by Fig. 1. For this paper, the reductions allow us to quantitatively explore possible unknown effects, even though they do not represent a ±1σ uncertainty.

IV. EVENT SAMPLE AND NEUTRON SELECTION
The sample of charged-current antineutrino events analyzed in this paper is the same inclusive sample from [8] with reconstructed three-momentum transfer q 3 less than 0.8 GeV/c. This sample has high neutron content but little other charged hadronic activity to complicate the analysis of neutron activity. This allows us to use a low threshold of 1.5 MeV for neutron candidates.

A. Selection of antineutrino events
Charged-current antineutrino interactions originate in the scintillator in the active tracker fiducial region. The resulting µ + must be fully tracked to the end of the MIN-ERvA detector and also reconstructed in the MINOS Near Detector [35] where its positive charge and momentum are analyzed. To ensure a region of good geometrical acceptance we require p µ > 1.5 GeV/c, and θ µ < 25 degrees. The selected reconstructed neutrino energy range is 2.0 < E ν < 6.0 GeV, so resulting data and simulated samples are the lower panels of Fig. 3 of [8] and also match the selection used for the related neutrino-mode analysis [9].
The muon energy and angle are combined with the observed hadronic energy deposits to form calorimetric estimates for the energy transfer q 0 (often called ω or ν by different groups in the literature), neutrino energy µ , and the magnitude of the three-momentum transfer q 3 = Q 2 + q 2 0 (often simply called q or |q|). The selection is inclusive because only the magnitude of reconstructed q 3 < 0.8 GeV/c enters the selection, not details of the number or type of hadrons observed.
We exploit the feature that this subsample can be divided into regions hereafter called QE-rich, dip, and ∆-rich. The QE-rich region has little or no observed hadronic energy and is predicted to be mostly 2p2h and QE. The ∆-rich region has the most energy transfer and is predicted to be mostly resonance production and some 2p2h. The so-called "dip" region in between these two is a mix of all three processes, but also has the highest predicted concentration of 2p2h events. Again referring to the lower panels of Fig. 3 in [8], boundaries are formed between the QE-rich, dip-region, and ∆-rich subsamples. The reconstructed "available energy" is an estimator for a quantity that includes proton and charged pion kinetic energy and the total energy of neutral pions, photons, and electrons produced by a neutrino interaction model. It explicitly does not include kinetic energy of neutrons nor the energy cost to remove nucleons from the nucleus, so is always lower than the true energy transfer, and for some antineutrino reactions will be zero. The boundaries are at reconstructed available energy of 0.06 and 0.12 GeV for 0.0 < q 3 < 0.4 GeV/c and 0.08 and 0.18 GeV for 0.4 < q 3 < 0.8 GeV/c. For brevity many distributions in this paper are not divided this way and are presented as two q 3 regions.

B. Selection of neutral particle induced candidates
Because the sample is limited to q 3 < 0.8 GeV/c, charged hadron activity is small and remains localized to the interaction point. The rest of the detector should have no hadronic activity except neutrons and photons, which this analysis considers signal. There are two backgrounds to neutron candidates to minimize: the "muon background" from electrons and bremsstrahlung photons, and the "accidental background" from activity induced by the neutrino beam or cosmic rays and unrelated to the antineutrino reaction being analyzed.
Analysis of hadronic activity is done with reconstructed objects formed of "clusters" of energy deposits in one or multiple adjacent strips within the same plane. Each cluster is assigned a calibrated energy based on the sum of the individual energy deposits. A cluster also has a two-dimensional position, one from the location of the plane in the detector and one from estimating the transverse position in its plane based on energy-weighted average positions of activity in multiple strips. The calibrated time is the weighted average of hit times, taking into account the measured correlation between the number of photoelectrons and the width of the timing distribution of single strip activity from muons in data.
The criteria for spatially-isolated clusters are designed to exclude four overlapping volumes in spatial proximity to the muon activity, the interaction point, other charged hadron activity, or to the outer boundary of the detector. An event display of the X-coordinate activity simulated event, Fig. 2, summarizes how activity in the first three categories is rejected and a single cluster neutron candidate is observed. Additional energy threshold and timing cuts complete the selection. Remaining clusters that are near each other, indicating they may be caused by the same neutral particle, are combined into multicluster candidates. The simulation is used to evaluate the effectiveness of these selections at reducing the muon and accidental backgrounds. interaction point, and other charged hadron activity (a π − in this simulated event) with the remaining activity promoted to a neutron candidate. The aspect ratio for this figure exaggerates the transverse dimension by almost a factor of two in order to emphasize the detail. Activity from the π − in adjacent U and V planes near the interaction point is not shown.

Muon exclusion zone
The muon itself is fully tracked for all events in this sample, and most clusters are already assigned to the reconstructed track. Additional activity near the muon track is very likely caused by photons from bremsstrahlung and knock-on electrons (delta-rays). From one module upstream of the interaction point to the back of the detector we exclude clusters within 17 cm (about 10 strips) of the muon from consideration, about one mean free path for photons with energy of a few MeV. The simulation predicts that anti-muon induced activity accounts for 90% of what would otherwise be candidates in this zone, but the simulation also underpredicts the data by 12%. This exclusion reduces the muon-induced background by a factor of 10.
This zone is increased to 24 cm starting 20 modules downstream from the interaction point. The additional volume is predicted to be 65% muon activity and 3% accidental backgrounds and excluding it further reduces the muon background by another factor of two. The simulation describes the activity in this outer zone well, contributing 3.5% fewer clusters than data. We assign twice this difference as the uncertainty on the muon contribution to the remaining selected clusters, which is negligible for this analysis.

Interaction point exclusion zone box
Charged hadrons will produce activity near the neutrino interaction point and the start of the muon track. The granularity of the detector imposes limits on separating muon, hadron, photon, and cross talk activity when some or all of it fails to meet tracking criteria. Clusters of activity within a transverse "box" around this point are excluded from further consideration. Because of the three X,U,V orientations, the box has a pinwheel shape, but can be defined and coded simply when three-dimensional reconstruction is not available. This exclusion is for clusters within 17 cm transverse from a horizontal line parallel to the detector Z axis through the start of the muon track from ten planes upstream to fifteen planes downstream. This region remains rich with neutron activity and will be discussed later with the results.

Charged particle exclusion zone
Some charged hadrons travel outside this box. A spatial algorithm considers seed clusters near the interaction point and adds additional clusters to the total charged particle activity if they are nearby. This algorithm will follow such spatially-connected activity arbitrarily far from the interaction point. This includes many π − and a few protons that are fully tracked, but also activity from reinteractions in the detector that do not satisfy stricter requirements to form or extend a track object. This is the first full deployment of such an algorithm in a MINERvA analysis.

Edge of detector, timing, and energy
Activity at the edges of the detector is excluded as follows: the first 20 planes and veto wall upstream, the last 10 planes in the downstream hadronic calorimeter, and the entire outer detector hadronic calorimeter. Most clusters are naturally found in the inner tracker region because it is near to the interaction point, has the largest fraction of the active scintillator in the detector, and the scintillator has the highest hydrogen content of all the detector materials.
A few additional selections reduce the already small accidental background.
Clusters not yet excluded must be within a time window from 20 ns before to 35 ns after the interaction time t 0 determined largely from the muon track timing information. The clusters must also have at least 1.5 MeV of energy. The accidental background overtakes the predicted signal processes at cluster energies below 1.2 MeV. This energy cut also eliminates photomultiplier tube cross talk effects in both selected data and simulated events. This version of the accidental background overlay technique is not perfect. For this early version it is checked against the data for accidentalrich subsamples leading to an uncertainty of 10%. This conservative uncertainty has negligible impact on the results.

Aggregating spatially nearby clusters into candidates
Some isolated clusters can be spatially connected. Two clusters within three modules of each other are merged into one neutron candidate. This merging continues if additional clusters satisfy this requirement. Because there is so little hadronic activity in these low q 3 subsamples this simple requirement is effective. In this sample, 62% are candidates made of a single cluster, 18% are made of exactly two clusters; the simulation predicts 61% and 19% respectively.

C. Signal and background
Neutrons that exit the nucleus where the neutrino interaction occurred are of the most interest. They are the aspect of the interaction model that has never before been directly tested. When referring to the GENIE simulation of this component, including direct production and via FSI, we will call them GENIE neutrons. Protons and charged pions produce real secondary neutrons as they travel through the detector. These are an irreducible background and no attempt is made to reject them. Table II quantifies the sources that cause the majority of neutron candidates.
Photons also produce energy deposits distant from the main parts of the event. Thus the decay products of π 0 are an electromagnetic component of candidates. The π − entries in Table II includes 14% (so 1.8 and 3.2% of the respective column totals) which charge exchanged to π 0 in the detector. With such little hadronic energy in the final state, it is convenient to treat these also as irreducible backgrounds, and make no selection to reduce them. Because of the q 3 selection, these photons are relatively low energy. MINERvA's γ and π 0 selection algorithms [67,78] have an efficiency of 15%, because π 0 identification requires two reconstructed photon candidates and an effective threshold of 25 MeV on the lower  energy photon. As will be seen, the electromagnetic component traveling at the speed of light will separate from slower neutrons in distributions that use time-of-flight information. Because no attempt is made to reject the π 0 backgrounds, the analysis of neutral particle activity remains inclusive of all hadronic final-states. The small "other" category originates with π + and kaons, mostly from DIS processes which feed down into the sample from higher q 3 reactions. This category also includes photons and electrons from η 0 production and ∆ → N γ decay. Photons from GENIE's deexcitation of residual nuclei would be present and isotropic in the data sample, but are only simulated for oxygen, not carbon or other nuclei. Radiative processes from the lepton exiting the nucleus might be co-linear with the muon and are also not simulated in GENIE.
The backgrounds from the muon and accidentals are small. The effectiveness of the cuts can be illustrated quantitatively, relative to the base selection. Including the very ends of the detector doubles the accidental background. Not excluding the second, outer volume around the muon doubles the muon background. Using an incident neutrino energy range up to 20 GeV increases the statistics of the total data and simulated sample by 17% and also all the simulated signal and irreducible background subcomponents. Compared to the base selection, these higher energy reactions have 30% more in the other category from feed-down of higher q 3 DIS events. The muon component (because Bremsstrahlung is more likely) is also higher by 22%.
Three other cuts could be relaxed to extend the physics reach, but doing so would increase the backgrounds. These samples are consistent with the main sample but do not further enhance the conclusions. Reducing the energy threshold for candidates from 1.5 MeV down to 1.2 MeV increases the predicted rate by 10% overall, but the data rate increases only by 8%. The predicted muon and accidental background rate increase by 30% and 50% respectively, collectively accounting for 20% of the additional candidates. Allowing the timing to go out to 100 ns adds events that are predicted to be half from the accidental background.

D. Efficiency
Focusing again on the most interesting signal process, the probability that a neutron from the GENIE model produces a cluster of activity is high because of the large volume of fully-active scintillator. The probability to survive the selection rises with kinetic energy from few to 60% . This is shown in Fig. 3, which also illustrates the predicted neutron kinetic energy spectrum for this antineutrino q 3 < 0.8 GeV/c sample. Details of the se- lection process sculpt the distribution. The region near the interaction point especially is a place where 45% of neutrons leave activity, as one would expect from a meanfree-path process. Some lower energy neutrons are effectively below threshold after just one interaction, and the lowest energy neutrons do not travel very far and thermalize locally. Higher energy neutrons may still interact again and be visible. The efficiency for detecting neutrons using plastic scintillator has been studied since the late 1950s. The 1962 measurement [79] in scintillator 15 cm thick with a simple discriminator-based 1.5 MeV electron-equivalent threshold is very similar to the conditions of this MINERvA measurement. They observed an efficiency of 20% at 76 MeV and 30% at 10 MeV. Below this the threshold reduces the efficiency significantly. These data have been compared to increasingly sophisticated simulations over the years, such as [80]. GEANT4 produces a similar efficiency for neutrons above 20 MeV, but predicts closer to 40% efficiency at 10 MeV.
The mix of hydrogen and carbon in the detector affects the efficiency for different ranges of neutron kinetic energy. A special purpose GEANT4 neutron simulation in the MINERvA detector with the same reconstruction as this analysis predicts that hydrogen produces more candidates up to neutron kinetic energy of 20 MeV where it is equally probable as detectable scatters from carbon. The hydrogen:carbon ratio remains around 1:6 up to 100 MeV, then falls further past 1:12 at higher neutron kinetic energies. Neutron inelastic scattering in carbon (including single nucleon knockout like the process called quasielastic electron scattering) begins around 10 MeV of neutron kinetic energy. At lower kinetic energy, neutrons are not able to transfer enough energy to a proton to remove it from the nucleus, yet elastic scattering off hydrogen can produce a proton above detection threshold. A liquid argon detector would have significantly less acceptance below 20 MeV kinetic energy for the same threshold; a liquid scintillator (CH 2 ) detector could have more.
Identifying the presence of a neutron is easier than guessing what the energy of the neutron was. The spectra of energy deposits for neutrons in three different kinetic energy ranges are shown in Fig. 4 (the same figure appears in [8]). The reconstructed energy of a single cluster, or the sum of two or more clusters aggregated into a single candidate is mostly uncorrelated with the kinetic energy of the neutron. The most likely neutron candidate for all kinetic energies considered in this sample is just at the 1.5 MeV threshold. As the neutron energy rises, the probability it will produce an energetic proton that travels several planes and leaves tens to hundreds of MeV also rises, and populates the distribution at and beyond the right edge of the plot up to hundreds of MeV.
The converse is also relevant. Neutrons in the range 10 to 20 MeV can only make the smallest energy deposits. The presence of a small energy deposit does not automatically indicate a low-energy neutron. But the presence of many low-energy neutrons can only enhance the rate of the smallest energy deposits.
Overall, the presence of neutrons down to 50 MeV kinetic energy is determined with good efficiency and low backgrounds. Though the efficiency continues to fall, neutrons down to 10 MeV are predicted to cause a substantial part of the sample.

V. MEASUREMENTS OF NEUTRON ACTIVITY
One observation already stands out from Table II. There is an overall overprediction of neutron candidates in the simulation, further emphasized in Table III. It appears in all subdivisions of the sample, despite different amounts of QE, 2p2h, and ∆ resonance, and their varying neutron final-state content. Because the efficiency is so high and the backgrounds are so low, either GENIE is producing too many neutrons per event, or the GEANT4 neutron propagation plus detector response simulation is making them more visible than in the data. These motivate the analysis strategy reported here. Compared to Table III  In this section, distributions of deposited energy, time, and position upstream or downstream are shown. Then neutron speed (actually 1/β) and multiplicity per event are used to draw final conclusions. Starting with the time distributions, the low-energy candidates and high-energy candidates are separated, and the oversimulation persists in the former.
In all cases, the reconstructed data distribution is shown with statistical uncertainties only (often too small to see) and the simulation is shown with model and detector systematic uncertainties. Some physics effects being tested with these data are not included in the systematic uncertainties and are described explicitly. Each figure is neutron candidates per event like Table III, reducing systematics that affect the numerator and denominator equally, such as the flux and some cross section uncertainties. Discussion of specific uncertainties and resolutions are provided with each new distribution, if not previously described.
In each plot, the simulation is broken down into GE-NIE neutrons, other sources of neutrons, and the electromagnetic component. The muon and accidental backgrounds are often too small to see and are described in the text instead of being plotted. Each figure has a lower, ratio subpanel where the data to simulation ratio is compared to the reference model and its uncertainties, to emphasize the magnitude and location of discrepancies.

A. Spectrum of candidate deposited energy
The energy deposit spectrum E dep shown in Fig. 5 highlights that the extra neutron candidates in the simulation are limited to those with less than 10 MeV. They are unambiguously from neutron production. The neutron components of the spectrum peak much more strongly at threshold than the electromagnetic component.
We use these data to probe for model features beyond the standard systematic uncertainties included in the error band. The ratio subpanel includes two modified models as simple benchmarks. For the solid "modified GEANT4" line, we eliminate a random 50% of neutron candidates with less than 10 MeV energy deposit and originating from neutron-producing GENIE particles (neutrons, protons, π ± ). Without regard to the energy of the particle, this mimics moving part of the GEANT4 or detector response below our detection threshold, making these elastic scatters invisible, or reducing the number of photons produced by nuclear deexcitations in carbon following the nucleon knockout process. The dotted line implements a neutrino interaction model change; 50% of candidates caused by GENIE neutrons below kinetic energy of 50 MeV are removed. This makes the prediction more like NuWRO in Fig. 1. These benchmark modifications are chosen empirically to have about the right size and allow the analysis to track their effects across the rest of the distributions.
Because the simulated overprediction is so prominent in the first two bins, the following distributions will be separated into subsamples below and above 10 MeV of candidate energy deposit E dep . Considering Fig. 4, the latter sample is necessarily higher kinetic energy neutrons, while the lowest energy candidates are a mix of everything. The backgrounds make a 12±1% contribu- tion (not shown) in the first bin. There will be essentially no backgrounds in the E dep > 10 MeV subsamples.

Uncertainties
Signal response uncertainties that affect the probability a neutron will produce a cluster near the 1.5 MeV threshold are important. Neutron elastic scatters produce a low-energy proton whose detector response is affected by a scintillator quenching effect parameterized using Birks' Law [81,82]. Our parameter and its uncertainties are calibrated using test beam data from stopping protons [14] especially the last 40 MeV of their energy deposits. The uncertainty is doubled for neutroninduced candidates with less than 5 MeV, which are not well constrained by the test beam data. The total Birks' uncertainty changes the response as much as 25% for candidates near threshold. Increasing the Birks suppression migrates events down in these distributions and some fall below the 1.5 MeV threshold. Modern studies of Birks' quenching in liquid and plastic scintillators [83,84], made for supernova and solar neutrino detection, are also at these energies and confirm the predicted scintillator response. This uncertainty is the largest single contribution in the first bin of Fig. 5, but is only 4%.
The GEANT4 cross section model uncertainties play a lesser role in this distribution than they do in the time and spatial distributions later. An increase in the cross section (decrease in the mean free path) makes candidates in and near the interaction exclusion region more likely and so fewer candidates. It contributes 4% to the uncertainty but only 2% in the first bin. The GEANT4 energy deposit model and photon yield are further explored using the benchmark modification shown with the solid black line.
Large rate uncertainties on the QE, 2p2h, and resonance models combine for 3% uncertainty on this distri-bution, but less than 1% uncertainty in the first two bins, and are the largest source for most bins for 0.4 < q 3 < 0.8 GeV/c. All these antineutrino processes produce some neutrons with similar energies. Distortions of the energy transfer spectrum for these processes, such as the uncertainty assigned to the RPA screening effect for QE used in the reference model [51,53] are at 4% and only important for the q 3 < 0.4 GeV/c subsample beyond the first two bins.
Effects related to the hadronic energy scale act in a special way and are significant for energy deposits above 10 MeV. They cause a migration of events up and down the range of q 3 . This migration effect is most significant and dominates parts of the error band for q 3 < 0.4. The subsample is lower statistics, has fewer candidates per event overall, and does not have compensating event feedin/out at its lower boundary. The hadronic energy-scale uncertainties are assigned based on test beam data with an enhanced uncertainty for the neutron response that was not directly tested. The GENIE FSI uncertainties also play this same migration role, in addition to directly changing the number of candidates in each event, and contribute a similar amount to total uncertainty.

B. Time-of-flight distribution
The time difference between the neutron candidate and the interaction time, t n −t 0 shown in Fig. 6 produces separation of the prompt electromagnetic component from the slowest neutron component. The simulation's overprediction for candidates with reconstructed energy <10 MeV in the upper panels appears roughly uniform across all high-statistics bins of time-of-flight. For higher energy candidates in the lower panels there are trends beyond the edge of the error band to relatively overpredict the latest times, either from neutrons that traveled the furthest or the slowest. The modified GEANT4 response model describes the data better than the modified GE-NIE model for the 0.4 < q 3 < 0.8 GeV/c top-right panel. By construction the GEANT4 benchmark has no effect on the lower panels, though the GEANT4 cross section uncertainties are significant in the error band. The GENIE modification reduces the slowest neutrons in the sample and has a slight shape effect smaller than the predicted error band.
Fluctuations to negative times are likely because the resolution on the timing for single-cluster candidates t n is 4.5 ns (shown later in Fig. 9): about one bin in these histograms. Systematic uncertainties directly from the measurement of time of flight for an individual event contribute negligibly. The simulation of the timing distribution is taken from a separate, in situ muon sample. The lack of bias is independently confirmed using clusters on the muon tracks of interactions in the selected sample and reconstructed the same as neutron clusters.
The same systematics described in the previous section are evaluated for this distribution. The GEANT4 neutron cross section model uncertainties enhance or reduce the appearance of neutrons that travel the furthest, and so have the longest times. It dominates the error bar in all bins beyond 15 ns. The other uncertainties described previously contribute roughly equal amounts in the center of the distribution.

C. Position upstream or downstream
The overprediction of candidates with energy deposits less than 10 MeV appear broadly around the interaction point, shown in Fig. 7. Far downstream the muon and accidental background dominate the E dep < 10 MeV sample, visualized as the difference between the total and the GENIE neutron lines, and the data are better described. For higher energy candidates, the simulation does well overall except for two underpredicted bins in the backward direction of the q 3 < 0.4 GeV/c panel (lower left).
There are more neutron candidates in the forward direction, where the QE process is especially relevant. In contrast, candidates from any process involving multiple particles can end up going backward from the interaction point. This includes multibody reactions 2p2h, ∆, and FSI in the nucleus, and neutrons produced when protons and pions and neutrons reinteract in the detector.
The cross section for neutron reinteractions is the direct and dominant uncertainty in the downstream region, as it was for long times. Again, the other uncertainties contribute roughly equally in the peak of the distribution.

D. Particle speed
We estimate the apparent speed of the particle as a fraction of the speed of light, β. Because the timing resolution plays a crucial role, and for better treatment of zero and negative times, it is more clear to present the apparent 1/β = speed of light × time/2D distance The result is shown in Fig. 8.
The estimate of a 2D distance can be made of the hypotenuse of the distance in Z shown Fig. 7 and the distance from the event axis (not the muon) along the one X, U, or V transverse direction measured for single-cluster candidates. If the neutron candidate is made of activity in more than one plane, the longest transverse position is used. This distance estimator is necessarily smaller than the true distance the neutron traveled, because it is missing the third of three coordinates, and because some neutrons bounce and take an indirect path to the point where an energy deposit is observed. This distribution has similar properties to the one in Fig. 7, and is not shown.
The systematic underestimate of the 2D distance means a systematic overestimate of 1/β of about 0.8 and an RMS resolution between 2 and 3, driven largely by the timing resolution. The resolution for the slowest particles with true 1/β > 5 is the worst because they do not travel very far and are observed closest to the interaction point. They have a resolution of around 4 and a bias of -0.8. The detector-only (without the effect of neutron multiple scatter) time and 2D distance resolutions are shown in Fig. 9. For neutrons, 1/β = 5 implies 20 MeV and 1/β = 10 implies 5 MeV, however the latter are expected to rarely produce candidates (see again Fig. 3) and the population beyond 1/β = 10 must be from fluctuations in the assigned time and distance. The resolutions and thresholds are such that the apparent 1/β is not usefully transformed into a kinetic energy distribution.
The GEANT4 cross section uncertainty, prominent in the time and z distance distributions separately, is much reduced and has little shape. A smaller or larger mean free path in the simulation affects both the time and the distance simultaneously. Other uncertainties contribute similarly across these distributions. The hadronic energy scale and FSI uncertainties (migrations in q 3 ) are especially significant in the first four bins of both the E cand > 10 MeV (lower plots), which is where discrepancies remain.
The electromagnetic component peaks near 1/β = 1.8. The neutron component peaks instead near 1/β ∼ 4 for candidates from the lowest kinetic energy neutrons and 1/β ∼ 2.5 for faster neutrons. The GENIE benchmark modification produces a reduction in the slowest neutron component that goes in the direction of the data, but again does not match the overall magnitude of the discrepancy in the upper right panel.

E. Neutron multiplicity
The overprediction of neutron candidates per event in the simulation also distorts the multiplicity of candidates per event, shown as percent of total in Fig. 10. The difference in percent (not the ratio) between the data and the reference simulation is shown in the first panel below the main distributions for compact comparison. The overprediction of neutrons masks the presence of the 2p2h component and RPA screening because both also enhance the percent of events with neutron candidates. Again the GENIE and GEANT4 serve as benchmark modifications showing the (data -modified simulations in the lower two rows, but now to expose these multinucleon features of the data. The two regions of q 3 shown in the previous figures are subdivided according to hadronic energy to produce distributions for QE-rich, dip, and ∆-rich subsamples, six in total. In the top panels the oversimulation of neutron candidates is most evident in the dip and ∆rich subsamples where the simulation significantly over predicts how many events have three or more neutron candidates and underpredicts how many have none.
The QE-rich subsets uniquely offer predicted sensitivity to the multinucleon effects of removing the 2p2h component and then the QE RPA screening effect leaving a GENIE model without either. One of the additional lines shows that removing the 2p2h component removes events that preferentially had multiple neutrons in the first place, increasing the bin with zero candidates. An-other line shows that removing also the QE RPA screening effect adds back events whose outgoing neutron was lower energy and less likely to make a candidate, also increasing the bin with zero candidates. In contrast to the QE-rich panels, there is no sensitivity to these multinucleon effects in the other panels; all reactions produce similar numbers of neutrons after FSI. Variations of RPA and 2p2h processes are the same size as the uncertainty bands in those panels.
A different way to summarize the subdivision of the data: neutrino model details in Fig. 10 are orthogonal to the neutron details in the previous figures. Modifications to the QE and 2p2h models show up in QE-rich region here while the excess of neutrons distorts all six panels similarly. The opposite happens in the previous figures; distortions of the spectra due to neutron production details are evident, but modifications to the 2p2h and QE models are largely flat with neutron candidate time, position, and speed.
What is desirable is to tune the neutron model to the dip and ∆-rich regions, a common technique when there are sidebands to a signal selection. Such a tune would correct and constrain the mis-modeled neutron effects in the multinucleon sensitive distributions. Though we do not directly have tunable parameters, a simplified version is obtained by remaking the distributions while applying the benchmark modifications to GEANT4 (third row) and GENIE (bottom row). The resulting dip and ∆-rich regions are now consistent with uncertainties for FIG. 9: Detector resolutions on the time and distance inputs to 1/β from the simulation, only for neutron candidates whose origin was a neutron from the GENIE simulation. There is no timing bias and the Gaussian fit to the timing resolution has σ = 4.5 ns. At the speed of light, the mean and RMS of the distance distribution correspond to 0.15 and 0.3 ns. Neutron multiple scattering effects are not included in the lower plot. 0 < q 3 < 0.4 GeV/c. The GEANT4 modification produces better distributions for 0.4 < q 3 < 0.8 GeV/c, perhaps overcorrecting, while the GENIE modification produces mild improvement that does not go far enough. Both roughly mimic the behavior of these benchmarks in the energy, time, position, and speed distributions.
In the QE-rich signal region, these benchmark modifications extrapolate the overall preponderance of neutrons in order to interpret the presence of multinucleon effects. Especially in the left-most two bins, the modified simulations now have an underprediction of events with one neutron candidate and an overprediction of events with none. In both cases, the resulting new predictions hint the data want even more 2p2h interactions or RPA screening than the reference MnvGENIE-v1.1 simulation. Alternatives are to add an additional process like deexcitation photons from carbon, or a more nuanced version of either of the benchmark modifications. In all the previous distributions, the current error bands preclude further detailed tests of the magnitude of RPA, the relative proton+neutron and neutron+neutron content of the 2p2h process, and the need for a low Q 2 suppression of resonances [85]. The sensitivity would be limited even if there were no large discrepancy, but modified RPA, 2p2h, or resonances would not explain the E dep < 10 sample.
Prior MINERvA measurements show distributions with sensitivity to the RPA and 2p2h multinucleon models that may also be sensitive to neutron effects. In Fig. 3 of [8] the reconstructed hadronic energy for this same sample is improved with the addition of RPA and a tuning of 2p2h to the neutrino data in [9], but the antineutrino agreement is not perfect. The untracked energy within 100 mm from the interaction point of antineutrino reactions in Fig. 25 of [68] is effectively in the excluded region of this analysis. That distribution is also not as well described by MnvGENIE-v1 compared to the equivalent Figs. 35-36 of [70]. In both examples, and based on the neutron observations in this paper, the description of the antineutrino data could be improved with reduction of the neutron component of the reconstructed energy in the simulation, while having a smaller effect on the neutrino mode data. Such a mechanism supposes the reduced neutron energy goes missing rather than being offset by additional charged hadron energy.

F. Protons
This final study is the analog to the untracked protons reported for the related neutrino interaction sample of [9], where the proton multiplicity with an untuned Valencia 2p2h process was at the edge of the error band of a prediction without it. Protons were counted by observing single strips with at least 20 MeV near the interaction point, indicating the Bragg peak at the end of the proton range. Because the single strip could be at the interaction point itself, the threshold is effectively just above 20 MeV. T2K has recently presented results for proton multiplicity [86] using a tracking threshold kinetic energy of 100 MeV Protons are also more common in the final state of a 2p2h antineutrino reaction, relative to antineutrino QE reactions. Repeating that earlier strategy reveals marginal sensitivity. The proton sample is more consistent with 2p2h in the QE-rich, zero proton bin (not shown, like the far-left bins of Fig. 10 but for proton activity), but the prediction without 2p2h is within the error band. This result, and the prominence of the uncertainty on the Birks' effect is the same as in the neutrino [9] case. In the dip region, the data are more consistent with 2p2h in one case and without 2p2h in the other case; the no-2p2h hypothesis is also just beyond the error band.
Liquid argon experiments [87][88][89][90] have shown more low-energy, charged proton tracks in the simulation compared to data. Under ideal circumstances, this detector technology permits tracking of protons with as little kinetic energy as 21 MeV. The GENIE model also produces more low-energy protons than other neutrino event generators, correlated with its behavior for neutrons. This supports that the GENIE benchmark modification may be part of resolving these discrepancies. In MINERvA, protons under 20 MeV would not meet the threshold for detection. They would usually deposit all their energy in the same scintillator strip the reaction occurred, and it takes 100 MeV before protons start to be trackable. So unlike neutrons and unlike protons in liquid argon detectors, multiple low-energy protons in a single scintillator strip in MINERvA would be counted only once, if at all.

VI. CONCLUSION
We have obtained the first time-of-flight, spatial, and speed (1/β) distributions of neutrons from antineutrino interactions. The reference simulation, whose components are widely used by neutrino oscillation experiments, overestimates the number of neutron candidates by 15% overall but by 25% for energy deposits less than 10 MeV, shown in the upper figures of Figs. 6-9. A possible interpretation is the GENIE neutrino event generator and its FSI model are overproducing the lowest energy neutrons. Also likely, the GEANT4 and detector models turn too many neutron interactions into measurable activity. Combinations and variations of these two benchmark modifications are paths forward. The discrepancy is around two standard deviations from the combination of the other sources of uncertainty.
Additional distortions may be present for candidates with energy deposits more than 10 MeV. The MC overestimates long time of flights relative to short in Fig. 6, far and forward relative to near and backward in Fig. 7, and slow relative to prompt in Fig. 8. These discrepancies are just beyond the error band, suggesting one or more of the uncertainties come close to accounting for these data.
It is a reasonable assumption that similar overproduction of small energy deposit neutron candidates is present for all subcomponents of the sample. In this case, the multiplicity distribution in the QE-rich subsamples indicate a preference for a model that has a combination of RPA screening and a 2p2h component, both of which reduce the relative proportion of events with zero neutron candidates.