This is a repository copy of Characterization of nuclear effects in muon-neutrino scattering on hydrocarbon with a measurement of final-state kinematics and correlations in charged-current pionless interactions at T2K

eprints@whiterose.ac.uk https://eprints.whiterose.ac.uk/ Reuse This article is distributed under the terms of the Creative Commons Attribution (CC BY) licence. This licence allows you to distribute, remix, tweak, and build upon the work, even commercially, as long as you credit the authors for the original work. More information and the full terms of the licence here: https://creativecommons.org/licenses/

results are compared to many neutrino-nucleus interaction models which all fail to describe at least part of the observed phase space. In case of events without a proton above a detection threshold in the final state, a fully consistent implementation of the local Fermi gas model with multinucleon interactions gives the best description of the data. In the case of at least one proton in the final state the spectral function model agrees well with the data, most notably when measuring the kinematic imbalance between the muon and the proton in the plane transverse to the incoming neutrino. Within the models considered, only the existence of multinucleon interactions are able to describe the extracted cross-section within regions of high transverse kinematic imbalance. The effect of final-state interactions is also discussed.

I. INTRODUCTION
Neutrino interactions with nuclei are the experimental tool exploited to provide evidence of neutrino oscillations [1][2][3][4][5][6][7][8] and to search for leptonic CP-symmetry violation [9][10][11][12]. In long-baseline accelerator-based neutrino oscillation experiments, neutrino beams are produced with energies in the range of hundreds of MeV to a few GeV. The produced neutrinos interact then with the bound nucleons of nuclei in the detectors via reactions such as quasi-elastic scattering (QE), resonant production (RES), and deep inelastic scattering (DIS). A precise measurement of the oscillation parameters relies on the understanding of the incoming neutrino beam flux, of the scattering of neutrinos with nucleons, and of the nuclear medium effects in the nucleus. The systematic uncertainties arising from neutrino-nucleus interactions, especially those related to nuclear effects, are currently one of the limiting factors for oscillation measurements [13] in T2K [14] and NOvA [15], and will become the dominant uncertainties for future long-baseline experiments, such as DUNE [12] and Hyper-Kamiokande [16].
Neutrinos of such energies can probe nuclear structure at the nucleon level and therefore an accurate description of the nucleus in terms of nucleonic degrees of freedom is essential. To a first approximation, in the independent particle model (IPM), each nucleon is subject to Fermi motion (FM) and a mean-field potential. It is then common to factorise neutrino-nucleus interactions into an interaction with such a bound nucleon (the impulse approximation), leaving the remaining nucleus in a one-particle-one-hole (1p1h) excitation state, and a separate description of the subsequent final state reinteractions inside the nucleus [17]. Driven by precision measurements of electron-nucleus scattering and first large statistics neutrino-nucleus scattering measurements [18,19], various theoretical developments beyond these approximations have been proposed. In the random phase approximation (RPA) approach [20][21][22][23][24], collective excitations approximated as a superposition of 1p1h excitations are calculated. This particular medium effect is parametrised as a correction factor to the interaction cross section as a function of the squared four-momentum transfer Q 2 . In addition to such long-range correlations, short-range correlations (SRCs) are also captured by the spectral function (SF) approach [25][26][27][28], which accounts for nucleon-nucleon correlations beyond the meanfield dynamics. These correlations produce an enhancement in the ground-state nucleon momentum distribution beyond the Fermi momentum, and can lead to twoparticle-two-hole (2p2h) excitations of the nucleus (and, more in general, to npnh excitations with n > 1). Formalisms developed for electron-nucleus scattering have been adapted to describe neutrino data, proposing that 2p2h contributions, notably due to meson-exchange currents (MEC), might be significant in neutrino-nucleus interactions [24,[29][30][31][32][33][34].
Among the reactions relevant for GeV energy neutrinos, the charged-current (CC) QE, is of primary importance for neutrino detection in oscillation experiments, where ν and are the neutrino and the corresponding charged lepton, N and N are the initial-and final-state nucleons. Embedded in a nucleus, the final-state nucleon propagates through and interacts with the nucleus remnant. These final-state interactions (FSI) could be highly inelastic, causing energy dissipation which can prevent hadrons escaping the nuclear medium or alternatively stimulate additional hadrons to be emitted. As a result, the QE reaction in Eq. 1 is not directly accessible. What can be measured are the CC interactions without pion in the extra-nucleus final state (CC0π). This process includes not only other reactions such as pion production, in which the pion is absorbed inside the nucleus, but also 2p2h excitations involving two-nucleon knockout. CC0π (sometimes called "CCQE-like") interactions have been extensively measured [19,[35][36][37][38][39][40][41][42][43][44][45], yet the unambiguous identification of various nuclear effects has proved difficult. This is primarily because the often measured single-particle finalstate kinematics, such as momentum and angular distributions, are determined by both the intrinsic dynamics of Eq. 1 and by nuclear effects. This paper reports measurements of muon-neutrino CC0π interactions with the T2K beam, which has a peak energy of around 600 MeV. The multi-differential cross section using muon and proton kinematics, their correlations, and the final-state multiplicity of protons (above a threshold energy) are measured. These measurements are performed using the T2K near detector (ND280), on a plastic scintillator (C 8 H 8 ) target, with approximately 6 × 10 20 protons on target (POT). The main aim of such new measurements is to improve the understanding of nuclear effects in neutrino interactions, notably with a view to minimizing the corresponding uncertainties in neutrino oscillation measurements. In oscillation measurements neutrino-interaction models are used to infer the neutrino energy from the final state particles and to extrapolate the near detector constraints to the far detector. To test the correctness of such inference, detailed comparisons of the measured cross sections with the most recent neutrino-nucleus interaction models are reported in this paper. The modelling of neutrino energy reconstruction in the CC0π sample, exploited for neutrino oscillation measurements, is affected by large uncertainties due to nuclear effects: even when protons can be in principle detected, the detector response depends on the actual kinematics of the outgoing protons. In the absence of a robust model prediction on the hadronic final state, a multi-differential measurement of single-particle kinematics and nucleon multiplicity, provides valuable input for the modelling of neutrino energy reconstruction and detector response. Furthermore, measurements of proton kinematics from neutrino-nucleus scattering may be used to infer neutron multiplicity and kinematics in the corresponding antineutrino reaction. While the single-particle kinematics and the multiplicity measurements provide a comprehensive description of the CC0π final state, the measurement of muon-proton correlations in the final state provides a powerful probe of nuclear effects. Considering the dynamics of Eq. 1 in the case of scattering on a free nucleon νn → p in the absence of nuclear effects, the final-state proton kinematics can be uniquely determined by that of the muon. In a CC0π measurement the deviation of proton multiplicity and kinematics from what is expected in the simple process of Eq. 1 originates solely from nuclear effects. Such deviations can be characterised using so-called transverse kinematic imbalances (introduced for the first time in Ref. [46]) and proton inferred kinematics, which are measured in the analyses presented here. This paper is organised as follows. After a short description of the T2K experiment in Sec. II, the measurements presented in this paper and the new variables are introduced in Sec. III. Sec. IV describes the analysis procedure, including the simulations used, the event selection and the method for cross section evaluation. Following this the results are reported for each of three analyses: one using proton and muon kinematics, another using transverse kinematic imbalances and a third using proton inferred kinematics. The interpretation of the results is discussed in Sec. V, followed by conclusions in Sec. VI.

II. THE T2K EXPERIMENT
The Tokai-to-Kamioka (T2K) experiment [14] is an accelerator-based long-baseline experiment which measures neutrino oscillations in a ν µ (ν µ ) beam [9]. The T2K neutrino beam is produced by the Japan Proton Accelerator Research Complex. A 30 GeV proton beam collides with a graphite target producing positive and negative pions and kaons which are focused and chargeselected by three horn magnets. The positive (negative) hadrons decay to produce a flux highly dominated by ν µ (ν µ ) [47].
The Super-Kamiokande far detector is located 295 km away from the production point and sits 2.5 o off the beam axis. T2K is further equipped with two near detectors: ND280 and INGRID. INGRID [48] is designed to monitor the direction of the neutrino beam whilst ND280 is dedicated to the study of the un-oscillated spectrum of neutrinos at 280 m from the production target and is the detector used by the analyses presented here. ND280 is positioned off-axis so that it has the same peak neutrino energy as Super-Kamiokande. Such configuration ensures a narrow energy spectrum of the beam centred around 600 MeV, in correspondence with the oscillation maximum. It also suppresses the intrinsic ν e and the non-QE interactions, which are primarily produced by the high-energy tail of the neutrino flux. ND280 is composed of an upstream π 0 detector (P0D) [49] and a central tracker region, described below, surrounded by an electromagnetic calorimeter (ECal) [50], consisting of interleaved layers of lead and scintillator, which itself is all contained within a magnet, providing a 0.2 T dipole field. The magnet is instrumented with the side range muon detector [51]. A schematic of ND280 is shown in Fig. 1.
The primary component of ND280 used in the analyses presented here is the central tracker region, comprising of three time projection chambers (TPCs) [52] and two fine grained detectors (FGD1 and FGD2) [53]. The FGDs are both instrumented with finely segmented scintillating bars which provide both charged particle tracking as well as a target mass for neutrino interactions and, whilst FGD1 is fully active, FGD2 also contains inactive water layers. In these analyses only FGD1 is used as a hydrocarbon (C 8 H 8 ) target. Events leaving the FGDs can be tracked into the TPCs, which provide high-resolution tracking and thereby allow the curvature of charged particles to be used to make accurate measurements of their momenta (the TPCs provide an inverse momentum resolution of 10% at 1 GeV). This can then be combined with measurements of particle energy loss for charged particle identification (PID). If charged particles stop before leaving the FGD1, their momentum is determined by their length. In this case the PID is performed using both track length and the total energy-deposition. Muons and FIG. 1. An exploded view of the ND280 off-axis near detector labelling each sub-detector. Adapted from [14].
pions can also be identified by searching for delayed signal at the track end due to the Michel electron from the decay of muons (including muons from pion decay).

A. Observables
This paper presents three different analyses which study the kinematics of the outgoing muon and protons in charged-current events without pions in the final state (CC0π). Each of these analyses measure differential cross sections as a function of different observables and with a slightly different selection, optimized to the observables being measured.
The first 'multi-differential' analysis measures the differential cross-section as a function of the momentum and angle of the particles in the final state. This approach minimises the dependency of the result on the input neutrino-nucleus scattering simulations, as will be described later, and provides the most complete information to characterise the final state. Such results can therefore be compared with present and future models of CC0π processes, even if their direct interpretation in terms of different nuclear effects is not straightforward. This multi-dimensional analysis simultaneously measures the cross section of events with and without detected protons in the final state, allowing a complete description of CC0π events and, due to improved constraints on systematic uncertainties, surpasses the accuracy of results previously reported by the T2K collaboration in Ref. [44]. Since this analysis classifies events based on the num-ber of reconstructed protons, it is also able to measure a cross-section as a function of the multiplicity of protons above detection threshold. The other two analyses require the presence of at least one proton and, in the case where multiple protons are reconstructed, only the most energetic one is used to form the measured observables.
The second 'single transverse variables (STV)' analysis measures the cross-section of CC0π events with (at least) one proton in the final state as a function of the STV, which are defined in Ref. [46]. The MINERvA experiment are also measuring transverse kinematic imbalances with a ∼ 3 GeV peak neutrino beam energy [54]. These variables are built specifically to characterise, and minimise the degeneracy between, the nuclear effects most pertinent to long-baseline oscillation experiments. In particular, the STV facilitate the possible identification of: Fermi motion of the initial state nucleon, final state re-interactions of the nucleons in the nucleus and multinucleon interactions (2p2h). As shown in Fig. 2, the STV are defined by projecting the lepton and proton momentum on the plane perpendicular to the neutrino direction. In the absence of any nuclear effects, the proton and muon momenta are equal and opposite in this plane and therefore the measured difference between their projections is a direct probe of nuclear effects in QE events: where p N T is the initial state nucleon transverse momentum and ∆ p T is the modification due to final state effects. δ p T can be fully characterized in terms of the vector magnitude (δp T ) and the two angles (δα T and δφ T ): where p l T and p p T are, respectively, the projections of the momentum of the outgoing lepton and proton on the transverse plane. Different nuclear effects alter the distributions of such STV in different and predictable ways. Measurements of the STV therefore have a unique sensitivity to identify nuclear effects, as will be exploited in Sec.V. This allows cross sections extracted using these observables to act as a powerful tool to tune and distinguish nuclear models. Furthermore, in case of disagreement the STV distributions provide useful hints on the possible causes of the discrepancies.
The third 'inferred kinematics' analysis utilises a similar kinematic imbalance to the STV analysis to probe nuclear effects in CC0π interactions by comparing the measured proton momentum and angle with the proton kinematics which can be inferred from the measured muon kinematics in the simplified QE hypothesis. Such The left side shows an incoming neutrino interacting and producing a lepton ( ) and a proton p, whose momenta are projected onto the plane transverse to the neutrino (ν). The right side then shows the momenta in this transverse plane and how the STV are formed from considering the imbalance within it. Taken from Ref. [55] inferred proton kinematics are estimated as follows: where the z axis corresponds to the neutrino direction, n, p, µ and ν denote the neutron, proton, muon and neutrino and E b is the nuclear binding energy. The value of E b used in the definition of these variables is 25 MeV for carbon, but this may be different from the event-byevent "physical" value of E b . The cross section for events with a muon and (at least) one proton in the final state is then measured as a function of three observables: These observables are built such as to enhance nuclear effects which manifest themselves as deviations from zero imbalance. The STV depend only on transverse components of muon and proton momentum vectors with respect to the neutrino direction, while the variables of Eq.8 depend also on the longitudinal components of both vectors. As such, there is no trivial relation between the two sets of variables such that each gives complimentary information about the nuclear effects involved in neutrino interactions. As can be seen in Eq.7, the definition of the inferred proton kinematics relies on the same QE formula as is used in the estimation of neutrino energy in oscillations measurements at T2K. Therefore the observed deviations from the expected proton inferred kinematic imbalance provide hints of the biases that may be caused from the mismodeling of nuclear effects in neutrino oscillations measurements at T2K. The measurement of the differential cross section as a function of these proton inferred kinematic variables is performed separately in bins of muon kinematics. This can highlight the possible mismodeling of nuclear effects in different regions of the muon kinematic phase-space and is also essential in order to mitigate the model dependence in the efficiency corrections (this will be further discussed in Sec. III B). Once de-convoluted from detector effects, this analysis measures how the true particle kinematics deviate from their inferred values under a QE approximation.

B. Minimisation of input-model dependence
In all three analyses extensive precautions are taken to ensure that the results are minimally dependent on the signal model used in the reference T2K simulation (this model is detailed in Sec.IV). This is particularly important for these analyses since the predictive power of available interaction models for the outgoing proton kinematics, and the relative kinematics between muon and protons, is poor. One crucial way to minimise such model-dependence is to ensure that the analyses' signal definition is only reliant on observables which are experimentally accessible at ND280. As such, the signal is defined as all events with no pions in the final state (CC0π) without correcting for FSI pion absorption. Moreover, for the analyses which integrate over large regions of kinematic phase space or do not estimate the efficiency as a function of all relevant kinematic variables, it is also absolutely necessary to apply phase space restrictions in the signal definition in order to avoid model dependence in the efficiency correction. The phase space re-strictions used in the analyses presented here are shown in Tab. I. Since the efficiency of detecting muons and protons in ND280 is not flat as function of the particles angle and momentum, the efficiency correction should be made as a function of the momentum and angle of both the outgoing particles. The relative angle between the outgoing particles is also important but, due to the magnetic field and the very good spatial resolution of the TPCs, this has only a second-order effect on the efficiency. The multi-differential analysis performs a complete multi-dimensional efficiency correction and therefore only a loose phase space restriction on the proton momentum is applied. The STV analysis may, in principle, be the most affected by this issue since each bin of the STV integrates over all possible muon and proton kinematics. As a consequence, the STV measurements use the most stringent restrictions in the signal phase space, selecting only regions of flat and/or well understood efficiency. Finally, the inferred kinematics analysis performs a measurement binned in muon momentum and angle and thus it requires only restrictions on the proton phase space. It should be noted that the restrictions listed in Tab. I are applied in the signal identification at generator level, therefore the multiplicity of the protons is defined counting only protons above the thresholds in the table. The final measurements do not correct for protons which cannot be detected efficiently and therefore the same restrictions have to be applied to any model in order to compare with the results presented in this paper.
To further alleviate model-dependence, the measured differential cross sections are flux-integrated, normalising all the bins of the measured variables to the same flux: where N CC0π i is the measured number of signal events in the i-th bin, i is the efficiency in that bin, Φ is the overall flux integral, N F V nucleons is the number of nucleons in the fiducial volume and x is the measured variable. The analyses can be further affected by modeldependent assumptions in the process of correcting for detector effects. The multi-dimensional and the STV analyses use a binned likelihood fit, similar to that used for 'Analysis I' in Ref. [44]. The results of this method, when unregularised, are completely independent on the nominal model used to create the reference templates for the signal. The STV analysis also provides results after applying a regularisation method which has been tuned and thoroughly tested in order to minimize the dependence on the signal model. The third analysis exploits the D'Agostini unfolding procedure [56,57], also described for 'Analysis II' in Ref. [44].
To additionally reduce model-dependence, and to minimise systematic uncertainties related to background modelling, each analysis employs dedicated control regions to achieve a data-driven background estimation and subtraction. Since the control regions chosen and the background subtraction method differs slightly between analyses, these will be discussed in the details of the strategy for each of the analyses which will be reported in Sec.IV.
Despite the many aforementioned precautions, it is still possible that residual model-dependence can bias analysis results. To ensure this does not happen, a comprehensive set of studies with mock datasets has been performed. A first set of mock datasets is created by modifying systematic parameters of particular interest within the reference model (for example 2p2h normalisation or M QE A ). The cross section extraction methods must be able to recover the truth when each mock dataset is treated exactly as real data. However, this only tests that the methods can extract the truth from mock data which are systematic variations of the input model, and so is more of a closure test than a true evaluation of possible bias. For a more rigorous test, alternative Monte Carlo event generators, which employ some entirely different signal and background models, are used to produce mock data. Moreover, some of these mock data is specialised to specifically modify the models of the nuclear effects that the analyses wish to characterise, namely modifying 2p2h shape, Fermi motion and FSI models. Using such mock data as an input, it has been verified that, even in the case of extreme deviations from the input signal model, the cross section extraction machinery for each analysis can recover the truth such that it is always well within the uncertainties on the extracted result and also produces a small χ 2 when the full resultant covariances are considered. Some examples of such studies can be found in [58].
Finally, it should be noted that the three analyses exploit the same data and rely on similar selections. The systematic uncertainties are also evaluated in similar ways, for instance relying on the same data in control regions. As a consequence, it is a very good approximation to assume all the uncertainties to be fully correlated between the different analyses thus the results of the analysis should not be used together in a joint fit. A full discussion on the interpretation of the results will be reported in Sec. V.

A. Simulation
The analysis of the neutrino data relies on simulation in order to correct the measured quantities for flux normalization, for detector effects and to estimate the systematic uncertainties.
The T2K flux simulation is based on the modelling of interactions of protons with a graphite target using the FLUKA 2011 package [59,60]. The modelling of hadron re-interactions and decays outside the target is performed using GEANT3 [61] and GCALOR [62] software packages. Multiplicities and differential cross sections of produced pions and kaons are tuned based on the NA61/SHINE data [63][64][65] and on other experiments [66][67][68], allowing the reduction of the overall flux normalisation uncertainty to 8.5%.
The neutrino interaction cross-section with nuclei in the detector and the kinematics of the outgoing particles are simulated by the T2K neutrino event generator NEUT 5.3.2 [69,70]. The final state particles are then propagated through the detector material using GEANT4 [71]. Various additional neutrino event generators are used in the analyses presented in this paper in order to both test the robustness of the results (as discussed in Sec. III B) and to compare the final measurements to different models. To this end NEUT 5.3.2.2, NEUT 5.4.0, GENIE 2.12.4 [72], GENIE 2.8.0, NuWro 11q [73] and GIBUU 2016 [74] are used.
NEUT version 5.3.2 utilises the Llewellyn-Smith formalism [75] to describe the CCQE neutrino-nucleon cross section and the spectral function (SF) from Ref. [76] is used as a nuclear model. The axial mass used for quasielastic processes (M QE A ) is set to 1.21 GeV, based on the Super-Kamiokande measurement of atmospheric neutrinos and the K2K measurement on the accelerator neutrino beam [18], while the resonant pion production process is described by the Rein Sehgal model [77] with the axial mass M RES A set to 1.21 GeV. The simulation of multinucleon interactions, when the neutrino interacts with a correlated pair of nucleons, also called 2p2h interactions, is based on the model from Nieves et. al in Ref. [78].
The deep inelastic scattering (DIS), relevant at neutrino energy above 1 GeV, is modeled using the parton distribution function GRV98 [79] with corrections by Bodek and Yang [80]. The FSI, describing the transport of the hadrons produced in the elementary neutrino interaction through the nucleus, are simulated using a semi-classical intranuclear cascade model. 12 GeV) and relies on a different nuclear model for CCQE events: a relativistic Fermi gas (RFG) with Bodek and Ritchie modifications [81]. A parametrized model of FSI is used (known as GENIE's "hA" model). Both GENIE 2.8.0 and 2.12.4 are used within the analyses, the latter facilitates the optional inclusion of 2p2h interactions using the so called 'empirical' MEC model alongside other improvements to the FSI model.
The NuWro 11q version is also used in these analyses. It simulates the CCQE process with the Llewellyn-Smith model, assuming an axial mass M QE A = 1.0 GeV, and the 2p2h process by the model in Ref. [78], similarly to NEUT. Different nuclear models are considered in the comparison to the data: SF, RFG and LFG. For LFG and RFG the effect of Random Phase Approximations (RPA) corrections, as computed in Ref. [32], is tested. RPA is not applied to SF since the model already partially contains the short-and long-range correlations between the nucleons in the nucleus. Similarly the 2p2h contribution should be different in SF with respect to what has been calculated in Ref. [32] for LFG. However, since a dedicated computation of the 2-body current for the SF is not yet available in simulations, the same 2p2h contribution as in LFG is added on top of the SF in both the NEUT and NuWro simulations. For pion production a single ∆ model by Adler-Rarita-Schwinger is used for the hadronic mass W < 1.6 GeV with M RES A = 0.94 GeV. A smooth transition to deep inelastic processes is made for W between 1.3 and 1.6 GeV. The total cross section for DIS is based on the Bodek and Yang approach, similarly to other generators. Like NEUT the FSI are simulated with a semi-classical cascade model.
The measurements presented in this paper are also compared to GiBUU 2016 where the Giessen-Boltzmann-Uehling-Uhlenbeck implementation of quantum-kinetic transport theory [82] is used. The nucleons are inserted in a coordinate-and momentum-dependent potential using the LFG momentum distribution. The CCQE process is modeled as in Ref. [83] with M QE A = 1.03 GeV. The 2p2h contribution is simulated by considering only the transverse contributions and translating to neutrino scattering the response measured in electron scattering [74]. In these comparisons the default GiBUU 2016 initial state isospin for 2p2h interactions is used (T = 1). The model used for single pion production [84] mostly differs from the other generators for the inclusion of medium effects on the ∆ resonance. The DIS is simulated with PYTHIA v6.4.
The comparison of the measurements presented in this paper to the various mentioned models is performed in the framework of NUISANCE [85].

B. Event selection
The three analyses presented herein share a common basic event selection, which aims to identify muon neutrino interactions with a hydrocarbon target producing one muon, no pions and any number of protons in the final state. Events are pre-selected by identifying a vertex in the most upstream fine-grained detector (FGD1) associated with either the highest momentum negative track in the central TPC or, if there is no negative track, the highest momentum positive track. If there is no such TPC track the event is rejected. This pre-selection is split depending on the charge of this primary track, as shown in Fig. 3.
If the primary track is negative then the track is required to be muon-like using the TPC PID. Extra tracks sharing a common vertex with the primary track must either have good quality measurement in the TPC, or be contained in FGD1 such that their kinematics can be reliably determined and they must be identified as protonlike by the TPC or FGD PID respectively. If there is more than one extra track sharing such a common vertex then it is required that at least one of these tracks enter the TPC but each must be identified as proton-like. Following this selection, any events with other tracks that are not muon-or proton-like are rejected. To reject events with low momentum charged or neutral pions, it is required that no Michel electrons (electrons from the decay of the muon that itself is from the pion decay) are tagged within the FGD and that there is no activity in the tracker ECal consistent with a photon. The selected events are then split into samples based on whether there was zero, one or more than one proton-like track and, if so, whether it left a track in the TPC.
If the primary track is positive (and there are therefore no identified negative TPC tracks) then the selection requires the identification of a single extra FGD track sharing a common vertex position with the primary track. This track must then either stop in FGD1 or in the surrounding ECal and be identified as muon-like by the FGD or ECal PID respectively. In the latter case, time of flight information between FGD1 and the ECal is used to ensure that energy depositions seen in the ECal are related to the same track that traversed the FGD.
Finally a last sample is selected with a single track travelling through the FGD before stopping in the ECal. This sample uses the measured time of flight between the track ends to verify propagation direction, and the ECal hit topologies to verify whether the track is muon-like. This is a small sample but all concentrated at high angle, therefore it is included only in the multi-differential analysis which measures the cross-section with finer binning in muon angle. Fig. 3 summarises the topology and the number of selected events within the six signal samples discussed while the number of selected events in each sample, broken down by true interaction topology, is shown in Fig. 4. Other samples are possible but typically with very poor efficiency, resolution and larger detector systematic uncertainties, for example events with a negative primary track and multiple FGD1 contained protons. Since such alternative samples are found to make up a very small number of selected events, less than 30 events in the available data, they are excluded. As discussed in Sec. III B, it is important not to attempt to correct for low efficiency in regions of kinematic phase-space that the detector is not sensitive to. This is particularly important when measuring a differential cross section in observables that do not well characterise a detector's acceptance such as the single-transverse and proton inferred kinematic observables. To avoid inputmodel bias from integrating over regions of changing efficiency, it is necessary to set appropriate limitations on the kinematic phase-space of the final state particles. In the analyses presented here, both muons and protons are identified and therefore ND280's acceptance is reasonably well characterised by the muon and proton momentum and angle 1 . Ideally, the selected phase-space restrictions should leave a flat efficiency within the four dimension regions of muon-and proton-kinematics which will be integrated over in the final measurement. This ensures that the efficiency corrections are independent of the distribution of kinematics which are not measured.
To determine the phase-space restrictions introduced in Tab. I, the efficiency and selected event distributions were studied in various projections of the underlying fourdimensional kinematics in order to find a suitable balance between efficiency flatness and the number of CC0π+Np events that fall out of the restricted phase space (which are then considered as background). The resultant impact of the phase-space restrictions is shown for both the ND280 NEUT 5.3.2 and GENIE 2.8.0 simulations in Fig. 5, which shows the efficiency after all the selection steps projected into the relevant kinematic variables, before and after phase-space restrictions are applied. In general it can be seen that the chosen phase space restrictions ensure a more flat efficiency within the regions of kinematic phase space that contribute most to the CC0π cross section, particularly in the poorly understood outgoing proton kinematics.
Following the event selection and the application of the phase-space restrictions in both true and reconstructed kinematics, the efficiency and purity of signal events for each analysis is shown in Tab. II. The reconstructed muon  Although the selection presented in this chapter identifies a high purity sample of CC0π+Np events, there are still non-negligible backgrounds. The majority of these come from CC1π + events, where the pion (and associated Michel electron) are missed, but there is also notable contribution from other (multi-pion) CC events. These backgrounds are constrained through dedicated control samples which allow an improved background estimation and thereby smaller background modelling uncertainties.
In the multi-differential and STV analyses two control samples are employed for the background constraint. Both require the identification of a negatively charged muon-like track and a positively charged pion-like track in the TPC and are split depending on whether there are any extra tracks sharing a common vertex with the identified muon and pion candidates. The vertex must be contained in the FGD1. These control regions will be referred to as CC1π + and CCOther respectively. An illustration of the topologies these aim to identify is shown in Fig. 8 while the distribution of the data and simulated events within each control sample are shown in Figs. 9 and 10.
These figures highlight an initial large discrepancy between the NEUT prediction and the data, particularly in the CC1π + sample. This is understood to primarily come from an over-estimation of the contribution from neutrino induced coherent pion production, as demonstrated in Ref. [86]. However, the likelihood fit used in the multi-differential and the STV analyses allows to adapt the NEUT model to the data within the control regions. The postfit NEUT prediction from the likelihood fit performed in the δp T measurement is shown in the figures to be in much better agreement with the data (similar results are also obtained in the other STV and multi-differential analyses).
In the inferred kinematics analysis a control sample is built by inverting the cut on Michel electrons in the signal samples. The resultant kinematic distributions of the selected data and MC events are shown in Fig. 11. This control sample is then unfolded simultaneously with the signal regions to constrain the background.

C. Sources of systematic uncertainties
The measurements presented in this paper account for the following systematic uncertainties: • neutrino flux uncertainty. The flux simulation is tuned using external hadron-production measurements and INGRID monitoring, as discussed in Sec. IV A. The residual flux uncertainties affect the cross section measurements presented in this paper, mainly through an overall normalization uncertainty of approximately 8.5%.
• detector effects (efficiency and resolution) which are not perfectly reproduced in the simulation. To evaluate such uncertainties the simulation is compared to the data in dedicated and independent control samples, any observed bias is corrected and the statistical uncertainties in such data and simulated samples are used as residual uncertainties.
• modelling of the signal and background interactions, including nuclear effects. As previously dis-cussed in Sec. III B, a possible model-dependent bias may be introduced in the multi-differential and STV analyses through efficiency corrections, while the measurement of proton inferred kinematics is also affected through the unfolding procedure and the simulation-based background corrections. Such effects are covered by dedicated systematic uncertainties which are quantified by evaluating the variation of the measured cross section using modified simulation models; the theory parameters describing the signal and the background, including proton and pion FSI, are varied inside their prior uncertainty, based on theory expectations and comparisons to external data.
Such uncertainties are implemented in the cross section extraction in different ways in each of the analyses presented in this paper, as will be described in the following sections. In general the systematic parameters considered and their variation is similar to that used for the near detector fit of T2K oscillation analyses described in Ref. [9]: the most notable differences being the inclusion of proton FSI uncertainties and the usage of Gaussian  . The reconstructed kinematics of the muon and of the highest momentum positive (HMP) hadron for events selected within the CC1π + control region from both data and NEUT 5.3.2. The plots are broken down by interaction topology and the CC0π contribution is further split depending on whether the interaction falls within the multi-differential analysis phase space constraints from Tab. I. The postfit NEUT prediction from the likelihood fit to extract a cross section as a function of δpT is also shown.

D. Method of cross section evaluation
Each of the analyses take different approaches when extracting a cross section from the selected events detailed in Sec. IV B. All of these methods involve an effective background subtraction; an efficiency correction; and the deconvolution of detector effects either by a binnedlikelihood fit for the multi-differential and STV analyses, or an iterative unfolding procedure for the analysis of the inferred kinematics.

Binned likelihood fitting
In order to produce a data spectrum that is deconvoluted from detector smearing, the input simulation is varied via a set of parameters, such that a best fit set can be extracted once the simulation best describes the observed data. The signal is parametrised using "template signal weights" (c i ) which alter the number of selected signal events in bins (i) of some truth-level observable(s) with no prior constraint. Parameters describing plausible systematic variations of the flux, detector response and background processes can also be fit simultaneously to the signal parameters. The effect of these parameter variations is then propagated through to the number of selected events in reconstructed bins of the same observable (using the expected smearing due to detector resolution and efficiency), such that the updated simulation prediction can be compared to the data. The best fit set of parameters are chosen by minimising the following negative log-likelihood: Where: and The term in Eq. 11 is the Poisson likelihood, where N sim j and N obs j are the number of simulated and observed events in each reconstructed bin, j. The term in Eq. 12 characterises the prior-knowledge of the values of the systematic parameters ( a syst ) and their correlations, as a multi-variate Gaussian likelihood where a syst prior are the prior values of these parameters and V syst prior is a covariance matrix describing the correlations between them.
As described above, N sim j is described by alterations to the nominal input simulation based on the template signal weights and the systematic fit parameters: describe the alterations to the input simulation from the aforementioned systematic parameters; and U ij is the smearing matrix describing the probability of finding an event in true bin i in reconstructed bin j. This smearing matrix is also subject to change with the alteration of systematic parameters.
The result of the fit is the N CC0π i term from Eq. 9: the number of selected signal events deconvoluted from detector smearing in each analysis bin. As shown in Eq. 9, this must then account for the integrated T2K flux, the number of target nucleons and the bin width before being efficiency corrected to produce a differential cross section.
Such a method of deconvolution is entirely unregularised and is therefore equivalent to using D'Agostini iterative unfolding [56] with an infinite number of iterations or to simply inverting the detector response matrix providing this gives an entirely positive unsmeared spectrum. Provided that the analysis bins do not integrate over regions of phase space of rapidly changing efficiency, this method of unsmearing is completely unbiased but is susceptible to the so-called "ill-posed problem" of deconvolution -where relatively small statistical fluctuations in the reconstructed bins can cause large variations in the fitted contents of true kinematic bins [87]. These results are fully correct and perfectly suitable for further use, for example in fits to constrain parameters in model predictions or to compare the suitability of different models, but they cannot easily be interpreted "by-eye", since they often contain large anticorrelation between adjacent bins which causes the result to strongly 'oscillate' between such bins. Moreover, within the pertinent observables in these analyses, neutrino-interaction cross sections are not expected to follow such an 'oscillating' behaviour. These large variations between neighbouring bins can be suppressed by regularising the results, i.e. imposing smoothness of the fitted parameters c i , thus inducing a small overall reduction of the uncertainties and some dependence of the results on the input signal simulation model. As such, the STV analysis provides both regularised and unregularised results. To achieve this, a regularisation term is optionally added to the likelihood in Eq. 10: Here c i is the signal weight for the i th true bin and p reg controls the regularisation strength. It is clear that this implementation of regularisation adds a constraint which can bias the fit toward the shape of the signal model in the input simulation. However, the impact of the bias can be mitigated by the careful selection of an appropriate regularisation strength. A simple method of choosing p reg in such a regularisation scheme is the 'L-curve' technique presented in Ref. [88]. In this approach a compromise is found between the impact of the regularisation (defined by the normalised regularisation penalty: −2 log(L) reg /p reg ) and the goodness of fit (decreased log(L) reg ). One of the significant advantages of this method, over those typically used to choose the regularisation strength (like tuning the number of iterations) in iterative unfolding methods, is that it is "data-driven": the regularisation strength is determined from assessing the properties of real data and is not solely reliant on simulation studies. It is important to emphasise that the application of regularisation produces a result that is easier to interpret without statistical methods but is at least slightly biased. A regularised result is therefore particularly well suited for result-theory comparison plots but the unregularised result is likely more suitable for forming quantitative conclusions. For this reason unregularised results will be provided in both the multi-differential and STV analyses.

Iterative D'Agostini unfolding
Unfolding accounts for smearing between the true spectrum and reconstructed spectrum due to the detector efficiency and resolution. The relation between true and measured spectrum can be written as : where C i is a number of events in true bin i, E j is a number of events in measured bin j, S ji is a smearing matrix, and N t is the number of true bins. The smearing matrix is constructed from MC predictions which gives the information of event migrations. The Iterative unfolding, proposed by D'Agostini [56,57], uses Bayes' theorem to obtain an unsmearing matrix from the smearing matrix as : where P(E j |C i ) is a probability of the true events in bin i measured in bin j written as : where N ji is the number of true events in bin i measured in bin j. P ef f (E j |C i ) is defined as : where N m is number of measured bins. P 0 (C i ) is a prior probability representing the predicted number of events in bin i, written as : Therefore, the unfolded spectrum is : where N m is the number of bins of measured spectrum. After each iteration, P 0 (C i ) is updated with the posterior of the previous iteration. This method is regularised by choosing the number of iterations, inducing a bias toward the input simulation used. Such bias is tested through multiple mock datasets with alternative simulation models. The number of iterations was chosen by requiring the χ 2 values obtained between the unfolded result and the truth of these mock datasets to reach a stable value: 2-iterations for ∆p p , 6iterations for ∆θ p and 4-iterations for |∆ − → p p |. The bias in the results was shown to always be well within the uncertainties.
Overall this produces an efficiency corrected and unfolded distribution of signal events which must then account for the flux normalisation, the number of target nucleons and the bin width to form a differential cross section, as described by Eq. 9.

E. Multi-differential muon and proton kinematics
This analysis measures the multi-differential cross section of CC0π events as a function of the muon and proton kinematics and the proton multiplicity. As previously described, a multi-dimensional efficiency correction is applied, the cross section is evaluated with a binned likelihood fit and the background is constrained by using dedicated control regions. The binning, reported in Tab. III, is chosen to keep the systematic uncertainty smaller than the statistical uncertainty and to cope with the track reconstruction capabilities of the detector. Due to the small available statistics, the events with two or more protons are all collected in a single bin.
The statistical uncertainties are evaluated by fluctuating the total number of observed event in each bin with a Poisson probability and running the fit multiple times. The systematic uncertainties are evaluated by running the analysis on many toy datasets produced by varying the parameters describing the systematics effects detailed in Sec. IV C. The uncertainties are then found by computing the covariance of the resultant cross sections between every pair of analysis bins. The fractional uncertainties are shown in Fig. 12 for some representative bins. The different sources of systematic uncertainties are shown separately and the total systematic uncertainty is evaluated by simultaneously varying all the nuisance parameters corresponding to the different source of uncertainties. The flux uncertainty is the largest, followed by detector effects. The background modelling uncertainty is sizeable, i.e. of the same order of detector effects, only in the regions with high momentum forward going muons where the background is larger. Finally the signal modelling uncertainty is only non-negligible in the region of backward, low momentum muons where the detector reconstruction capabilities are limited and, due to low available statistics, the angular bin is large, averaging over angles with different reconstruction efficiencies. This is also the region where the backgrounds coming from outside the FGD1 fiducial volume are larger. All these effects tend to increase the dependence of the results on the signal modelling in this particular region of phase space. The statistical uncertainty dominates in most of the bins, except in the regions where the width of the bins is driven by the detector performances. For instance the bin of low proton momentum cannot be further subdivided due to the limited resolution for short tracks. Analogously, in the regions where the muon or the proton have an angle almost perpendicular to the neutrino direction, in order to match the reconstruction capabilities of the detector in absence of a TPC track, the measurement is reported in large bins, and thus systematic uncertainties are larger than statistical uncertainty. Fig. 12 shows the uncertainties on the measurement of proton multiplicity. In this case, integrating over all the bins of muon and proton kinematics, the statistical uncertainty is always smaller than the systematic ones. The dominant uncertainty is still due to the flux, followed by the detector effects, which become very important in the events with two or more protons where the outgoing nucleons have very low momentum and are thus difficult to reconstruct. The uncertainty due to background is completely negligible, while the effect of nucleon FSI becomes important at increasing proton multiplicity due to migration effects between different multiplicity bins. The uncertainty due to signal modelling is very small and more or less constant across all the multiplicities.
The extracted CC0π multi-differential cross section results are shown in Figs. 13 and 14, compared to a variety of different model predictions. Value of the comparisons' χ 2 are reported in the figures. Additional model comparisons to assess the suitability of RFG+RPA nuclear models (Figs. [24][25] and the impact of FSI (Figs. [26][27] are also shown in appendix A 1. The last bin of each of the momenta bins is shortened to improve the plot's readability but these bins are normalised by their total width (as specified in Tab. III) in accordance with Eq. 9. These results are discussed in details in Sec. V. The total CC0π cross section extracted is given in Tab. IV alongside a prediction from the default NuWro 11q simulation (which will be used as a standard reference for the comparison of the integrated measured cross section in all the analyses).
Cross section NuWro prediction 4.329 ± 0.502 3.669 TABLE IV. The total CC0π cross section extracted in units of 10 −39 cm 2 within the multi-differential analysis alongside the prediction from NuWro 11q using an SF nuclear model and with the 2p2h model of Nieves et. al. [78].

F. Single Transverse Variables
In this analysis a novel set of observables is used to directly probe nuclear effects in measurements of CC0π+Np interactions. As previously discussed, these observables exploit the kinematic imbalance between the outgoing lepton and highest momentum proton in the plane transverse to the incoming neutrino, which can act as a powerful probe of nuclear effects both in the initial state and in FSI.
To extract the CC0π+Np differential cross section in    the STV, a binned likelihood fit to the number of selected events in reconstructed STV bins is used, as described in Sec. IV D. The uncertainties are evaluated from the postfit covariance matrix, thereby assuming they are Gaussian distributed. The binning for each of the STV are shown in Tab. V but it should be noted that these STV bins are also restricted to the reduced muon and proton kinematic phase space discussed in Sec. III. The binning is chosen such that the statistical error is comparable to the systematic error and that the bin widths are comparable to the detector resolution.  As described in Sec. IV D, the fit must include effects from plausible variation of detector, flux, and neutrino interaction models. This is achieved by fitting the systematic parameters described in Sec. IV C alongside the signal weights in both the signal region and control samples simultaneously (as described in Eq. 13). The systematic uncertainty due to the impact of these parameters is assessed alongside the statistical uncertainty from a post-fit covariance matrix of the parameters, which itself is constructed from the shape of the likelihood surface close to the best-fit point. All parameters are then marginalised in order to project the uncertainty onto the true number of selected signal events in bins of the STV. For a more conservative error estimation, the prefit uncertainty on the detector and flux systematic parameters is considered when evaluating the uncertainty on the subsequent flux normalisation and efficiency correction. Parameters describing the signal which are overly degenerate with the signal weights (such as M QE A ) are not fit. In such cases their uncertainty is taken into account, without any constraints from the data, by computing the effect of their variation on the efficiency corrections. Such effect is found to be relatively small (less than 2% in all but the last bin of δp T and δφ T where it is about 4%). Further details regarding the uncertainty calculation and the handling of systematic uncertainties can be found in Ref. [58].
The extracted cross sections from the regularised fit in each of the STV are compared to the latest predictions from the current state-of-the-art models from the NEUT 5.3.2.2, NEUT 5.4.0, GENIE 2.12.4, NuWro 11q and GiBUU 2016 neutrino interaction simulations in a variety of configurations. The χ 2 of each comparison is reported within the figures. As discussed in Sec. IV D 1, it can be more useful to instead consider the χ 2 formed from the unregularised result but in this case the conservative data-driven regularisation strength means that the difference in the regularised and unregularised χ 2 is marginal. This is demonstrated and discussed in appendix B (which provides the χ 2 from the comparisons with the unregularised result).
The left plots of Fig. 15 compares the results to a variety of different initial state models whilst a shape-only comparison to a subset of these shown on the right alongside the GiBUU 2016 prediction. The contribution from each interaction mode predicted from NEUT 5.3.2.2 is shown in Fig. 16, alongside the impact of altering the simulations with and without a 2p2h contribution and modifying the FSI strength by varying the mean free path of nucleons within the NEUT FSI cascade model. Finally a similar breakdown by interaction mode is then made for GiBUU and GENIE in Fig. 17. In GENIE the empirical 2p2h contribution, used for instance in the neutrino-interaction model of the NOνA experiment [89], is enabled. Figures to evaluate the impact of RPA and the role of regularisation of the cross section extraction are shown in appendix A 2. These results are discussed in details in Sec. V.
The total CC0π+Np cross section extracted (within the phase space constraints listed in Tab. I) for each of the STV is given in Tab. VI alongside a prediction from NuWro 11q. These total cross sections are not identical since the best-fit parameters are altered slightly depending on which projection of the event selection is used as an input.

G. Proton inferred kinematics
As outlined in Sec. III, this analysis uses the inferred kinematic imbalance between measured proton kinematics and what would be inferred from the measured muon kinematics under a QE approximation, which can act as a metric for the extent to which the QE approximation is reliable for events that are approximately characteristic of the dominant sample used in T2K neutrino oscillation analyses.  In this analysis the muon phase-space is divided into multiple bins in order to correct for the different selection efficiencies and CC-non-QE contributions across the phase-space: • Bin 0 : cos θ µ < −0.6.
Within each muon angular bin, the same binning in the inferred kinematic variables is used, which is shown in Tab. VII. The binning is chosen to ensure that the efficiency is suitably flat within each bin and that the bin width is not less than the detector resolution. To estimate the systematic uncertainties the entire unfolding procedure is repeated for a comprehensive set of plausible variations of the T2K reference model according to the systematic uncertainty sources discussed in Sec. IV C. The covariance of the ensemble of results from the different pseudo-experiments is then taken to characterise the uncertainty: The total CC0π+Np cross section extracted (within the phase space constraints listed in Tab. I) for each of the inferred kinematic observables is given in Tab. VIII alongside a prediction from NuWro 11q. The total cross sections and uncertainties are not identical since the unfolding method couples the systematic parameters between the signal and control regions.

V. DISCUSSION OF THE RESULTS
In all of the results shown in Sec. IV χ 2 statistics are quoted to indicate the agreement between each model and the result, calculated using the full covariance matrix from the cross-section extraction. However, since no model describes the results well over the entire phase space, these χ 2 statistics should be treated carefully: such quantitative estimation of the global data-model agreement is reliable only when a model is capable of describing the whole phase-space of the measurement and ideally when it is fit to the result using a well-motivated parametrisation of the signal predictions. Moreover, since multiplicative normalisation uncertainties make a significant contribution to the covariance of the results, such χ 2 statistics can also suffer from Peelle's Pertinent Puzzle [90] and may therefore not accurately characterise the agreement (this does not effect shape-only χ 2 ). For this reason these statistics should only be taken as an approximate metric and this section will mainly focus on discussing and interpreting discrepancies between the model and the simulations in specific regions of kinematic phase-space, rather than on the overall concurrence indicated by the χ 2 . However, when doing this it is important to be aware of the significant correlations between each bin in the extracted results. The significant detector smearing in the STV leads to adjacent bins being fairly anti-correlated (up to ∼35%, with δp T the most affected) whilst the large flux normalisation uncertainty correlates all other bins (by ∼10% to 35%). Due to the larger regularisation strength in the inferred kinematics analysis and the coarser binning in the multi-differential analysis,   FIG. 22. The extracted differential cross section as a function of the inferred and true proton three-momentum difference in different muon kinematic bins, within a restricted proton kinematic phase space, compared to the NEUT 5.3.2.2 simulation with various scalings of the mean free path of nucleons undergoing FSI processes to simulate different FSI strengths. A comparison of the NEUT prediction without a 2p2h contribution is also shown. 2p2hN indicates the 2p2h model is an implementation of the Nieves et. al. model of Ref. [78]. More details of these models can be found in Sec. IV A. Note that the last bin in each plot is shortened for improved readability. The inlays show the same comparisons on a logarithmic scale.  [78]. More details of these models can be found in Sec. IV A. Note that the first and last bin in each plot is shortened for improved readability. The inlays show the same comparisons on a logarithmic scale. the anti-correlations are generally much less prominent and most bins are positively correlated by the flux. The full covariance matrices are available in the data release for these results. A summary of the full and shape only χ 2 for the regularised and unregularised STV results are provided and discussed in appendix B.
The measurement of proton multiplicity and proton and muon kinematics from the multidifferential analysis (Figs. 13, 14) shows the phase space regions where the present models fail to describe the data. From this measurement it can be seen that, when there is no proton above threshold in the final state, while the SF prediction gives a reasonable agreement with the extracted result in the region (0.6 < cos θ < 0.8) , it clearly underestimates the cross section in the region of backward muon angle and overestimates it in the region with forward muons for intermediate muon momentum (0.5-0.8 GeV). Conversely, the SF model describes very well the rate of events with one proton above the momentum threshold in the final state, where a slight preference for the presence of 2p2h is observed in the region with forward muons (cos θ > 0.8).
The LFG predictions from NEUT 5.4.0 and NuWro 11q differ when describing events without protons above threshold in the final state in the region with muons at high angle. Here the NEUT implementation describes the results well, while NuWro underestimates the cross section, in a similar manner to its SF prediction. In the phase space with intermediate and low muon angle, both LFG implementations describe the result reasonably well. However, both also overestimate the cross section with one above-threshold proton in the final state. Indeed, the two LFG models predict a tendency in proton multiplicity to have a larger rate of events with one proton with respect to events without protons, while the extracted result has the opposite behaviour.
Finally, in the bin with two or more protons in the final state, the result prefers the SF with 2p2h over the case without 2p2h. Here it can also be seen that the two implementations of LFG give very different predictions, thereby demonstrating the importance of 1p1h modelling also for the events with multiple protons.
It is also interesting to compare the effect of FSI and 2p2h in these distributions, as shown in Figs. 26-27 in appendix A 1, where SF with and without 2p2h is compared together with different strengths of FSI. A larger FSI strength tends to redistribute events between the bins of proton multiplicity: larger FSI increases the rate of events without protons above the momentum threshold and decreases the number of events with one or more protons above it. On the other hand, 2p2h tends to increase the cross section for all proton multiplicities. Since the measurement of the shape of the cross section in proton multiplicity is well known (the uncertainties are dominated by effects that fully correlate all three bins, mostly due to the flux) these results may offer an interesting capability to separate the effects of proton FSI and 2p2h. However, more robust predictions of the outgoing proton kinematics in 1p1h and 2p2h events after FSI would be needed in order to exploit the proton multiplicity measurement to expose a possible significant 2p2h excess in the result. It is particularly striking that no model is able to simultaneously describe accurately the events with and without above-threshold protons in the final state, the LFG being more in agreement with the result in former and the SF being in better agreement in the latter.
The results are also compared in Figs. 24-25 with RFG models as implemented in NuWro 11q and NEUT 5.3.2.2 and with the RFG with Bodek-Ritchie corrections as implemented in GENIE 2.12.4. GENIE overestimates the cross section in most of the phase space, except for backward-going muons. It also reproduces the same trend as the extracted result in the multiplicity plot, showing more events without a final state proton above the threshold. The other RFG implementations in NuWro and NEUT behave similarly to the NuWro LFG model.
It also is interesting to consider the impact of the RES pion-production contribution to the signal (where the pion is absorbed inside the nucleus). In general NEUT 5.3.2.2 predicts that the RES contribution to the cross section is always less than about 5% except when there is one above-threshold proton in the final state and cos(θ µ ) > 0.8, where it reaches around 15%, and when there is more than one proton in the final state where the cross section is dominated almost equally be RES and 2p2h interactions. Although there is large theoretical uncertainty on this RES contribution (notably from nuclear effects in pion production and in pion FSI), within most bins unrealistically large changes would be required to alter the interpretation of the results.
In general, the interpretation of the aforementioned discrepancies between the result of the multidifferential analysis and different simulations is not straightforward since the measured cross section is affected by multiple initial state and final state nuclear effects which cannot be easily separated in the momentum and angular kinematic distributions. The STV are expressly designed in order to unambiguously distinguish the impact of different nuclear effects, their measurement therefore offers a more transparent interpretation of such discrepancies.
The comparison of the STV distributions to different CCQE models, as implemented in NuWro, are shown in the left plots of Fig. 15. Given the definition of δ p T (Eq. 2), the dominant contribution below the Fermi surface (δp T ∼ 230 MeV) is CCQE with limited FSI strength (as can be seen in the left plots of Fig. 16), thus the Fermi motion determines the bulk structure of δp T , thereby allowing it to act as a probe of the initial-state nucleon. The measured δp T distribution strongly disagrees with the RFG prediction: the prominent imprint of the cliff at the Fermi surface, a characteristic of RFG, is firmly disfavoured by this result.
It is interesting to note that both Fermi gas models (LFG and RFG) exhibit similar excesses over the result, but at different kinematic regions. Indeed, considering shape-only comparisons of data to various simulations in the right plots of Fig. 15, it can be seen that the LFG predictions well describe the differential distribution, but are plagued by an overestimation of the overall cross section, even if RPA corrections are applied. Such a normalisation discrepancy could come from a general overly large CCQE cross section or weak proton FSI keeping too many protons above signal threshold. The latter is particularly further supported in the proton multiplicity plot in Fig. 14, where an increase in proton FSI would migrate events from the 0 proton bin to the 1 proton bin thereby bringing the prediction into better agreement with the results. The fact that the GiBUU (LFG) CCQE prediction in Fig. 17, which largely differs from the NuWro and NEUT LFG models through its FSI modelling, seems to provide a normalisation in good agreement with the results adds additional evidence that the normalisation discrepancy seen in the NEUT and NuWro LFG normalisations is at least partially related to FSI modelling.
In general it seems that the nucleon dynamics for δp T 400 MeV, are better described by SF than Fermi gas models. The consistency between SF and the result at δp T ∼ 300 MeV suggest that the nucleon-nucleon correlations captured by SF are required. Future measurements of the STV with higher statistics may allow further exploration of the nature of such correlations.
Above ∼ 400 MeV, δp T is driven by nucleon-nucleon correlations and FSI effects, so it is not surprising that the predictions from the Fermi gas and SF models become more similar. The SF model in this region, as it is implemented in the simulations, is not fully consistent since a 2p2h contribution computed for an LFG model is added on top of the CCQE SF. The SF model without a 2p2h contribution is also shown for comparison: within the hard tail of δp T and δφ T the result clearly indicates the need for additional strength, consistent with that from a 2p2h contribution, beyond the nucleon-nucleon correlations already included in the SF model.
Both RFG and LFG models have consistent predictions regarding the total cross section and the δα T distribution, which represents to good approximation the direction of the initial nucleon momentum p N . Furthermore, the distributions of δα T show a significant difference between Fermi gas and SF models in the shape. In fact, in NuWro predictions the discrepancy at low δα T between Fermi gas models and the data is caused by RPA (see the left plots of Fig. 28), without which the shape would be consistent.
The left hand plots of Fig. 16 show the comparison of STV with NEUT 5.3.2.2 model, demonstrating in more detail how the 2p2h contribution is clearly located in the tails of δp T and δφ T , where the agreement with the data is good and the 2p2h contribution seems essential. It also highlights the CCQE dominance in the bulk of the δp T distribution, where the model tends to overestimate the data.
The distribution of δp T beyond 400 MeV and the shape of δα T are also sensitive to intra-nucleus momentum ex-change such as 2p2h, as already discussed, as well as FSI. As can be seen in right plots of Fig. 16, in order to bring the NEUT 5.3.2.2 model in to agreement with the data, the proton FSI strength must be increased by reducing the mean free path between re-interactions inside the nucleus by a a factor of two. Although it is challenging to draw firm conclusions from external measurements of electron-nucleus scattering, this appears to be around the maximum plausible variation based on the data [91]. Moreover, in the present semi-classical model of FSI implemented in all the simulations, the FSI mainly affect the probability of observing the outgoing proton (here defined as proton momentum above 450 MeV, see Tab. I) thus changing the integrated cross section. Only small modifications to the shape of the STV distributions are visible. As can be seen in Fig. 16, as the FSI strength increases, both δp T and δφ T spectra become harderwith depletion and enhancement in regions of small and large imbalances, respectively-as is expected from the intra-nucleus momentum transfer during FSI. Nevertheless, the enhancement in this particular FSI model is much smaller than that caused by the presence of 2p2h in the region of high transverse kinematic imbalance, so it is far too small to invalidate the evidence for a 2p2h contribution.
The GENIE predictions in the left plots of Fig. 17 strongly overestimate the data in the collinear regions where the proton momentum aligns with the threemomentum transfer in the transverse plane (see p p T and q T in Fig. 2), i.e. in regions where δp T → 0, δφ T → 0, and δα T → 0 and 180 degrees. Such overprediction originates from the elastic interaction of GENIE's widely used hA FSI model [46]. Moreover, compared to other 2p2h models, the empirical MEC model in GENIE features a much stronger enhancement in regions of large imbalances, where the overall predictions clearly overestimate the data.
The GiBUU predictions (the right plots of Fig. 17) provide an integrated cross section in good agreement with the data but it is characterised by one of the hardest δp T and δφ T distributions of all the simulations, as can be seen in right plot Fig. 15, which is in disagreement with the measured shape of the observables. However, it should be noted that there is theoretical motivation to reduce GiBUU's 2p2h model strength by a factor of two [74] and an exploration of the results presented here within a more recent version of GiBUU with such reduced 2p2h has recently shown much better agreement [92]. Fig. 17 and the left hand plots of Fig. 16 also demonstrate that the contribution from RES pion-production within the STV restricted phase-space is universally small relative to the 2p2h contributions (except perhaps in the final bins of δp T and δφ T ). For this reason even relatively large changes in the prediction for the RES contribution to the result would not invalidate the conclusions regarding 2p2h.
Finally Figs.18-20 show the results of the cross section measurement as a function of the proton inferred kinematics (∆p p , ∆θ p , |∆p p |, as in Eq. 7), compared to different models. The most precise measurements come from the region with largest statistics: 0.0 < cos θ < 0.8, p µ > 250 MeV, which corresponds to the region of intermediate Q 2 . Similarly to what was observed in the STV analysis, this region is best described by the SF model and in the high |∆ − → p p | tail a net preference for the presence of 2p2h contribution is visible. Such indication is independent of the strength of FSI effects, as shown in Figs. 21-23, where there is also a small preference for larger FSI. The other regions of muon kinematics are not all consistently well described by any of the models, as can be seen by the high χ 2 values. Depending on the kinematic region and on the observable considered, the LFG or SF may better describe the data. It is also interesting to note that both the very forward-going and low momentum muon kinematic bins suggest a multitude of regions of ∆p p and |∆p p | which are largely dominated by the 2p2h contribution and largely independent of FSI variation, as can be seen in figures Figs. [18][19] respectively. However whether the result favours these large 2p2h contributions depends on the kinematic bin: for example both the NEUT and NuWro SF and 2p2h predictions largely agree with the result in the 0.8 < cos θ < 1.0, p µ > 750 MeV bin, but are quite different in the 0.8 < cos θ < 1.0, 250 < p µ < 750 MeV bin. This difference is understood to stem from the substantial RES pion-production contribution to these bins, which differs considerably in NuWro and NEUT. These bins therefore offer a powerful probe of the CCnonQE contribution, but specific conclusions are difficult due to the large uncertainties in both 2p2h and nuclear effects in RES interactions, making these results complementary to the multidifferential and STV analyses that mostly have a very small RES contribution.
In appendix A 3 a comparison to RFG models is reported which, consistently with what is observed in the STV analysis, generally gives worse agreement with the result than either LFG or SF. In particular the results clearly show that both the NEUT 5.3.2.2 RFG+RPA model and the GENIE 2.12.4 BRRFG model, which are similar to the nominal models used in recent T2K [9] and NOνA [2] oscillation analyses respectively, do not describe the result well. This conclusion is supported by also considering the poor agreement between the multidifferential results and these models seen in appendix A 1, and even more so considering the aforementioned strong contention between the STV analysis results and the RFG and empirical MEC models.

VI. CONCLUSIONS
This paper has presented the measurement with ND280 data of muon and proton multi-differential cross section, as well as their kinematic correlations, for charged-current neutrino-nucleus scattering without pion in the final state. The muon-proton correlations in the final state are measured through two new sets of variables never used before as observables in neutrino cross section measurements: the STV and the proton inferred kinematics. The analysis selection separates events without protons, with 1 proton or with more than 1 proton in the final state above a 500 MeV momentum threshold, thus enabling also the measurement of proton multiplicity. Particular care has been taken to minimise the model-dependence of the measurement in efficiency corrections, background subtraction and cross section evaluation, thus enabling an unbiased and large set of model comparisons with the results. For the first time in neutrino scattering measurements the concept of data-driven regularisation is introduced to achieve a result that is easy to interpret but contains minimal bias. Overall the results offer a powerful new probe of nuclear effects and that their exploration of kinematic imbalances facilitates a method of separating, at least partially, the different contributions of a CC0π measurement. Since prediction power of proton kinematics in neutrino-interaction simulations is still poor, it remains challenging to draw firm quantitative conclusions. Nevertheless, an extensive comparison with generator predictions has allowed interesting qualitative conclusions to be drawn. As briefly summarised below, the three analyses suggest similar conclusions.
The RFG model is able to describe only a very limited region of phase space (and only when there is no abovethreshold proton in the final state) and is categorically disfavoured when considering the result in δp T . The LFG prediction shows slightly better agreement with data than RFG when considering interactions with abovethreshold protons, especially considering the distribution of STV, but it still overestimates the soft part of the STV spectrum. A more consistent LFG implementation, such as the one in NEUT 5.4.0, gives improved results. It provides better agreement with both the STV distributions and in the region without above-threshold protons in the final state and with large muon angle, where no other model is able to describe the result. For the events with one or more above-threshold protons in the final state, the best description of the data is given by the SF model. Beyond the nucleon-nucleon correlations already included in SF, a clear requirement for a 2p2h contribution is visible in the hard tail of STV distributions (δp T and δφ T ) and of the |∆ − → p p | distribution. The requirement for a large 2p2h contribution to the result remains even with dramatic variations of the semi-classical FSI models available in the simulations. On the other hand, GENIE's "empirical MEC" 2p2h model appears to substantially overemphasize the hard tail of the STV. The prominent features of the STV predicted by GENIE's hA FSI model (driven by its elastic component) are in very poor agreement with the results.
The results with one proton in the final state, when compared to the SF model, suggest the need for stronger FSI effects, at the limit to what is allowed by external data of proton-nucleus scattering. The measurement of proton multiplicity can in principle disentangle 2p2h from FSI effects, with the former increasing the cross section in all bins of proton multiplicity while the latter redistributing events between different bins. Currently the primary limitation is the absence of a model able to properly describe both the events with and without above-threshold protons in the final state.
The measurement of neutrino-nucleus interactions with a pionless final state with protons clearly shows the potential to provide an even more detailed characterisation of nuclear effects in neutrino-nucleus scattering in the future. To this aim, larger statistics are needed, alongside more robust predictions of outgoing proton kinematics in 2p2h and FSI models.
As discussed in Sec. IV F this appendix provides additional comparisons of the single transverse analysis results to assess the impact of RPA on the data-simulation comparisons and also to identify the role of regularisation in the cross section extraction procedure. This is shown in Fig. 28 which first shows NuWro 11q predictions with and without RPA compared to the same regularised results. Alongside this, the figure also shows the predictions of different nuclear models, already shown in Fig. 15, compared to the unregularised results. In this final figure it should be noted that, while the nominal result is different than the regularised result presented in Sec. IV F, the physics conclusions remain the same with a similar goodness of fit.

Inferred proton kinematics measurement
As discussed in Sec. IV G this appendix provides additional comparisons of the proton inferred kinematics analysis results to RFG nuclear models. This is shown in Fig. 29 As discussed in Sec. IV D 1, the regularisation of crosssection extraction methods allows results that are easy to interpret at the cost of some bias. Although the regularised STV results presented in Sec. IV F minimise this using a data-driven method, the unregularised results are also produced to guarantee no unfolding bias. To demonstrate that the application of regularisation does not alter the physical interpretation of the results, Tab. IX to XI show a summary of the χ 2 agreement between the various models considered in Sec. IV F and both the regularised and the unregularised results. The shape-only χ 2 is also shown for each model as this does not suffer from Peelle's Pertinent Puzzle discussed in Sec. V and so may be the most useful quantitative metric for comparing model agreement. This shape-only χ 2 is formed by numerically decomposing the full covariance matrix into a shape and normalisation component. Occasionally it was found that it was necessary to add a small component to the diagonal of the shape-only covariance matrix to make the matrix invertible, thereby slightly overestimating the error (this is responsible for the two instances of seeing a shape-only χ 2 greater than the full χ 2 in Tab. XI).
Overall the χ 2 tables show very little difference between the regularised and unregularised result, demonstrating that the cautious data-driven choice of regularisation strength means that only a very small bias is added. The tables echo the conclusions drawn in Sec. V: a 2p2h contribution and a non-RFG nuclear model is shown to be a hallmark of the lower χ 2 in the comparisons.  [81]. The NEUT and GENIE predictions shown here are similar to those used as a starting point for the T2K and NOνA experiment's oscillation analyses respectively. More details of these models can be found in Sec. IV A. Note that the last bin in each momentum plot is shortened for improved readability. The same models are also compared to the cross section as a function of proton multiplicity. More details of these models can be found in Sec. IV A. Note that the last bin in each momentum plot is shortened for improved readability.    [78]. More details of these models can be found in Sec. IV A. Note that the last bin in each momentum plot is shortened for improved readability.    FIG. 29. The extracted differential cross section as a function of the inferred and true proton momentum difference in different muon kinematic bins, within a restricted proton kinematic phase space, compared to GENIE 2.12.4, NuWro 11q and NEUT 5.3.2.2 predictions which utilise an RFG or RFG+RPA nuclear model. GENIE's RFG model also includes the empirical correction from Bodek and Ritchie (BRRFG) [81]. The NEUT and GENIE predictions shown here are similar to those used as a starting point for the T2K and NOνA experiment's oscillation analyses respectively. More details of these models can be found in Sec. IV A. The 2p2h subscript indicates whether it is an implementation of the Nieves et. al. model of Ref. [78] (N) or the GENIE empirical 2p2h model (E). More details of these models can be found in Sec. IV A. Note that the first and last bin in each plot is shortened for improved readability. The inlays show the same comparisons on a logarithmic scale. FIG. 30. The extracted differential cross section as a function of the modulus of the inferred and true proton three-momentum difference in different muon kinematic bins, within a restricted proton kinematic phase space, compared to GENIE 2.12.4, NuWro 11q and NEUT 5.3.2.2 predictions which utilise an RFG or RFG+RPA nuclear model. GENIE's RFG model also includes the empirical correction from Bodek and Ritchie (BRRFG) [81]. The NEUT and GENIE predictions shown here are similar to those used as a starting point for the T2K and NOνA experiment's oscillation analyses respectively. The 2p2h subscript indicates whether it is an implementation of the Nieves et. al. model of Ref. [78] (N) or the GENIE empirical 2p2h model (E). More details of these models can be found in Sec. IV A. Note that the last bin in each plot is shortened for improved readability. The inlays show the same comparisons on a logarithmic scale. FIG. 31. The extracted differential cross section as a function of the inferred and true proton outgoing-angle difference in different muon kinematic bins, within a restricted proton kinematic phase space, compared to GENIE 2.12.4, NuWro 11q and NEUT 5.3.2.2 predictions which utilise an RFG or RFG+RPA nuclear model. GENIE's RFG model also includes the empirical correction from Bodek and Ritchie (BRRFG) [81]. The NEUT and GENIE predictions shown here are similar to those used as a starting point for the T2K and NOνA experiment's oscillation analyses respectively. The 2p2h subscript indicates whether it is an implementation of the Nieves et. al. model of Ref. [78] (N) or the GENIE empirical 2p2h model (E). More details of these models can be found in Sec. IV A. Note that the first and last bin in each plot is shortened for improved readability. The inlays show the same comparisons on a logarithmic scale.