We report the first double-differential neutrino-argon cross section measurement made simultaneously for final states with and without protons for the inclusive muon neutrino charged-current interaction channel. The proton kinematics of this channel are further explored with a differential cross section measurement as a function of the leading proton's kinetic energy that extends across the detection threshold. These measurements use data collected with the MicroBooNE detector from 6.4 × 10 20 protons on target from the Fermilab booster neutrino beam with a mean neutrino energy of ∼0.8 GeV. Extensive data-driven model validation utilizing the conditional constraint formalism is employed. This motivates enlarging the uncertainties with an empirical reweighting approach to minimize the possibility of extracting biased cross section results. The extracted nominal flux-averaged cross sections are compared to widely used event generator predictions revealing severe mismodeling of final states without protons for muon neutrino charged-current interactions, possibly from insufficient treatment of final state interactions. These measurements provide a wealth of new information useful for improving event generators which will enhance the sensitivity of precision measurements in neutrino experiments. DOI: 10.1103/PhysRevLett.133. 041801 Neutrino experiments that measure flavor oscillations as a function of neutrino energy aim to determine the neutrino mixing parameters and search for new physics beyond the standard model [1][2][3][4][5]. This requires precise mapping between reconstructed and true neutrino energy. The inclusive muon neutrino charged current (ν μ CC) interaction channel, ν μ N → μ -X, where N is the struck nucleus and X is the hadronic final state, is important for these measurements because it identifies the neutrino flavor with high purity and efficiency due to the sole requirement of detecting the muon.

A number of these neutrino experiments utilize liquid argon time projection chambers (LArTPCs) [5,6]. These tracking calorimeters have low detection thresholds and excellent particle identification (PID) capabilities [7][8][9][10][11]. LArTPCs enable the inclusive ν μ CC channel to be divided into subchannels based on the composition of the final state, each having a different mapping between true and reconstructed neutrino energy. This improves the energy reconstruction and increases the sensitivity of precision measurements [12].

Maximizing the physics reach of LArTPCs requires neutrino-argon interaction modeling capable of describing all final state particles. Existing models are unable to describe data with such detail, necessitating large interaction modeling uncertainties [12,13]. This is unsurprising; theoretical models attempting to describe experimental observables must simultaneously account for multiple scattering mechanisms [14], in-medium nuclear Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI. Funded by SCOAP 3 .

modifications to the fundamental neutrino interactions [15,16], and final-state interactions (FSI) involving the hadronic reaction products as they exit the nucleus [17]. The prominence of nuclear effects grows with the size of the target nucleus, further complicating the modeling of scattering for heavy nuclei like argon.

Efforts to simulate ν μ CC interactions benefit from measurements that simultaneously probe the leptonic and hadronic kinematics. Building on previous MicroBooNE work [18][19][20], and analogous to a similar measurement from T2K [21], we report a double-differential measurement of the muon energy, E μ , and muon scattering angle with respect to the neutrino beam, cos θ μ , for the ν μ CC channel split into final states with one or more protons ("Np" where N ≥ 1) and without protons ("0p"). An event is only included in the Np signal if the leading proton exiting the nucleus has kinetic energy above 35 MeV, which roughly corresponds to the proton tracking threshold in MicroBooNE [22]. The proton kinematics are further explored with a differential cross section measurement of the inclusive ν μ CC channel ("Xp" where X ≥ 0) as a function of the leading proton's kinetic energy, K p , that extends across the tracking threshold via the inclusion of a 0-35 MeV bin that includes events without a final state proton. A more expansive set of measurements employing the same analysis strategy can be found in [22].

Additional motivation for these measurements comes from the fact that LArTPCs utilize the gap between the neutrino and shower vertices to differentiate electrons from photons. The absence of additional vertex activity, usually from protons, makes it difficult to determine if a gap is present. This impacts ν e CC selections, which are the signal in many oscillation measurements, through lower efficiencies and purities for ν e CC 0p events than ν e CC Np events [23,24]. Since the ν μ CC channel is essential in constraining the ν e CC prediction, improved modeling of 0p and Np final states for ν μ CC is important. This need is highlighted by prior MicroBooNE results [20], which observed an excess of ν μ CC events at low reconstructed neutrino energies, potentially indicative of poor 0p cross section modeling.

We utilize data collected with the MicroBooNE detector [25] from an exposure of 6.4 × 10 20 protons on target (POT) from the booster neutrino beam (BNB) at a mean neutrino energy of ∼0.8 GeV [26]. The detector is comprised of a TPC volume with an active mass of 85 tonnes of liquid argon, and an array of 32 photomultiplier tubes (PMTs). When an interaction occurs in the detector, scintillation light and ionization electrons are produced by the charged particles emanating from the interaction. The light is recorded by the PMTs, providing ns-scale timing used to reject background cosmic ray events that are out of time with the beam. The ionization electrons drift in a 273 V=cm electric field to three wire readout planes which record charge distributions used for calorimetry and threedimensional (3D) mm-scale imaging.

Event reconstruction, calorimetry, and PID are performed with the Wire-Cell topographical 3D image processing algorithm [27]. Wire charge distributions are first deconvolved from the detector response by a TPC signal processing algorithm [28][29][30]. Wire-Cell uses the deconvolved readouts to reconstruct 3D images without topological assumptions about the source of activity as "tracks," which leave continuous energy depositions, or electromagnetic "showers" [31], which deposit more charge perpendicular to their trajectory. A many-tomany TPC-charge to PMT-light matching algorithm is used for further cosmic ray rejection [32].

Particle identification starts with finding kinks in the selected group of charge activity to identify tracks [31]. Candidate neutrino vertices and final state particles are identified concurrently based on dQ=dx, topology, and allowed particle relationships. A final neutrino vertex is chosen by a SparseConvNet deep neural network [33]. Proton and muon candidates are distinguished based on characteristic differences in their dQ=dx profile [22] using a test statistic constructed from a Kolmogorov-Smirnov shape comparison score and the normalization to the median dQ=dx of protons.

Two methods are used to calculate the energy of tracklike particles; range and summation of dE=dx. Between the two methods, the energy resolution is ∼10% for muons at all energies and ∼8% for protons below 200 MeV, above this the resolution degrades to ∼25% due to the increased probability of rescattering [22]. The range method is used for stopping tracks and is based on the NIST PSTAR database [34]. This method estimates the energy of tracks with minimal bias [22]. Summation of dE=dx is used for tracks that are shorter than 4 cm, exit the active volume, have a "wiggled" topology indicative of many small angle deflections [31], or emit δ rays. To calculate the kinetic energy, this method converts dQ=dx to dE=dx with an effective recombination model then sums dE=dx for each ∼6 mm segment of the track. The dE=dx method underestimates the energy by ∼10% [31,35], but this is incorporated into the detector model and appears consistent in data and simulation [20].

The ν μ CC event selection is identical to that of previous MicroBooNE work [18][19][20]22]. It utilizes the "generic neutrino selection" [36] as a preselection, which reduces cosmic ray contamination down to 15%. Backgrounds are further rejected with a boosted decision tree (BDT) trained using the XGBoost library [37] on a set of background taggers, which are variables designed to characterise nonν μ CC events. The efficiency (purity) of the fully inclusive ν μ CC selection is 68% (92%), with backgrounds predominantly coming from neutral current π AE events. The selection is further divided into 0p and Np selections based on the reconstructed leading proton's kinetic energy, K rec p . The Np selection contains events in which there is at least one proton with K rec p > 35 MeV. The 0p selection contains all other events. True Np events are analogously defined as PHYSICAL REVIEW LETTERS 133, 041801 (2024) 041801-3 having a proton with true K p > 35 MeV and true 0p events are defined as having either zero final-state protons, or no proton with K p > 35 MeV. The Np selection has a 49% efficiency for true Np events and a high purity of 95% due to low contamination from non-ν μ CC and 0p events. The 0p selection has an efficiency of 53% for true 0p events. The larger number of Np events, which outnumber 0p events ∼7∶1, increases the prominence of Np events in the 0p selection, reducing the fraction of true 0p events in the 0p selection to 32% [22].

The 0p and Np cross sections are extracted simultaneously. This allows the number of true Np events in the 0p selection to be predicted based on the observation of the Np selection (and vice versa). The same strategy is employed in other MicroBooNE work [22,24]. Simultaneous extraction requires unfolding

where M is the reconstructed distribution, S is the differential cross section to be extracted, B is the distribution of background events that are not part of the inclusive ν μ CC channel, and R is the response matrix describing the mapping between the true and reconstructed distributions.

The first (second) index on R corresponds to the reconstructed (true) proton multiplicity. The Wiener-SVD unfolding technique [38] is used to extract nominal flux-averaged cross section results [39]. This method returns a regularized unfolded cross section and corresponding covariance matrix, V S , which describes the uncertainties and bin-to-bin correlations on the result. The form of Eq. ( 1) allows the unfolding to account for correlations between the 0p and Np channels [22]. An additional smearing matrix, A C , that captures the bias induced by regularization is also obtained in the unfolding. Cross section predictions should be multiplied by A C when compared to the unfolded result. The extracted cross sections, A C , and V S which is obtained via blockwise unfolding to preserve intervariable correlations [22,40], can be found in the Supplemental Material [41].

Monte Carlo (MC) simulations are used to estimate R and part of B. The neutrino flux is modeled with the GEANT4 simulation of the BNB from MiniBooNE [26,42]. Neutrino-argon interactions are simulated with the G18_10a_02_11a configuration of the GENIE v3.0.6 event generator [43] tuned to ν μ CC data without final state pions from T2K [44] by reweighting based on two CC quasielastic and CC meson-exchange-current parameters [45]. The resulting prediction is referred to as the "MicroBooNE Tune." Final state particles are propagated through a detector simulation using the GEANT4 toolkit v4_10_3_03c [42] and processed using the LArSoft [46] framework. The simulated TPC and PMT waveforms are overlaid on beam-off data to provide an accurate description of cosmic ray activity. These overlaid MC samples are processed like real data and used to estimate R and B.

Uncertainties are estimated with covariance matrices calculated from a multi-universe approach as in [22]. The total covariance matrix,

given by the sum of the covariance matrices calculated for each systematic uncertainty described below.

Uncertainty in the neutrino flux [26] is contained in V flux . The flux contributes 5%-10% uncertainty to the cross section results and is often the dominant systematic for the Np cross section measurements. Neutrino-argon cross section uncertainties are accounted for in V xs [45] and contribute ∼5% uncertainty. In V reint , uncertainties on reinteractions of final state particles outside the nucleus are accounted for. These are estimated using the GEANT4REWEIGHT [47] package and are relatively small except at high K p where reinteractions occur for > 75% of protons. The flux, cross section, and reinteraction uncertainties are estimated with the multi-sim technique [48].

Detector response uncertainties [49] are contained in V det . These are the biggest sources of uncertainty for 0p events, typically ranging from 5%-15% compared to ∼5% for Np, and are larger at high energies and backward scattering angles. To evaluate these uncertainties, a detector model parameter is varied by 1σ and bootstrapping is used to estimate the impact of this variation and form V det [20,22]. A Gaussian process regression smoothing algorithm [50][51][52] is implemented to prevent an overestimation of detector systematics due to statistical fluctuations [19,22].

Flat 50%, 2%, and 1% uncertainties on neutrino interactions outside the detector, POT counting, and the number of target nuclei are contained in V dirt , V POT , and V Target , respectively. Their impact on the extracted cross sections is small.

A data-driven model validation procedure is employed to detect mismodeling that may bias the extracted cross section results. As described in [22], this relies on the conditional constraint formalism [53] to increase the stringency of the validation. These constraints leverage correlations between variables and channels arising from shared physics modeling to update the model prediction and reduce its uncertainties based on data observations. The unfolding does not utilize these constraints; they are only for model validation. To be validated, the model is required to describe the data at the 2σ level. This is evaluated with χ 2 goodness of fit (GOF) tests interpreted by using the number of degrees of freedom, ndf, which corresponds to the number of bins, to obtain p-values.

A kinematic variable that is relatively well understood and reconstructed can be validated by directly comparing the model prediction to the data in the phase space relevant to the unfolding. The muon kinematics for events fully contained (FC) within the detector fit this criteria, and are validated with a GOF test on the E rec μ distributions in cos θ rec μ slices used for the cross section extraction. These tests yield p-values of 0.45 and 0.98 for the 0p and Np distributions, respectively. The model passes validation in these tests. For the partially contained (PC) muon kinematics, the modeling of activity outside the active detector volume, which cannot be reconstructed, must also be validated. This is done by using the 0p and Np FC muon kinematics to constrain the PC distributions. These tests, found in the Supplemental Material of [22], result in pvalues of 0.84 for 0p and 0.99 for Np. This indicates that the overall model adequately describes the 0p and Np E rec μ distributions in cos θ rec μ slices for PC events. When the K rec p distribution is constrained by the 0p and Np muon kinematics the resulting p-value is 3 × 10 -5 , indicating that the model is unable to describe the observed proton kinematics. The discrepancy appears at low K rec p where modeling becomes challenging due to the prominence of FSI and nuclear effects. Appropriate uncertainties on K rec p are important for the K p differential cross section measurement and the split into 0p and Np subchannels. Thus, this shortcoming of the model motivates expanding the uncertainties to mitigate the possibility of extracting biased results. The additional uncertainty was obtained empirically by estimating a true K p distribution for FC signal events passing the generic neutrino selection by unfolding the FC K rec p data using the constrained prediction and only statistical uncertainties. The ratio of the unfolded distribution to the constrained signal prediction defines a reweighting function. This is applied to all ν μ CC signal events and is treated as a 1σ deviation from the original prediction due to a cross section effect allowing V rw to be calculated identically to V xs [22]. Including V rw reduces the p-value for the K rec p distribution to 0.82 after constraint, enabling the extraction of the desired cross sections.

The extracted nominal flux-averaged cross sections are compared to event generator predictions from GENIEv3.0.6G18_10a_02_11a (GENIE) [43], the MicroBooNE tune GENIE configuration (μBooNE tune) [45], NuWro21.02 (NuWro) [54], GiBUU 2023 (GiBUU) [55], and NEUT 5.4.0.1 (NEUT) [56]. These were processed with NUISANCE [57], do not include theoretical uncertainties, and are smeared with A C obtained from unfolding. Agreement between the data and each prediction is quantified by χ 2 =ndf values calculated with uncertainties according to V S and ndf corresponding to the number of bins.

The double-differential cos θ μ and E μ cross section results are shown in Fig. 1 as a function of bin index, which are in angular slices ranging from backwards on the left to forward on the right and increasing in energy along each slice. Binning details are in the Supplemental Material [41]. In all three channels, the generators tend to underpredict the peak of the E μ distribution, with the only exception being for GiBUU in the 0p channel, where its prediction shows good normalization agreement around the peak.

Figure 1(a) shows the fully inclusive Xp result. The χ 2 values indicate that NEUT best describes this data. The 0p result is shown in Fig. 1(b). At forward angles, GiBUU describes this result better than the other generators, which significantly underpredict the cross section. Agreement at backwards angles is more comparable, but GiBUU still has the lowest χ 2 =ndf over all bins. NEUT, though offering the best description of the inclusive channel, shows the largest discrepancy with the 0p data. The differential cross section for the inclusive channel as a function of K p is shown in Fig. 2. The first bin extends from 0-35 MeV and includes events without a final state proton; it is equivalent to the 0p signal definition. The generator predictions diverge at low energies, particularly, for the 0p bin where only GiBUU is able to describe the data. This gives GiBUU the lowest χ 2 despite its underprediction of the data at moderate-to-high energies. Similar underprediction is seen for GENIE and NuWro in this region. NEUT and the μBooNE tune describe moderate-tohigh energies well, but NEUT also significantly underpredicts the 0p bin, increasing its χ 2 .

Comparing the 0p and Np results to the Xp result in Fig. 1 demonstrates how a model that does well for inclusive scattering may not also be able to describe the hadronic final state. NEUT stands out in its relatively good description of Xp but not of 0p due to a significant underprediction of the 0p cross section. This is possibly attributable to the way NEUT treats binding energy for nucleon FSI. NEUT assigns nucleons an effective mass when propagating them through the nucleus and only allows interactions if the total energy is twice the energy of the free nucleon mass [56]. This reduces the strength of FSI for low energy nucleons [17] leading to a sharp dropoff in the cross section at low K p and the low 0p cross section prediction that agrees poorly with the data.

Unlike NEUT, GiBUU describes the data relatively well when the ν μ CC channel is split into final states with and without protons. This can possibly be attributed to its implementation of FSI with a transport model, which includes a binding potential that treats "target" and "ejected" nucleons identically. The "ejected" nucleons are propagated on realistic trajectories through the residual nucleus according to a potential consistent with the initial interaction [55,58]. This is absent in other generators, which propagate ejected nucleons on straight lines and decouple the initial interaction from the FSI. After the initial interaction, an outgoing nucleon may repeatedly collide with other nucleons, depleting the initial nucleon of its energy and shifting the K p distribution towards smaller values [17,[58][59][60][61]. It is plausible that the more selfconsistent treatment of FSI in GiBUU better captures this effect. This hypothesis is consistent with the way GiBUU better describes the 0p final states in Fig. 1, and has the only prediction that mirrors the sharp peak in the data at the lowest K p in Fig. 2. Additional 0p and Np results in [22] and measurements of transverse kinematic imbalance variables [62] sensitive to FSI modeling in [63,64] support this hypothesis. Of course, different modeling of the initial neutrino-nucleon interaction could also play a role, especially given that none of the generators adequately describe the inclusive measurement shown in Fig. 1(a).

In summary, we report differential cross-section measurements of the ν μ CC channel that probe the phase space of lepton and hadronic kinematics. An underprediction of the cross section for final states without protons is observed for all event generator predictions except GiBUU, which offers a significantly better description of the data possibly due to its more sophisticated treatment of final state interactions. These measurements provide new information to stimulate further improvement of models and generators to match the precision required for future neutrino oscillation measurements and beyond the standard model searches.

Acknowledgments-This document was prepared by the MicroBooNE Collaboration using the resources of the