First combined measurement of the muon neutrino and antineutrino charged-current cross section without pions in the final state at T2K

This paper presents the first combined measurement of the double-differential muon neutrino and antineutrino charged-current cross sections with no pions in the final state on hydrocarbon at the off-axis near detector of the T2K experiment. The data analyzed in this work comprise 5.8$\times$10$^{20}$ and 6.3$\times$10$^{20}$ protons on target in neutrino and antineutrino mode respectively, at a beam energy peak of 0.6 GeV. Using the two measured cross sections, the sum, difference and asymmetry were calculated with the aim of better understanding the nuclear effects involved in such interactions. The extracted measurements have been compared with the prediction from different Monte Carlo generators and theoretical models showing that the difference between the two cross sections have interesting sensitivity to nuclear effects.


I. INTRODUCTION
Current and future long-baseline neutrino oscillation experiments have as primary goals the measurements of the Charge-Parity (CP) violating phase (δ CP ), the neutrino mass ordering and the octant determination of the mixing angle θ 23 [1][2][3][4]. To this end, the associated systematic error must be minimized. At present, the limited knowledge of (anti)neutrino-nucleus interactions dominates the uncertainties [5,6]. The main obstacles behind a better understanding of such interactions are a result of limited modeling of the nuclear dynamics and the difficulties in measuring its effect on the cross section. Despite theoretical and experimental efforts in investigating the (anti)neutrino-nucleus cross section during the last decade, a comprehensive picture has not yet emerged [7,8].
Measured values of the muon neutrino and antineutrino charged-current quasi elastic scattering (CCQE) cross sections at K2K [9], MiniBooNE [10,11], MI-NOS [12] and SciBooNE [13], and more recently by T2K [14][15][16][17] and MINERVA [18][19][20][21][22][23][24][25][26] were found to be higher than predictions obtained using the Relativistic Fermi Gas (RFG) nuclear model. The results favored * also at INFN-Laboratori Nazionali di Legnaro † also at J-PARC, Tokai, Japan ‡ affiliated member at Kavli IPMU (WPI), the University of Tokyo, Japan § also at National Research Nuclear University "MEPhI" and Moscow Institute of Physics and Technology, Moscow, Russia ¶ also at the Graduate University of Science and Technology, Vietnam Academy of Science and Technology deceased * * also at JINR, Dubna, Russia † † also at Nambu Yoichiro Institute of Theoretical and Experimental Physics (NITEP) ‡ ‡ also at BMCC/CUNY, Science Department, New York, New York, U.S.A. a higher value of the nucleon axial mass (M QE A ) than those previously measured in bubble chamber experiments using deuterium as targets and pion electroproduction data [27][28][29]. Furthermore, the CCQE muon neutrino and antineutrino cross sections measured by the NOMAD Collaboration at energies above 3 GeV are in agreement with a value of M QE A around 1 GeV/c 2 [30]. This discrepancy highlighted the need for a more detailed description of the (anti)neutrino-nucleus scattering in the few-GeV energy region.
In a muon neutrino CCQE interaction a negatively charged muon and proton are produced via W exchange with a neutron, while in an antineutrino interaction of the same type a positively charged muon and neutron are produced via W exchange with a proton: ν µ +n → µ − + p ν µ +p → µ + + n. (1) As modern long-baseline neutrino experiments use relatively heavy nuclei as targets, nuclear dynamics plays an important role in the interpretation of the (anti)neutrino oscillations. Several theoretical models have proposed that these effects may explain the observed anomalies between data and Monte Carlo (MC) simulations [31][32][33][34][35][36][37][38][39][40][41][42]. If nuclear effects are considered, the particles produced in the (anti)neutrino-nucleus interaction can interact with other nucleons before exiting the nucleus. These, so-called Final State Interactions (FSI), can alter the type, number and kinematics of particles that exit the nucleus after such interactions. For example, in a muon neutrino resonant pion production interaction a pion, a proton and a muon are produced. The pion could be reabsorbed by the nucleus, with the result that only the muon and the proton are observed to exit the nucleus. This would be indistinguishable from a CCQE interaction. Anyway, the observed discrepancies between data and models cannot be explained by FSI alone. Martini et al. [31] indicates that further contributions to CCQE-like processes arise from two (or more) interacting nucleons, referred to as 2p2h excitations or multi-nucleon knockout. Such interactions eject low-energy nucleons (200-500 MeV), which cannot be easily detected. Multi-nucleon knockout is expected to be less significant for antineutrinos relative to neutrinos; in particular it has a different role in the vector-axial interference term which differs by a sign for the neutrino and antineutrino cross section [43]. Studying differences between CCQE-like cross sections for neutrino and antineutrino interactions could provide information about the role of multi-nucleon knockout in (anti)neutrino-nucleus interactions. The sum and the difference of the neutrino and antineutrino CCQE-like cross sections could yield this information. In the sum, the axial-vector interference term is eliminated whereas the difference isolates this term [29]. In Ref. [44], the predicted sum and difference of the neutrino and antineutrino CCQE-like cross-sections are compared with the equivalent values computed using the CCQE-like double-differential cross sections obtained by the MiniBooNE experiment. The analysis found that additional nuclear effects, other than FSI, would be needed to explain the discrepancy between the observed and predicted values of the sum and difference. However, the analysis in Ref. [44] was limited as the neutrino and antineutrino beams peaked at different energies and the two cross sections were measured independently, implying that correlations between the two data sets were not taken into account. A more rigorous analysis can be performed at the T2K near detector complex. Data has been taken with neutrino and antineutrino beams, both of which peak at the same energy. Combining the two data sets can exploit the correlation between them leading to a more precise crosssection measurement. This paper reports the first combined measurement of the double differential neutrino and antineutrino charged current cross sections on hydrocarbon without pions in the final state. This CCQElike cross section will include contributions from CCQE events as well as events in which a pion was produced and then reabsorbed by the nucleus and multi-nucleon knockout events. The neutrino and antineutrino cross sections were used to compute the sum, difference and asymmetry. The neutrino-antineutrino cross-section asymmetry, which is the ratio between the difference and the sum, is a crucial quantity to control in order to avoid biases in the search for CP violation in neutrino oscillation. All these quantities have been compared with predictions made using several MC generators and models, which are discussed in this paper.
The paper is organized as follows. The T2K experiment is described in Section II. The data and MC samples used in this analysis are reported in Section III. Section IV describes the analysis procedure, including the event selection and the cross-section extraction method. Finally the results and their interpretation are discussed in Section V, followed by conclusions in Section VI.

II. THE T2K EXPERIMENT
The Tokai-to-Kamioka (T2K) experiment [1] is a longbaseline experiment that studies neutrino oscillations in an accelerator-produced ν µ (ν µ ) beam. The neutrino beam, produced by the J-PARC facility, utilizes a 30 GeV proton beam. A proton spill consisting of 8 bunches with 580 ns spacing is produced every 2.48 s. At a beam power of 430 kW, this spill and repetition rate correspond to 2.25 × 10 14 protons on target (p.o.t) per spill. Secondary hadrons, mainly pions and kaons, are produced when the proton beam interacts with a graphite target. Three magnetic horns are used to perform focusing and charge selection of the pions and kaons. The polarity of the magnetic horns can be changed to select positively (forward horn current) or negatively (reverse horn current) charged pions and kaons to produce a beam that is predominantly made of ν µ in the forward horn current case orν µ for the reverse horn current. The selected hadrons decay in a 96 m long decay volume, to produce a (anti)neutrino beam whose direction is parallel to that of the initial proton beam. Both the neutrino and antineutrino beam consist of a mixture of ν µ ,ν µ , ν e andν e . The compositions of the neutrino and antineutrino beams are shown in Fig. 1. In the neutrino beam mode, the "right-sign" ν µ (and ν e ) flux is around 15% higher around the flux peak when compared with the right-signν µ (andν e ) flux in the antineutrino mode. The background antineutrino flux is also lower in the neutrino mode compared with the neutrino flux in the antineutrino mode, especially at high energy. These differences can be attributed to the higher production multiplicities of positively, rather than negatively, charged parent particles.
The Super-Kamiokande far detector is located 2.5 • off the beam axis, at a distance of 295 km from the production point. The near detector complex, located 280 m downstream from the production target, contains two sets of detectors: INGRID and ND280. INGRID [45] is on-axis and monitors the flux and direction of the neutrino beam. The ND280 detector is positioned 2.5 • offaxis and is used to study the unoscillated beam. At an off-axis angle of 2.5 • , the energy spectrum of the beam is narrowed and centered around 600 MeV, which corresponds to the oscillation maximum for a baseline of 295 km. In addition, this narrow energy spectrum suppresses the intrinsic ν e (ν e ) and non-quasi-elastic interactions, leading to lower intrinsic backgrounds to the ν e (ν e ) appearance search at the far detector. This work has been performed using the off-axis near detector, ND280. Fig. 2 shows a schematic of such detector. The ND280 detector is formed from five subdetectors; an upstream π 0 detector (P∅D) [46], two Fine-Grained Detectors (FGDs) [47], three Time Projection Chambers (TPCs) [48], Electromagnetic Calorimeter (ECal) [49] and a Side Muon Range Detector (SMRD) [50]. The P∅D, FGDs, TPCs and ECal are contained within a magnet that provides a 0.2 T field, whilst the SMRD is embedded in the magnet. The measurements reported in this paper used the FGDs, TPCs, ECal and SMRD to select charged-current ν µ andν µ interactions. The most upstream FGD (FGD1) is formed from layers constructed from polystyrene scintillator bars. The scintillator layers are perpendicular to the beam's direction and alternating layers are orientated orthogonal to each other. The FGD is composed of 86.1% carbon, 7.4% hydrogen, 3.7% oxygen, 1.7% titanium, 1% silicon and 0.1% nitrogen by mass. The active region of FGD1 consists of scintillator layers only, whereas the downstream FGD (FGD2) has alternating layers of scintillator and water. The drift gas mixture used in the TPCs is Ar:CF 4 :i C 4 H 10 (95:3:2). The TPCs (TPC1 the most upstream, TPC2 the central and TPC3 the most downstream) provide excellent particle identification and accurate measurement of momentum. Together the TPCs and FGDs form the tracker region of ND280. The ECal surround the tracker and consists of 13 modules made up of plastic scintillator bars alternating with lead sheets. SMRD consists of 440 modules of plastic scintillator counters.

III. DATA AND MONTE CARLO SAMPLES
The studies reported in this paper use 5.80 × 10 20 p.o.t forward horn-current (ν-mode) data and 6.27×10 20 p.o.t of reverse horn-current (ν-mode) data broken into run periods shown in Table I. The MC simulation used for this analysis consist of a sample corresponding to ten times the data p.o.t. It is performed generating (anti)neutrino interactions according with the flux predicted at ND280. The simulation of the ν µ andν µ fluxes reaching the near detector are described in detail in Ref. [51]. The neutrino and antineutrino interactions in the ND280 sub-detectors, as well as events inside the magnet yoke and in the rock surrounding the ND280 pit, were simulated using the Neut MC generator version 5.3.2 [52]. The CCQE neutrino-nucleon cross section is simulated according to the Llewellyn-Smith formalism [53] with a dipole axial form factor and BBBA05 vector form factors [54]. The nuclear model uses a spectral function (SF), developed in Ref. [55] with an axial mass M QE A = 1.21 GeV/c 2 based on the K2K measurement of the ν µ CCQE cross section [9]. It utilizes the multi-nucleon interaction model (2p2h) from Nieves et al. [56] to simulate interactions with nucleon pairs. The model for resonant pion production (RES) is based on the Rein-Sehgal model [57] with updated nucleon form factors [58] and an invariant hadronic mass W ≤ 2 GeV. The DIS interaction is calculated for W > 1.3 GeV, using GRV98 parton distribution functions [59] with Bodek-Yang corrections [60]. Single pion production through DIS is suppressed for W ≤ 2 GeV to avoid double counting with RES and it uses a custom hadronization model. For values of the invariant hadronic mass W > 2 GeV, Pythia/JetSet [61] is used for hadronization. FSI, i.e. interactions of the hadrons produced by neutrino interactions with the other nucleons before leaving the nuclear environment, are simulated using a semiclassical intranuclear cascade model [62,63].
The propagation of the final state particles through the ND280 sub-detectors is simulated using the package Geant4 version 4.9.4 [64] as detailed in Ref. [1] employing the following physics lists: QGSP BERT for the hadronic physics, emstandard opt3 for the electromagnetic physics and G4DecayPhysics for the particle decays.

IV. ANALYSIS STRATEGY
A joint measurement of neutrino and antineutrino cross-sections, fully accounting for correlations in the systematic uncertainties, has been performed. Given the relatively large background of neutrino interactions in the antineutrino sample, such a joint analysis is mandatory for a robust antineutrino cross-section measurement. Indeed, since the neutrino and antineutrino cross-sections are largely driven by the same underlying physics, it would be inconsistent to assume to know the former while measuring the latter. A further advantage of a joint measurement, is that it exploits the full, high-statistics, neutrino sample minimizing the correlated detector and flux systematic uncertainties and thus resulting in a more precise antineutrino measurement. Finally, a joint analysis enables interesting measurements, as explained in Section I.
An unregularized binned likelihood fit with control sample to constrain the background is performed as in Ref. [15,17,65]. This analysis method guarantees a negligible dependence on the signal model used in the simulation for the correction of detector effects, provided that a too coarse binning is not used. A simultaneous fit is applied to the antineutrino sample and the neutrino sample, the former being further sub-divided in different signal and background samples depending on the direction of the outgoing muon, while the latter depending on the kinematics of the outgoing muon and proton. The number of selected events in each bin of reconstructed kinematics (j) for each signal and background sample (s) is computed as where N MC i is the true number of events in MC with the superscript indicating which interaction type they correspond to. The index i runs over the bins of the "true" muon kinematics prior to detector smearing effects, k runs over the possible background reactions and n runs over the neutrino or antineutrino energy bins. Both c i νµ and c iν µ are the parameters of interest which adjust the ν µ andν µ CC-0π number of events in MC, in order to match the observed number of events in data. The transfer matrix t det ij , relates the true (i) and reconstructed (j) muon kinematics bins and d j represents the nuisance parameters in the fit describing the detector systematics which are constrained by a prior covariance matrix. The flux parameters f n and weights w i n , describe the neutrino energy distribution for each bin of p true µ , cos θ true µ . The f n are nuisance parameters in the fit constrained by a prior covariance matrix. The product x,b runs over the systematics related to the theoretical modeling of the interaction channels contributing to the signal (x) or the background (b). Each w(x) νµCC-0π i and w(b) k i term is a weighting function describing how the generated muon kinematics change (in bins i for each signal and background process) as a function of the value of a particular theoretical parameter. All the parameters x and b are nuisance parameters in the fit and are constrained by a prior covariance matrix. Signal modeling parameters x are not fitted to avoid model dependence but they must be included to account for their effect on the uncertainty of the efficiency corrections. The parameters of interest c i νµ reweights the neutrino signal in neutrino mode and neutrino background in the antineutrino mode, whereas c iν µ reweights the antineutrino background in neutrino mode and antineutrino signal in antineutrino mode. The nuisance parameters may be different in each sample and their correlations between samples are fully taken into account.
The best fit parameters are those that minimize the following log-likelihood: where N νµ j (Nν µ j ) is the expected total number of events in the neutrino (antineutrino) sample and N νµ obs j (Nν µ obs j ) is the observed number of events. χ 2 syst is a penalty term for the systematics, where p are the parameters that describe the effect of nuisance parameters, p prior are the prior values of these systematic parameters and V syst cov is their covariance matrix which describes the confidence in the nominal parameter values, as well as, correlations between them.
To minimize the dependence of the results on the signal model used in the simulation, two-dimensional differential cross-sections are extracted as a function of muon momentum and angle. Those are kinematic quantities directly observable in the detector and they represent all the relevant variables to characterize the detector acceptance and efficiency. The signal is defined by the absence of pions in the final state, avoiding model-dependent corrections for pion re-absorption in the nucleus.
The flux-integrated cross-sections are evaluated per nucleon and for each bin i of detector unsmeared muon momentum and angle: where N νµCC-0π i and Nν µCC-0π i are the number of neutrino and antineutrino CC-0π events respectively evaluated by the fit, νµ i and νµ i are the efficiency evaluated from MC, N FV nucleons is the number of target nucleons in the fiducial volume, Φ νµ and Φν µ are the integrated fluxes for neutrino and antineutrino, ∆p µ and ∆ cos θ µ are the bin widths of the muon momentum and cosine of the muon scattering angle w.r.t. the bean direction. The number of nucleons, computed using the areal density of the different elements composing the fiducial volume [47], is equal to 5.9×10 29 and it is used to extract both cross sections. The cross sections are normalized in all bins of muon kinematics with the same integrated flux to avoid a model-dependent mapping of such bins into energy intervals of the incoming neutrino.
The binning of the true muon kinematics has been optimized to reduce the bin-by-bin fluctuation derived by the extrapolation of an unsmeared quantity, as the crosssection is, and also to ensure that the systematic uncertainty are smaller than the statistical uncertainty. If the binning is too coarse, the results do not give much information about the shape of the cross section, while on the other hand if the binning is too fine, some bins could be empty causing problems with the minimization algorithm. The best binning lies in between these extreme cases and requires that the bin width is always greater than the resolution of the muon kinematics. A MC sample simulated using the version of Neut described in Section III has been used as the prior of the fitting algorithm. This choice does not introduce model dependencies as extensively demonstrated in previous analyses [17,65]. The stability of the results has been confirmed by using alternative models in the fitting framework. To this end, a set of mock data samples has been created by modifying the amount of 2p2h interactions, the nuclear or the background model, and the input MC. Through them, it has been verified that the extracted cross section is always in agreement, within the uncertainties, with the mock data set predicted cross section and also produces a small χ 2 computed considering the final cross section covariance matrix.

A. Event selections
The event selections developed for this analysis aim to select ν µ andν µ CC-0π interactions in the FGD1 and to provide appropriate control samples to constrain the main background sources.
In previous analyses, the selection criteria were optimized to select forward going muons (with respect to the beam direction) originating from FGDs [14,15,[66][67][68]. For this analysis, the phase space of the muon kinematics was enlarged, including also high-angle and backwardgoing tracks. The acceptance has been increased using all the ND280 sub-detectors and the time of flight (ToF) of the particles between different sub-detectors which gives information about the direction of the track, i.e. if it is forward or backward with respect to the beam direction, following the same strategy described in Ref. [65].
In addition to the common goal of enlarging the acceptance, the event selections have several common features: • The selection criteria have been optimized by employing a MC sample simulated using the version of Neut described in Sec. III; • Particles that enter the TPCs or are fully contained in FGD1 are identified through the TPC or FGD particle identification (PID), based on dE/dx measurements; • ECal PID is performed if there is an associated ECal segment, which reduces the shower-like contamination (mostly π 0 ); • The ratio between the track length and the electromagnetic energy associated with the track is used to reduce the proton contamination; • Particles stopping in the SMRD are identified as muons, since most likely this is the only particle that will reach this detector.
Each selection applies a set of cuts which have been optimized to give the best signal efficiency and purity. Two requirements are common to both selections: • Events must occur within the time window of one of the eight beam bunches of the spill structure of the beam and when all ND280 sub-detectors are functioning correctly; • The interaction vertex, defined as the starting position of the muon candidate, must be inside the FGD1 fiducial volume (FV). Compared with the previous analyses where both a true and a reconstructed vertex in the first two scintillator layers were rejected [15,66], in this analysis the full span has been taken as the FV. Depending on the direction of the muon, the events with a reconstructed vertex in the first (forward-going muon) or the last (backward-going muon) layer have been rejected.
In the following sections, the selection strategy is discussed in detail.

νµ CC event selection
The selection described in this section is an improved version of the one used in Analysis I in Ref. [15], and similar to that detailed in Ref. [17] where it has been extensively described.
The target for ν µ interactions is FGD1. This is used also as a tracker with TPC1, TPC2, ECal and SMRD. After the first requirements on the data quality and the position of the vertex are fulfilled, the selection requires tracks with a TPC segment with good reconstruction quality. For such tracks, the negatively charged one with the highest momentum, and compatible with the muon hypothesis according to the TPC PID is identified as a muon candidate. Tracks fully contained in the FGD and compatible with the energy loss by a muon have also been selected.
Protons are selected by looking for a track which starts in the FGD1 FV. The track should be identified as a positively charged in a TPC, and passes both the TPC track quality cut and PID criteria. Alternatively, if the track stops within the FGD it is identified as a proton if the track is consistent with the FGD proton hypothesis. To ensure the cross section is fully inclusive in terms of numbers of protons, events without a reconstructed proton are also included. Proton selection helps in further enlarging the phase space to high-angle and backward muons, as shown in Analysis I of Ref. [15]. The selected events are divided into five signal samples: Sample I: characterized by events with only a muon candidate in one of the TPCs (TPC2 if the muon is going forward and TPC1 if it is going backward), Sample II: a muon candidate in one of the TPCs and one proton candidate in TPC2, Sample III: a muon candidate in one of the TPCs and a proton candidate in FGD1, Sample IV: a muon candidate in FGD1 and one proton in TPC2; Sample V: only a muon candidate in FGD1 that reaches the ECal or SMRD.
Events with a muon candidate in FGD or TPC and more than one proton in the final state, with the leading proton in TPC, have been selected as well. As these events only accounts for 0.8%, they have been added to the signal samples II-IV, accordingly with the muon candidate position (track in FGD only or in TPC). Fig. 3 summarizes the signal samples described above. The kinematics of the muon candidate in each sample for the CC-0π signal and the various backgrounds are  shown in Fig. 4 where the MC is broken down by true topologies. The selection is dominated by events with one reconstructed muon and no other tracks. The signal samples where the muon is reconstructed in the TPC have very similar momentum distributions, although events with a reconstructed proton tend to have muons at slightly larger angles. The sample with the muon in the FGD and the proton in the TPC have muons with much smaller momenta and larger angles. The ν µ CC-0π cross section is extracted by adding together the contributions from all the samples, but it is important to keep the events with and without protons and with the muon in different sub-detectors separated in the analysis because they are affected by different systematics and backgrounds.
The main background arises from charged-current events with one true positively charged pion (CC-1π + ), or any number of true pions (CC-Other) which are misidentified or not reconstructed, neutral current interactions (NC) and interactions that occurred outside the FV (out FV) but were reconstructed inside consti-tute a smaller background. Two control samples were selected to constrain charged current event rates with single-pion and multiple-pion production: the CC-1π + is made up of events with exactly two tracks, one negatively charged muon and one positively charged pion, and the CC-Other, made of events with more than one pion in the final state. Pions have been identified in different ways according to their charge. A π + can be identified by looking at the curvature of the track in the TPC and by requiring that the energy loss in this detector is consistent with a pion. π − are only identified by looking at the curvature of the tracks while π 0 are identified by looking for tracks in the TPC with charge depositions consistent with an electron from a γ conversion. The kinematic distributions of the control samples are shown in Fig. 5. The data-MC discrepancy in the kinematic distributions of the CC-1π + control sample was already observed in previous analysis [15,17] and is corrected by the likelihood fit as shown in Section V.

2.νµ CC event selection
After the first two common requirements described at the beginning of the Sec. IV A are fulfilled, the events are divided in four samples depending on the length of the muon candidate track in the TPCs and its direction, following the same strategy described in a recent T2K publication [65]: • If the muon candidate travels forward w.r.t. the beam direction and the associated track has more than 18 hits in TPC2 then the event belongs to the forward (FWD) sample; • If it travels backward and the associated track has more than 18 hits in TPC1 then the event belongs to the backward (BWD) sample; • If the muon candidate travels forward but the track has less than 19 hits in TPC2 then the event belongs to the high angle forward (HAFWD) sample; • if it travels backward and the associated track has less than 19 hits in TPC1 then the event belongs to the high angle backward (HABWD) sample.
For each sample, different sets of cuts have been developed to reduce the background as much as possible without decreasing the efficiency.
One of the main backgrounds is caused by interactions that occur outside the FV (out FV) but are incorrectly reconstructed as starting inside the FV. This can be due to a failure of the reconstruction algorithms or a scattering of the particle which can lead to two unmatched track segments, one of which may start in the FV. The ratio of the momentum of the muon candidate to the other unmatched track and also the minimal distance between the tracks are used to reduce this background. The ratio should be lower than 1 if the two segments belong to the same track. Different cut values have been chosen for the event falling in the FWD, HAFWD and HABWD samples. These cuts are not applied in the selection of the BWD sample since signal events could be rejected.
Another misreconstruction pathology can break a single track into two segments, with a reconstructed vertex inside the FV and a forward-going track into the downstream TPC. This often happens near the downstream edge of the FGD and the second track is considered as the muon candidate. Therefore, events for which the start position of the track associated to the muon candidate is in one of the last two layers of FGD1 are rejected.
The muon candidate is identified as the highest momentum track that is consistent with the muon PID. If the muon candidate enters a TPC, the track must pass the TPC muon PID. If the track does not enter a TPC, the ECal portion of the reconstructed object must be consistent with the ECal muon PID. In the case where the muon candidate enters a TPC, the charge of the track will be included in the selection. For particles entering ECal the information on the charge is not available, therefore this sample of events presents a high contamination of negatively charged muons that is constrained by measuring at the same time the ν µ cross section.
In the selection of the FWD sample two additional cuts have been applied to reduce the pion and proton contamination: if the muon candidate stops in FGD2 and has a momentum greater than 280 MeV/c, the candidate is most likely a pion or a proton and the event is rejected; if it reaches the ECal it must have an ECal PID compatible with the muon hypothesis.
The described cuts select samples of muon antineutrino CC events with muons in every direction. Every sample is then split in three sub-samples according to the event pion multiplicity: events without a reconstructed pion, with one negatively charged pion, or with more than one pion in the final state. The TPC pion selection is similar to the one described previously in Section IV A 1. If the pion-candidate track is contained in the FGD1, pions are identified in two ways: by requiring a charge deposition consistent with a pion, or using the delayed energy deposition in the FGD due to a decay electron coming from π → µ decay. In the latter case, the pion is tagged as positively charged since the π − are more likely to be absorbed.
The kinematics of the muon candidate for the signal   constituting an irreducible background. In the high-angle selection, the charge is not reconstructed, therefore negatively charged muons are also selected (Section IV). The statistics of the BWD sample is limited as the antineutrino cross section is suppressed for backwards going muons. Fig. 6 summarizes theν µ signal samples.
The background mostly arises from events with one true negatively charged pion (CC-1π − ), any number of true pions (CC-Other) and out FV, that contributes more to the BWD and HABWD samples. The CC-1π − and CC-Other samples identified through the pion tagging are employed to constrain such backgrounds. For the out FV background there is not a dedicated control sample. The majority of them are ν µ CC interactions that are constrained by the existing control samples. An uncertainty on the prediction of the out FV interactions is taken into account as reported in Section IV B. The kinematic distributions of the control samples are shown in Fig. 8. As shown in the legend, the purity is lower than the ν µ selection, at 48% for the CC-1π − sample and 24% for the CC-Other sample. Indeed, positively charged pi-ons generated in ν µ interactions are mis-identified as positively charged muons decreasing the purity. This difference with the ν µ selection is caused by the higher ν µ contamination of the antineutrino beam compared with the smallerν µ component in the neutrino beam. The data-MC discrepancy observed for the CC-1π − control sample is mainly due to an overestimation of the antineutrino coherent pion production cross-section as is implemented in NEUT version 5.3.2 [69]. Also in this case, the discrepancy is corrected by the likelihood fit (see Section V).

B. Sources of uncertainties
The uncertainties can be split into the following categories: statistical uncertainty, flux uncertainty, detector systematic uncertainties, uncertainty on the modeling of signal and background processes.
Statistical Uncertainty. To compute the data statistical uncertainty, the nominal MC was normalized to the number of protons on target in the data. Then 1,000 toy samples were generated where the MC was varied in each reconstructed bin according to a Poisson distribution and the fit was performed on each toy. The statistical error is taken as the width of the variation induced on the cross sections distribution in each true bin.
Flux Uncertainty. The evaluation of the uncertainties on the flux prediction are described in detail in Ref. [51]. It is around 10% at the energy peak and is dominated by the hadron production model and is evaluated using data published by the NA61/SHINE experiment using a thin Carbon target [70][71][72]. The flux covariance matrix was used to generate many toy MC sets. The flux bins include separate bins for the "right-sign" and "wrong-sign" components of the flux in both neutrinomode and antineutrino-mode. The fit includes 32 nuisance parameters for the fluxes which are constrained by the fit, reducing the flux uncertainties by around 60%. In previous analyses the flux was not constrained by the fit since this procedure could introduce a model dependencies [15]. For this reason dedicated mock data studies have been performed as described in Section IV.
Detector Systematic Uncertainties. Detector uncertainties can be grouped into three categories depending on the way they are propagated: efficiency-like, observable variation and normalization systematics. The systematics belonging to the first group have been propagated by applying a weight that depends on one or more observables, the second by adjusting the observables and reapplying the selection, the last by applying a single weight applied to all events. Dedicated data and MC samples have been used to quantify the detector uncertainties in the modeling of FGD and TPC responses, of neutrino interactions outside of the FGD1 FV, pion and proton secondary interactions. The differences between data and MC observed in control samples have been applied as correction factors to the nominal MC to take into account the observed discrepancies, while the error on these factors has been taken as detector systematic uncertainty. The ν µ andν µ selections are affected by the same detector uncertainties since they employ similar features of the sub-detectors. The dominant systematics are due to the uncertainties on the amount of background from interactions occurring outside of the fiducial volume, the modeling of the pion secondary interactions and the TPC PID. The detector systematics have been stored in a covariance matrix corresponding to the uncertainties on the total number of reconstructed events in each bin and in each signal and control regions. The systematics on the cross sections have been propagated by repeating the fit over many toy MC data sets where the detector parameters have been varied according with their covariance matrix but have been kept fixed in the fit. This choice has been driven by the necessity to ensure the stability and convergence of the fit by reducing the number of nuisance parameters. The uncertainty associated with the number of targets has been computed separately. A covariance matrix has been formed from the uncertainties on the areal densities of all the elements present in FGD1 and, the uncertainty calculated from toy experiments sampling such covariance matrix.

Parameter
Prior Error Modeling of Signal and Background. The signal efficiency and number of background events in each bin are affected by uncertainties in our cross-section model. Table II summarizes         ization of such process; NC Coherent and NC Other the normalization of NC interactions. Pions that are produced in neutrino interactions can be affected by FSI as they leave the nuclear medium, changing their kinematics, charge and multiplicity. Dedicated systematic parameters have been included in the fit to describe the pion production, absorption, charge exchange and quasielastic scattering of the exiting pions. Again these modify not only the selected number of events, but the selection efficiency as well. Similarly protons are also subjected to FSI: the uncertainty is evaluated by comparing two dif-ferent NuWro [73] MC simulations 1 , with and without FSI. The difference in the efficiency as a function of the muon kinematics between the two simulations has been taken as the uncertainty due to the proton FSI. To be conservative, it has been added in quadrature to the other efficiency uncertainties. Since inν µ interactions proton FSI has a negligible impact, this uncertainty has been added to ν µ cross section only. The cross section parameters have been propagated by throwing from Gaussian distributions that have as mean and sigma the prior and error values reported in Table II. All of the systematic errors in each bin are summarized in Figs. 9 and 10. In most bins, the dominant uncertainty is due to the statistical error on data. The systematic uncertainties are typically dominated by the detector systematics. The modeling errors are generally subdominant and smaller than 10% and closer to 1% in regions of high purity. Uncertainties related to the flux, modeling of the background and pion FSI have been propagated together since they are anti-correlated. The errors on the sum, difference and asymmetry has been computed numerically from toy experiments sampling the covariance matrix that includes the uncertainty and correlations between the two cross sections.

V. RESULTS AND COMPARISONS WITH MODELS
The distribution of reconstructed events in bins used to evaluate the cross sections and in the background control samples is shown in Figs. 11 and 12. The data are compared to MC predictions before and after the fit. The fit is able to reproduce the observed distributions in data by varying the parameters of interest to describe the signal cross section and the nuisance parameters describing the systematic uncertainties, as explained in Section IV. As Figs. 13 and 14 show, the large discrepancy in the prefit MC prediction in the CC1π ± control region is well corrected by the fit by varying the nuisance parameters describing the pion production cross-section listed in Table II.
In the following, the measured cross sections, and their combinations, are compared to different (anti)neutrinointeraction models using the framework Nuisance [74] and the agreement is quantified by the χ 2 statistic. Since the ν µ and theν µ cross sections are extracted simultaneously, a global χ 2 computed using the full covariance matrix, i.e. the one containing the correlation between the two cross sections, is reported. It should be noted that, especially for the neutrino and antineutrino cross sections, as well as their sum, the overall normalization uncertainty (fully correlated between bins) constitutes a relatively large fraction of the uncertainty. In particular it contributes to 48% and 35% of the total systematic uncertainty for neutrino and antineutrino cross sections respectively and to 49% for the sum, while for the difference it decreases to 19% and for the asymmetry to 5%. Therefore the χ 2 statistics may suffer from Peelle's Pertinent Puzzle [75] and may not be a reliable estimation of the data-MC agreement. This issue does not affect the shape-only χ 2 which is reported as well.
The models considered for the data-MC comparisons are as follows: • NuWro Spectral Function (SF), as developed in Ref. [55], using the same 2p2h model as Neut; • GiBUU 2019 LFG in a coordinate-and momentum-dependent nuclear potential, as described in Ref. [76], using a 2p2h model based on Ref. [77] and further tuned in Ref. [78], which uses the T2K measurements of final-state muon and proton kinematics and correlations in charged-current pionless interactions discussed in Ref. [17]; • SuSav2 is a complete implementation of the SuSAv2 model [79][80][81][82] in GENIE, as described in [83], where 1p1h is based on the Relativistic Mean Field approach [84] and 2p2h is based on the calculation from Ref. [85]. The pion production and FSI models are the same as in Genie version 3.00.04, configuration G18 10b 000 00; The comparisons with the models described above are shown in Figs. 15 to 34. The full and shape-only χ 2 are reported in the legends (shape-only χ 2 is reported in parenthesis) and are summarized in Table III. In Table IV the reduced χ 2 is reported as well.
In order to evaluate the sensitivity to the 2p2h process, the measured cross sections, and their combinations, are compared in Figs. 15 to 19 to Neut LFG with and without 2p2h. The full and shape-only χ 2 show that the sensitivity is limited. Some conclusions can be drawn looking at each angular bin. In the intermediate and high-angle region both neutrino and antineutrino data tend to prefer the presence of 2p2h, as already shown in the previous T2K neutrino analysis [15]. The χ 2 in each angular bin has been computed, further confirming the preference for the presence of 2p2h in the intermediate and high-angle region. The effect is particularly evident in the sum of the neutrino and antineutrino cross sections, where the statistical uncertainty is smaller. For instance, in the angular bin 0.6 < cos θ µ < 0.7 the reduced χ 2 is 0.8 and 2.4 with and without 2p2h respectively. On the other hand, a clear overestimation of the cross section is visible in the forward region below 1 GeV, both for neutrinos and antineutrinos. This may point to incorrect 1p1h predictions, notably in the region of small energy transfer to the nucleus, where the treatment of various nuclear effects, like binding energy, is crucial. This issue is further discussed below, in the comparison to different 1p1h models. As expected the neutrino-antineutrino cross section difference emphasizes the 2p2h cross section, due to the change of sign of the axial-vector component. The statistical and systematic uncertainties, which are dominated by the flux, are still too large for a measurement of this component. Future foreseen reduction of such uncertainties, with more ND280 data and relying on NA61/SHINE T2K replica target data for flux tuning [86], will improve the sensitivity to the axial-vector 2p2h component. In some bins the difference is negative since antineutrinos can interact with the hydrogen of the hydrocarbon molecule, leading to a cross-section for antineutrino higher than for neutrino. The neutrinoantineutrino cross-section asymmetry shows a very small 2p2h dependence. The fractional change of the asymmetry with and without 2p2h is very small, except in the low momentum region where, at forward angle, it may reach 50%. The sensitivity to such observable is drastically limited by the statistical uncertainty. Despite most of the systematic uncertainties cancel out a residual detector systematic not correlated between neutrino and antineutrino dominate the systematic error.
A more sophisticated assessment of the 2p2h sensitivity is shown in Figs. 20 to 24, where the results are compared to different 2p2h models. The 2p2h model in Neut and the 2p2h model by Martini et al. [31] are both implemented on top of a similar 1p1h LFG model while the SuSav2 model includes different 1p1h [84] and 2p2h [85] predictions. For the comparison with the model from Martini et al. the number of degrees of freedom (ndof) have been reduced to 96 for the cross sections and to 48 for their combinations in terms of sum, difference and asymmetry (w.r.t. 116 and 58 respectively) because the model predicts the cross section only for muon momentum lower than 3 GeV/c. Thus 10 high-momentum bins have been removed from the covariance matrix to compute the full χ 2 . Similarly, a complete shape-only covariance matrix has been obtained and those 10 bins have been removed afterwards to compute the shape-only χ 2 . An extended implementation of this model would be crucial for a better comparison with other models. None of the model is able to well describe the measured neutrino and antineutrino cross-sections in the entire phase space. As previously mentioned, the disagreement with crosssection measurements can be interpreted both in terms of 1p1h or 2p2h mismodeling. On the other hand, the various 2p2h models have quite different predictions for the axial-vector component, making the measurement of the neutrino-antineutrino cross-section difference a powerful probe to test the physics implemented in the different 2p2h predictions.
To further investigate the dependence of the results on the 1p1h model, the measured cross-sections, and their combinations, are compared to different LFG implementations in Figs. 25 to 29. The Neut, NuWro and Genie LFG implementations differ mainly in the treatment of the nucleon binding energy. None of the generators is able to describe the measured neutrino and antineutrino cross-sections in the entire phase space. Among the different combinations the cross-sections difference show the lowest full χ 2 in the comparison with Genie.
The measured cross-sections, and their combinations, are also compared to a SF model in Figs. 30 to 34. The SF cross-section shows a different angular dependence than the LFG one: smaller for the backward and high-angle region and larger in the forward region. Interestingly, while SF is a more sophisticated model, the full χ 2 is the largest (see Table III). This may be due to an incomplete implementation of SF or to the merging with a 2p2h simulation modeled using RFG as nuclear model. The difference between LFG and SF tends to cancel in the neutrino-antineutrino cross-section difference and asymmetry. A more complete implementation of an SF model (including a 2p2h contribution) is likely required to investigate this further.
The integrated cross sections per nucleon and their combinations are reported in Table V and compared with the model described above. It is striking that the models which exhibit best agreement in shape and in normalization are different, calling for further measurements with smaller systematic uncertainties and further model development.
In summary, even if some conclusion can be drawn looking at the comparisons in some angular bins, none of the models is able to simultaneously describe ν µ and ν µ CC-0π cross sections in all the phase space. Among the different combination, the difference between neutrino and antineutrino cross sections shows interesting sensitivity to different 2p2h models, which is limited by large uncertainties.
The poor (anti)neutrino-nucleus interaction modeling highlighted by this analysis is a limiting factor for the future neutrino oscillation experiments that have as primary goal the measurement of the CP violation, calling for a deeper understanding of the underlying processes involved in (anti)neutrino-nucleus interactions and for new cross-section analyses with larger statistics and improved systematic uncertainties.

VI. CONCLUSIONS
The T2K experiment has measured the first combined double-differential ν µ andν µ cross sections with no pions in the final state in the full phase space using 5.8 × 10 20 POT of neutrino data and 6.2 × 10 20 POT of antineu-trino data. The inclusion of ToF, in the selection of backward-going and high-angle tracks, enabled the exploration of the full phase space with better efficiency over previously reported T2K measurements of neutrino cross sections [15]. The sum, difference and asymmetry of neutrino and antineutrino cross sections were measured, including full treatment of the correlations between the neutrino and antineutrino samples. Such observables have been compared with different models to shed light on the nuclear effects involved in the (anti)neutrinonucleus interactions. Although none of the models considered in this work are able to describe the full phase space of the neutrino and antineutrino CC-0π cross section, it is difficult to determine the source of the problem. A precise understanding of this mis-modelling may be of critical importance for the next generation of neutrino oscillation experiments. Further investigation would benefit from smaller uncertainties and a mitigation of some of the approximations built into generator implementations of the models.
This analysis opens the road to joint cross-section measurements putting together different samples to minimize systematic uncertainties and to account properly for correlated systematics, enabling more complete and precise tuning of neutrino-nucleus interactions. A promising observable, measured here for the first time, is the difference between neutrino and antineutrino cross sections which shows interesting sensitivity to different 2p2h models, that can be further explored with more statistics and improved systematics uncertainties.
The data release for the results presented in this analysis is posted at the link in Ref. [87]. It contains the ν µ andν µ double-differential cross sections central values, their combinations and associated covariance matrices.

VII. ACKNOWLEDGEMENTS
We thank the J-PARC staff for superb accelerator performance. We thank the CERN NA61/SHINE Collaboration for providing valuable particle production data. We acknowledge the support of MEXT, Japan;                      FIG. 34. Measured double-differential CC-0π cross-section asymmetry in bins of true muon kinematics with systematic uncertainty (red bars) and total (stat.+syst.) uncertainty (black bars). The result is compared with NuWro version 18.02.1 with LFG+RPA (green solid line) and with the SF nuclear model (green dashed line), both including 2p2h predictions. The full and shape-only (in parenthesis) χ2 are reported. The last bin in momentum is not displayed for readability.