Multi-Differential Cross Section Measurements of Muon-Neutrino-Argon Quasielastic-like Reactions with the MicroBooNE Detector

We report on a flux-integrated multi-differential measurement of charged-current muon neutrino scattering on argon with one muon and one proton in the final state using the Booster Neutrino Beam and MicroBooNE detector at Fermi National Accelerator Laboratory. The data are studied as a function of various kinematic imbalance variables and of a neutrino energy estimator, and are compared to a number of event generator predictions. We find that the measured cross sections in different phase-space regions are sensitive to nuclear effects. Our results provide precision data to test and improve the neutrino-nucleus interaction models needed to perform high-accuracy oscillation analyses. Specific regions of phase-space are identified where further model refinements are most needed.


I. INTRODUCTION
High-precision measurements of the neutrino mixing angles, mass differences, and charge-parity violating phase, and the search for physics beyond the Standard Model are the primary physics goals of many currently operating as well as next-generation neutrino experiments [1][2][3][4][5][6]. These measurements require reliable comparisons of measured and theoretically-expected neutrino interaction rates in the corresponding detectors. Thus, understanding the neutrino-nucleus scattering processes in detail is a prerequisite for these experiments to reach their discovery potential. A number of neutrino oscillation experiments employ liquid argon time projection chambers (LArTPCs) [3][4][5][7][8][9] to detect the particles produced in neutrino interactions. The ultimate goal of these efforts is both to reconstruct the energy of the neutrino based on the kinematics of the outgoing particles and to enable few-percent-level modeling of neutrino-argon interaction rates [10]. Therefore, highaccuracy modeling of neutrino-argon interactions is of the utmost importance [11][12][13].
This work presents the first measurement of fluxintegrated single-and double-differential cross sections for muon-neutrino-argon (ν µ -Ar) charged-current (CC) quasielastic (QE)-like scattering reactions as a function of kinematic imbalance variables [14][15][16][17][18]. Doubledifferential measurements as a function of a neutrino energy estimator are further reported for the first time in kinematic imbalance bins on argon. Motivated by a previous analysis with a similar signal event topology [19], we focus on reactions where a single muon-proton pair is reconstructed with no additional detected particles. The results reported here use the MicroBooNE detector [20] * microboone info@fnal.gov with an exposure of 6.79 × 10 20 protons on target from the Booster Neutrino Beam (BNB) [21] at Fermi National Accelerator Laboratory.
The experimental setup is presented in Sec. II, followed by the signal definition and event selection in Sec. III. The observables of interest are defined in Sec. IV. Section V describes the cross section extraction and systematics procedure and Sec. VI outlines the modeling configurations used for comparison to the data. The results are reported in Sec. VII and the conclusions are discussed in Sec. VIII.

II. EXPERIMENTAL SETUP
The MicroBooNE LArTPC has an active volume that contains 85 tonnes of argon. It is exposed to BNB neutrinos, with an energy spectrum that peaks around 0.8 GeV and extends to 2 GeV.
Charged particles are produced after the primary neutrino interaction with the argon nuclei in the LArTPC active volume. Scintillation light and electron ionization trails are produced while these charged particles travel through the liquid argon. In the presence of an electric field of 273 V/cm, the ionization electrons drift towards a system of three anode wire planes. Photomultiplier tubes (PMTs) are used to measure the scintillation light.
If the PMT signals are in time coincidence with the beam arrival time, then events are recorded. Trigger hardware and software selection criteria are designed to minimize the contribution from background events, which are primarily cosmic muons. After these are applied, enriched data samples are obtained in which a neutrino interaction occurs in ≈ 15% of selected beam spills [22].
Individual particle tracks are reconstructed with Pandora pattern recognition algorithms based on the measured ionization signals in the enriched data sam-ples [23]. Particles are identified based on the measured track energy deposition profile, while the particle momenta are obtained based on the track length [24,25].

III. SIGNAL DEFINITION & EVENT SELECTION
The QE-like signal definition used in this analysis includes all ν µ -Ar scattering events with a final-state muon with momentum 0.1 < p µ < 1.2 GeV/c, and exactly one proton with 0.3 < p p < 1 GeV/c. Events with finalstate neutral pions at any momentum are excluded. Signal events may contain any number of protons below 300 MeV/c or above 1 GeV/c, neutrons at any momentum, and charged pions with momentum lower than 70 MeV/c. We refer to the events passing this definition as CC1p0π. The aforementioned momentum ranges are driven by considering resolution effects, as well as regions of the phase space with non-zero efficiencies and systematic uncertainties that are well understood. This signal consists predominantly of QE events. More complex interactions as labeled at a generator level, namely meson exchange currents (MEC), resonance interactions (RES) and deep inelastic scattering events (DIS), can still yield CC1p0π events. That can be the case due to final-state interactions (FSI), such as pion absorption, or due to the presence of particles outside the momentum range of interest in the CC1p0π signal definition as defined above.
Events that satisfy the CC1p0π signal definition at a reconstruction level but not at truth level are treated as background events. We refer to these events as non-CC1p0π. Based on simulation predictions, we find that the dominant background contribution originates from events with two protons in the momentum range of interest, where the second proton was not reconstructed. These events are referred to as CC1µ2p0π and are the focus of a dedicated MicroBooNE cross section analysis that demonstrated good data-simulation agreement [26].
Candidate muon-proton pairs are isolated by requiring the existence of precisely two track-like and no showerlike objects, as classified by Pandora using a track-score variable [27,28]. The log-likelihood ratio (LLR) particle identification (PID) score [29] is used to identify the muon and proton candidates. Figure 1 shows the particle composition breakdown of the sample as a function of the LLR PID score.
Muons tend to have higher LLR PID score values than protons, thus the track with the highest score is tagged as the candidate muon. Meanwhile, the track with the lower score is treated as the candidate proton.
Cosmic muon and non-CC1p0π contamination backgrounds were significantly reduced by applying a requirement on the candidate proton LLR PID score. We studied the effect of cutting on different values of this quantity, which has a strong discrimination power for rejecting MC non-CC1p0π background, out-of-cryostat and cosmic events. That yielded an optimal cut on the pro- ton candidate LLR score of < 0.05, as shown in Fig. 2a. Figure 2b shows the corresponding muon candidate LLR score, which is peaked at values close to one. The uncertainty bands account for potential data-MC discrepancies observed for both particle scores. The particle composition of the panels included in Fig. 2 is shown in the Supplemental Material.
To further minimize the contribution of misreconstructed track directions, we took advantage of two muon momentum reconstruction methods available for contained tracks, namely the momentum from range [30] and the momentum from Multiple Coulomb Scattering (MCS) [31]. The range and MCS muon momenta needed to be in agreement within 25% and the improvement in the muon momentum reconstruction can be seen in Fig. 3. We required that the distance between the track start points and the vertex is smaller than the corresponding distance between the track end points and the vertex. We also demanded that the distance between the start points of the two candidate tracks is smaller than the distance between the two end points.
Further reduction of the cosmic tracks and minimization of bin-migration effects is achieved by considering only fully contained candidate muon-proton pairs within a fiducial volume of 10 cm inside the edge of the detector active volume. We retain 9051 data events that satisfy all event selection criteria.
In order to provide an accurate description of the dominant cosmic backgrounds pertinent to surface detectors, the full Monte Carlo (MC) simulation consists of a combination of simulated neutrino interactions overlaid on top of beam-off background data. This approach has been extensively used by MicroBooNE [19,[32][33][34]. The GENIE v3.0.6 event generator is used to simulate neutrino interactions with the G18 10a 02 11a configuration [35,36]. The CCQE and CCMEC predictions have been additionally tuned to T2K ν µ -carbon CC0π data with any number of protons in the final state [37,38]. The different target nuclei across T2K and MicroBooNE might result in particle reinteraction differences that can affect the reconstructed final state topologies, such as different absorption effects. Yet, the T2K data sets used for tuning are dominated by CCQE and CCMEC interaction processes, which are the main contributors to the CC1p0π topology presented in this work. Predictions for more complex interactions, such as RES, remain unaltered and no additional MC constraints are applied. We refer to the corresponding tuned prediction as G18. All the final state particles following the primary neutrino interaction are generated by GENIE. They are further propagated in GENIE through the nucleus to account for FSI. The propagation of the particles outside the nucleus is simulated using GEANT4 [39]. The MicroBooNE detector response is modeled using the LArSoft framework [40,41]. Based on this MC prediction, we obtain a purity of ≈ 70% and an efficiency for selecting CC1p0π events of ≈ 10%. The final efficiency is primarily driven by the demand for exactly two fully contained track-like candidates.

IV. OBSERVABLES
In neutrino-nucleus scattering events, there is an imbalance between the true initial neutrino momentum and the true sum of final-state lepton and hadron momenta as a result of nuclear effects [14]. A schematic representation of the kinematic imbalance variables of interest in this work is shown in Fig. 4.
Using the CC1p0π candidate muon-proton pair kinematics, the missing momentum in the plane transverse to the beam direction is defined as where p T µ and p T p are the projections of the momenta of the outgoing lepton and proton on the transverse plane, respectively. In the absence of nuclear effects, purely QE interactions would yield δp T = 0. In the presence of the dense nuclear medium, this variable encapsulates information related to the Fermi motion, but it is smeared due to FSI and non-QE interactions, as can be seen in Fig

FIG. 5.
Distribution of the selected CC1p0π events as a function of the transverse missing momentum δpT . Only statistical uncertainties are shown on the data. The interaction contributions are obtained from simulation and their separation in signal (S) and background (B) events is presented. The bottom panel shows the ratio of data to prediction.
The direction of the transverse momentum imbalance δp T is described by the angle which is uniformly distributed in the absence of FSI due to the isotropic nature of the Fermi motion. In the presence of FSI, the proton momentum is generally reduced and the δα T distribution becomes weighted towards 180 • , as can be seen in Fig. 6. The opening angle δφ T between the correlated candidate muon-proton pair on the transverse plane is given Distribution of the selected CC1p0π events as a function of the transverse missing momentum direction δαT . Only statistical uncertainties are shown on the data. The interaction contributions are obtained from simulation and their separation in signal (S) and background (B) events is presented. The bottom panel shows the ratio of data to prediction. by In the absence of nuclear effects, QE events would be concentrated at δφ T = 0. When nuclear effects are present, QE events can occupy wider angles. At the same time, non-QE events are dominant in the high δφ T part of the tail and their contribution is fairly flat across all angles, as can be seen in Fig. 7. The muon-proton momentum imbalances transverse and longitudinal to the transverse lepton momentum [17] are defined as and can also be written as These distributions can be seen in Fig. 8 and Fig. 9, respectively. The δp T,x distribution is symmetric around 0 GeV/c due to the presence of the sin δα T factor in Eq. 5 and the fact that δα T ranges from 0 o to 180 o . The width of the distribution is driven by the Fermi motion that affects the δp T magnitude. Unlike δp T,x , the δp T,y distribution is asymmetric with an enhanced contribution from negative values. The asymmetry is driven by the presence of the cos δα T factor in Eq. 5 and the fact that Distribution of the selected CC1p0π events as a function of the muon-proton transverse opening angle δφT . Only statistical uncertainties are shown on the data. The interaction contributions are obtained from simulation and their separation in signal (S) and background (B) events is presented. The bottom panel shows the ratio of data to prediction.
δα T is mainly peaked around 180 o . Given that the forward δα T peak is driven by FSI, the size of the δp T,y asymmetry is also sensitive to the FSI strength. Finally, the calorimetric energy reconstruction is investigated, where E µ is the muon energy, T p is the  proton kinetic energy and BE = 0.04 GeV is the average binding energy for argon [42]. This energy estimator, shown in Fig. 10, is an approximation for the true energy of the incoming neutrino and is used in oscillation searches.

V. CROSS SECTION EXTRACTION & SYSTEMATICS
The flux-averaged differential event rate as a function of a given variable x in bin i is obtained by where N i and B i are the number of measured events and the expected background events, respectively. T is the number of target argon nuclei in the fiducial volume of interest. Φ ν corresponds to the integrated BNB flux and ∆ i corresponds to the i-th bin width or area for the singleand double-differential results, respectively. We report the extracted cross sections for CC1p0π interactions using the Wiener singular value decomposition (Wiener-SVD) unfolding technique as a function of unfolded kinematic variables [43]. This unfolding procedure corrects a measured event rate for inefficiency and resolution effects. This is achieved by performing a minimization of a χ 2 score that compares data to a prediction and allows for a regularization term. A Wiener filter determines the level of regularization that is required to minimize the mean square error between the variance and bias of the result. In addition to the measured event rate, the method uses a covariance matrix calculated from simulated events accounting for the statistical and systematic uncertainties on the measurement as input. It also requires the construction of a response matrix describing the expected detector smearing and reconstruction efficiency.
The output of the method is an unfolded differential cross section, a covariance matrix describing the total uncertainty on the unfolded result, and an additional smearing matrix that we refer to as A C . The latter contains information about the regularization and bias of the measurement. The corresponding A C matrices have been applied to all the cross section predictions included in this work when a comparison to the unfolded data is performed. The A C matrix should be applied to any independent theoretical prediction when a comparison is performed to the data reported in this paper. The data release, the unfolded covariance matrices, and the additional matrices A C can be found in the Supplemental Material.
The total covariance matrix E ij = E stat ij + E syst ij includes the statistical and systematic uncertainties on the differential event rate associated with our measurement. E stat ij is a diagonal covariance matrix with the statistical uncertainties and E syst ij is a covariance matrix that incorporates the total systematic uncertainties both on the CC1p0π signal and on the non-CC1p0π background events as detailed below.
The neutrino flux is predicted using the flux simulation of the MiniBooNE collaboration that used the same beam line [44]. Neutrino cross section modeling uncertainties were estimated using the GENIE framework of event reweighting [35,36,38]. The rescattering uncertainties were obtained using GEANT4 and the relevant reweighting package [45]. For each of these sources of uncertainty, we use a multisim technique [46], which consists of generating a large number of MC replicas, each one called a "universe", where model parameters are varied within their uncertainties. The simultaneous varying of many model parameters provides a correct treatment of their correlations. A total of n such universes are used to construct a covariance matrix corresponding to each source of uncertainty, where R CV i (R CV j ) and R k i (R k j ) are the flux-averaged event rates for the central value and systematic universe k in a measured bin i (j), respectively. The resulting covariance matrices are summed together to estimate the relevant uncertainty from each source.
In order to account for potential biases due to the nominal MC modeling prediction used in the unfolding procedure and presented in Sec. VI, an additional cross section uncertainty using the NuWro v19.02.2 event generator prediction [47] as an alternative universe has been added. The relevant NuWro modeling is significantly different when compared to the nominal MC one, as detailed in Sect. VI. The flux-integrated NuWro cross sections are obtained using Eq. 7 and the corresponding covariance matrices are constructed using Eq. 8 and a single universe (n = 1).
For detector model systematic uncertainties, one detector parameter is varied each time by 1σ and is referred to as a "unisim". These include variations in the light yield, the ionization electron recombination model, space-charge effects, and waveform deconvolution [48]. We then examine the impact of each parameter variation on the MC event rates by obtaining the differences with respect to the central value on a bin-by-bin basis. We define the total detector 1σ systematic uncertainty by summing in quadrature the effect of m detector variations using the formalism introduced in Eq. 8, The full fractional uncertainty on the integrated total cross section is 11% and includes contributions from the neutrino flux prediction (7.3%), neutrino interaction cross section modeling (6%), detector response modeling (4.9%), beam exposure (2.3%), statistics (1.5%), number-of-scattering-targets (1.2%), reinteractions (1%), and out-of-cryostat (dirt) interaction modeling (0.2%). The Supplemental Material includes tables detailing all the cross section uncertainties used in this work. The main contributors are found to be the strength of the RPA correction and CCMEC cross section shape. The signal related cross section uncertainties are found to be 8.6%, while the background ones account for 6.3%. Note that the individual contributions are higher than the total cross section uncertainty of 6% due to correlations between the signal and background events, since the same interaction processes can contribute both as signal and background.
In the results presented below, the inner error bars on the reported cross sections correspond to the statistical uncertainties. The systematic uncertainties were decomposed into shape-and normalization-related sources following the procedure outlined in [49]. The crossterm uncertainties were incorporated in the normalization part. The outer error bars on the reported cross sections correspond to statistical and shape uncertainties added in quadrature. The normalization uncertainties are presented with the gray band at the bottom of each plot. Overflow (underflow) values are included in the last (first) bin. The relevant A C matrices have been applied to the theoretical predictions to account for regularization effects.
In addition to the alternative event generators, our results are compared to a number of different GENIE configurations.
Finally, the newly added theory-driven GENIE v3.2.0 G21 11b 00 000 configuration (G21) is shown. This includes the SuSAv2 prediction for the QE and MEC scattering parts [69] and the hN2018 FSI model [70]. The modeling options for RES, DIS, and COH interactions are the same as for G18.
To quantify the data-simulation agreement, the χ 2 /bins ratio data comparison for each generator is shown on all the figures and is calculated by taking into account the total covariance matrix. Ratios close to unity are indicative of a sufficiently accurate modeling performance. Theoretical uncertainties on the models themselves are not included.

VII. RESULTS
Along with the aforementioned kinematic imbalance and energy estimator results, the data are also presented as a function of the lepton angular orientation (Fig. 11). Previous MicroBooNE measurements using different signal definitions [19,71,72] showed discrepancies in that quantity, primarily in the forward direction. These analyses used an older simulation prediction, namely GENIE v2.12.2, to account for the efficiency corrections and beam-induced backgrounds. This work illustrates that all generator (Fig. 11a) and GENIE configuration ( Fig. 11b) predictions are in good agreement with the data when reported as a function of cosθ µ . Figures 12 and 13 show the measured single-differential cross sections as a function of δp T using all the events (panel a), as well as the double-differential results as a function of the same kinematic variable in δα T bins (panels b-e). In the presence of FSI, the proton can rescatter or be absorbed, yielding larger kinematic imbalances on the transverse plane and δp T values that extend beyond the Fermi momentum, as can be seen in Fig. 14. Furthermore, the same extended tail can be obtained when pions produced due to multi-nucleon effects (MEC or RES) are either absorbed or below the detection threshold. The single-differential result shows such a high-momentum tail that extends above 0.8 GeV/c. This picture is consistent with the results reported by the T2K and MINERvA collaborations [15,16,73]. Unlike the single-differential result, the double differential results with low δα T extend only slightly above 0.4 GeV/c. That indicates that this region contains minimal FSI and multi-nucleon effects and the δp T distribution is driven by the nucleon Fermi motion. On the other hand, the higher δα T values correspond to δp T distributions that extend beyond 0.8 GeV/c. This behavior is indicative of the presence of FSI and multi-nucleon effects that smear the δp T distribution to higher values. Future multi-differential results can help further disentangle the contributions from these effects. Figure 12 shows the comparisons to a number of available neutrino event generators with NuWro and G18 showing the best agreement over all events. Figure 13 shows the same results compared to a number of GENIE configurations illustrating that Gv2 is disfavored, an observation that is driven by the Gv2 low δp T behavior. Furthermore, Untuned shows a good χ 2 /bins performance across all slices but predicts lower values than data. Additionally, Fig. 14 shows the effect of final state interactions (FSI) on the CC1p0π selection using the G18 configuration of GENIE. The addition of FSI allows for more non-QE events to satisfy the CC1p0π signal definition that smear the δp T distribution to higher values. Figure 15 shows the double-differential results as a function of δp T in cosθ µ bins. In a factorized nuclear model such as the LFG, the Fermi motion part of δp T should stay constant in terms of the shape as a function of the outgoing lepton kinematics, since in such models the initial state nucleon momentum is a property of the nucleus that cannot be affected by the interaction momentum or energy transfer. That is indeed the observed behavior in the reported results across all event generators and configurations, where no evidence of the inadequacy of the factorization approach is observed. Figure 15 shows the comparisons to a number of available neutrino event generators, where the G18 prediction is favored based on the χ 2 /ndf results. Apart from the factorization, a better separation between QE and non-QE can be gained depending on the cosθ µ region. As can be seen in Fig. 16 for G18, MEC events play a more pronounced role for forward muon scattering and in the high δp T tail, as opposed to backward scattering angles, which are much more strongly populated by QE events. Furthermore, the G18 cross section prediction falls below the data in the -1 < cosθ µ < 0.5 region, as seen in Fig. 16a and Fig. 16b. That could indicate that additional contribution from the QE part of the G18 prediction is needed beyond the MicroBooNE tune. Figure 17 shows the same interaction breakdown for GiBUU. Unlike G18, GiBUU illustrates a peak shift to the right, which becomes more pronounced in the backward direction. This shift is driven by the enhanced MEC contribution in higher δp T values and the reduced QE contribution at smaller values. In the backward direction, GiBUU further shows a cross section excess driven by the MEC contribution. Figure 18 shows the same results compared to a number of GENIE configurations illustrating that Gv2 is disfavored due to the low δp T bin behavior. Figures 19 and 20 show the double-differential cross section as a function of δp T in cosθ p bins. The factorization of the nuclear motion is mostly preserved in cosθ p bins, analogously to the previous result in cosθ µ . Figure 19 shows the comparisons to a number of available neutrino event generators. The GiBUU prediction is significantly lower than the data in the backward proton angle for low δp T values, as shown in Fig. 19a. Figure 20 shows the same results compared to a number of GENIE configurations illustrating that Gv2 is disfavored across all cosθ p bins. As can be seen in Fig. 21, this particularly poor performance is driven by the QE contribution. For backward scattering events (panel a), the QE contribution predicted by Llewellyn Smith is significantly overestimated. For intermediate angles (0 < cosθ p < 0.5), the same QE model results in an unphysical double peak. For forward scattering (0.5 < cosθ p < 1), the Gv2 QE prediction yields a pronounced contribution at lower values of δp T compared to the data. tions as a function of δα T using all the events (panel a), as well as the double-differential results in the same kinematic variable in δp T bins (panels b-d). The single-differential results shown in panel a yield some interesting observations when compared to the relevant T2K and MINERvA results [15,16,73]. Our distribution illus- trates a slightly asymmetric behavior, similar to the one reported by the T2K collaboration at a comparable energy with MicroBooNE. Both the already-published T2K results on carbon and the ones presented in this work on argon mostly demonstrate data-MC agreement within the experimental uncertainties. Therefore, the mass- Shape 16. Comparison between the flux-integrated double-differential cross sections as a function of δpT for data and the G18 GENIE prediction in cosθµ bins. Inner and outer error bars show the statistical and total (statistical and shape systematic) uncertainty at the 1σ, or 68%, confidence level. The gray band shows the normalization systematic uncertainty. Colored stacked histograms show the results of theoretical cross section calculations using the G18 GENIE prediction for QE (blue), MEC (orange), RES (green), and DIS (red) interactions. ment by MINERvA reports a more pronounced asymmetry on hydrocarbon. The breakdown plots in Fig. 18 in Ref. [73] show that this behavior is driven by enhanced pion-production rates due to the higher average beam energy. Low δp T values result in a fairly uniform δα T distribution indicative of the absence of FSI effects in that part of the phase-space. On the other hand, higher δp T values result in a highly asymmetric δα T distribution, which is driven by the strength of the FSI interactions. Figure 22 shows the comparisons to a number of available neutrino event generators, where NuWro is the generator with the most conservative FSI strength. Figure 23 shows the same results compared to a number of GENIE configurations, where Gv2 yields the highest χ 2 /bins result, especially in the lowest δp T region. As shown in Fig. 24, this is driven by the Gv2 QE performance, which results in peaks at the edges of the distribution, unlike the data result. Additionally, Fig. 25 shows the effect of FSI on the CC1p0π selection using the G18 configura-tion of GENIE that introduces an asymmetric behavior in δα T . Figures 26 and 27 show the double-differential results as a function of δα T in cosθ µ bins. All the bins illustrate an asymmetric δα T distribution, with the exception of the region where cosθ µ ≈ 1, with the latter implying that this part of phase-space includes events with minimal FSI effects. Figure 26 shows the comparisons to a number of available neutrino event generators with GiBUU giving the best performance. Figure 27 shows the same results compared to a number of GENIE configurations, illustrating that Gv2 is disfavored in the region where cosθ µ < 0.75. Figures 28 and 29 show the double-differential cross sections as a function of δα T in cosθ p bins. The results in the region with 0 < cosθ p < 0.75 show a fairly flat distribution. The cross section distributions corresponding to forward and backward proton scattering exhibit an FSI-driven asymmetric behavior. Figure 28 shows the comparisons to a number of available neutrino event gen- erators, where NuWro yields a prediction that is disfavored for forward scattering. Figure 29 shows the same results compared to a number of GENIE configurations, illustrating that Gv2 is disfavored across all cosθ p bins. In the -1 < cosθ p < 0 region shown in Fig. 29a, all the predictions illustrate a peak close to 180 • with the exception of Gv2. The driving force for this difference is the Gv2 QE contribution, as can be seen in Fig. 30. This is indicative of potential modeling issues in the Llewellyn Smith QE cross section and of the hA FSI performance used in the Gv2 prediction. Unlike Gv2, the theory-driven GENIE v3 family of predictions (G18, Untuned, and G21) closely follow the data. Figures 31 and 32 show the single-differential cross sections as a function of δφ T using all the events (panel a), as well as the double-differential results as a function of the same kinematic variable in δp T bins (panels b-d). Figure 31 shows the comparisons to a number of available neutrino event generators, with all the genera-tors illustrating a fairly good performance. This result is consistent with the one reported by the T2K collaboration [15,73]. In the lowest δp T region shown in panel b, NuWro is the generator with the best performance. Figure 32 shows the same results compared to a number of GENIE configurations, where Gv2 is disfavored in all regions. At small δp T values the cross section is decreasing and zero above ≈ 40 • which indicates the absence of multi-nucleon and FSI effects. Higher δp T values lead to δφ T cross sections that extend up to 180 • . This behavior is primarily driven by multi-body effects with hadrons below the detection threshold that introduce large kinematic imbalances, as can be seen in panels c-d of Fig. 33. |δp T,y | < 0.15 GeV/c is dominated by QE interactions, while the broader distributions with |δp T,y | > 0.15 GeV/c are mainly driven by MEC events, as can be seen in Fig. 36. In the MEC dominated region of δp T,y < -0.15 GeV/c, all the generators, apart from GiBUU, seem to be lacking in terms of the peak strength. GiBUU seems to be overestimating that MEC contribution in the δp T,y < -0.15 GeV/c bin. With the exception of NEUT, all the event generators illustrate a good performance in the |δp T,y | < 0.15 GeV/c region. Figure 35 shows the same results compared to a number of GENIE configurations, where Gv2 shows the worst performance.
The aforementioned results in kinematic imbalance variables illustrate significant differences across the event generators and configurations used for comparison, especially in the case of the double-differential studies. Yet, the quantity that enters the oscillation probability is the true neutrino energy. Neutrino energy estimators, such as the calorimetric energy E Cal defined in Eq. 6, are used as a proxy for the true quantity. The studies reported next present the results as a function of E Cal in bins of the kinematic imbalance variables. Figures 37 and 38 show the single-differential cross sections as a function of E Cal using all the events (panel a), as well as the double-differential results in the same kinematic variable in δp T bins (panels b-d). Figure 37 shows the comparisons to a number of available neutrino event generators, where the E Cal distribution covers the same energy spectrum across all bins. All the event generators illustrate an equally good performance in the lowest δp T bin. NEUT and NuWro show a deficit relative to the data in the highest δp T bins. Figure 38 shows the same results compared to a number of GENIE configurations, where G18 illustrates the best performance. In the lowest δp T bin, the different configurations illustrate a shift to the left compared to the data, unlike G18, which drives the significantly higher χ 2 values. Interestingly, all the alternative GENIE configurations illustrate a plateau in the highest δp T bin that also yields high χ 2 /bins ratios. Figures 39 and 40 show the double-differential results as a function of E Cal in δα T bins. Figure 39 shows the comparisons to a number of available neutrino event generators. Once again, the E Cal distribution covers the same energy spectrum across all of our results and all the event generators show fairly good behavior. NuWro illustrates a mild deficit in the 135 • < δα T < 180 • bin, which is also reflected in the χ 2 /bins ratio. Figure 40 shows the same results compared to a number of GENIE configurations, where all the GENIE configurations except for G18 illustrate shape and strength differences. Figures 41 and 42 show the double-differential results as a function of E Cal in δp T,y bins. Figure 41 shows the comparisons to a number of available neutrino event generators. All event generators predict very similar cross sections for -0.15 < δp T,y < 0.15 GeV/c (panel a). Unlike this central region, the |δp T,y | > 0.15 GeV/c results yield a wide spread across the generator predictions (panels b-c). Furthermore, apart from GiBUU, all the predictions lack strength in the δp T,y < -0.15 GeV/c bin (panel b). Additionally, NEUT illustrates the same deficit in the δp T,y > 0.15 GeV/c bin (panel c). Figure 42 shows the same results compared to a number of GENIE configurations, where all the GENIE configurations but G18 illustrate a poor performance due to shape and strength issues.

VIII. CONCLUSIONS
This work reports on measurements of flux-integrated differential cross sections for event topologies with a single muon and a single proton detected in the final state using the Booster Neutrino Beam at Fermi National Accelerator Laboratory and the MicroBooNE detector. The data were studied for the first time in the form of single-differential cross sections in kinematic imbalance variables on argon. Furthermore, the first double-differential cross sections in these variables were reported on the same nucleus. Additionally, novel double-differential cross section measurements of a neutrino energy estimator in bins of these variables were presented. The results were compared to a number of event generators and model configurations. The predictions as a function of the energy estimator across all generators and model configurations remain mostly unchanged regardless of the kinematic variable used for the double-differential measurements. Based on the reported χ 2 / bins, the good agreement observed across the calorimetric energy distributions suggests that the energy dependence is largely well-modeled across most predictions. Unlike the energy estimator results, we found that the measured kinematic imbalance cross sections in different phase-space regions are sensitive to nuclear effects. The performance of the event generators and configurations varies depending on the observable of interest. Overall, the GENIE v3.0.6 G18 10a 02 11a cross section predictions with the MicroBooNE-specific tuning (G18) fit the data well. On the other hand, the GENIE v2.12.10 (Gv2) cross section predictions are systematically a poor fit to data with significant shape differences across all variables of interest. The GENIE v3.0.6 G18 10a 02 11a configuration without additional tuning (Untuned) shows a systematic deficit of ∼ 20% which necessitated the development of the aforementioned tune. The GENIE v3.2.0 G21 11b 00 000 configuration (G21) serves as an example of a theorydriven GENIE configuration that shows good agreement with data in most variables without the need for additional tuning. GiBUU 2021 (GiBUU) shows good agreement with data in most kinematic variables, with the exception of δp T , where a systematic shift to higher values of δp T has been identified. A potential source of this shift is due to the GiBUU MEC modeling. The NuWro v19.02.2 (NuWro) prediction falls bellow the data due to poor FSI modeling and shows significant shape differences in FSI-dominated parts of the phase-space. NEUT v5.4.0 (NEUT) also results in predictions mostly falling below the data points. This mismodeling remains largely unnoticed when combined into the calorimetric energy estimator. Yet, future neutrino oscillation measurements will rely on accurate cross section predictions and a precise mapping between measured and true neutrino energies. Therefore, such mismodeling effects might impact their experimental sensitivity. The reported results both provide precision data to benchmark neutrinonucleus interaction models and establish phase-space regions where precise reaction modeling is still needed.   Overflow (underflow) values are included in the last (first) bin. The additional smearing matrix A C should first be applied to the event distribution of an independent theoretical prediction when a comparison is performed to the data reported herein, and then divided by the bin width. The A C matrices are dimensionless. The double-differential cross sections include correlations between the phase-space slices. The data release with the data results, the covariance matrices, and the additional smearing matrices is included in the DataRelease.root file. Instructions on how to use the data release and the description of the binning scheme are included in the README file. Candidate muon-proton pairs are isolated by requiring the existence of precisely two track-like and no shower-like objects, as classified by Pandora using a track-score variable [1,2]. The loglikelihood ratio (LLR) particle identification (PID) score [3] is used to identify the muon and proton candidates. Figure 1 shows the particle composition breakdown of the sample as a function of the LLR PID score. Protons tend to have lower LLR PID score values than muons, thus the track with the lowest score is tagged as the candidate proton, as can be seen in Fig. 1a. Meanwhile, the track with the higher score is treated as the candidate muon, as shown in Fig. 1b 1. (Left) the proton candidate LLR PID score distribution, illustrating the particle composition of the variable. (Right) the muon candidate LLR PID score distribution, illustrating a peak close to one. Only statistical uncertainties are shown on the data. The bottom panel shows the ratio of data to prediction.

UNCERTAINTY COMPONENTS
Below we present the contribution of each one of the sources of systematic uncertainties for each of the variables of interest. "Stat" refers to the statistical uncertainty due to the data, cosmic, MC non-CC1p0π background, and out-of-cryostat samples. "LY" refers to the systematic uncertainty due to light yield detector variations [4]. "TPC" refers to the systematic uncertainty due to wire modification TPC detector variations [4]. "SCERecomb2" refers to the systematic uncertainty due to the space charge (SCE) and recombination detector variations [4]. "XSec" refers to the systematic uncertainty due to the cross-section related uncertainties [5][6][7]. "G4" refers to the systematic uncertainty due to the GEANT4 reinteraction related uncertainties [8]. "Flux" refers to the systematic uncertainty due to the flux related uncertainties [9]. "Dirt" refers to the systematic uncertainty due to the out-of-cryostat strength related uncertainties. "MC Stat" refers to the systematic uncertainty due to the MC CC1p0π uncertainties obtained with a bootstrapping technique. "POT" refers to the proton-on-target uncertainty of our measurement. "NTarget" refers to the number-of-targets uncertainty of our measurement. "Total" refers to all the previous uncertainties combined. 14

GENIE Cross Section Uncertainties
The systematic uncertainty due to the cross section modeling with the GENIE framework is calculated using the relevant EventWeight weights detailed in [7]. The breakdown of the relevant uncertainties are shown in the tables below, both for the total as well as the signal (Sig) and background (Bkg) events.