First measurement of muon neutrino charged-current interactions on hydrocarbon without pions in the final state using multiple detectors with correlated energy spectra at T2K

This paper reports the first measurement of muon neutrino charged-current interactions without pions in the final state using multiple detectors with correlated energy spectra at T2K. The data was collected on hydrocarbon targets using the off-axis T2K near detector (ND280) and the on-axis T2K near detector (INGRID) with neutrino energy spectra peaked at 0.6 GeV and 1.1 GeV respectively. The correlated neutrino flux presents an opportunity to reduce the impact of the flux uncertainty and to study the energy dependence of neutrino interactions. The extracted double-differential cross sections are compared to several Monte Carlo neutrino-nucleus interaction event generators showing the agreement between both detectors individually and with the correlated result.


I. INTRODUCTION
In the last decade, neutrino oscillations experiments have continued to collect more statistics and reduce their systematic uncertainties, moving forward to the precision era of neutrino oscillation physics [1][2][3][4][5][6][7].For this reason, the next generation of long-baseline (LBL) neutrino oscillations experiments, DUNE (Deep Underground Neutrino Experiment) [8] and Hyper-Kamiokande [9], require systematic errors reduced to a few percent to achieve their physics goals, including precise measurements of the neutrino mass hierarchy and leptonic CP-violation * also at Université Paris-Saclay † 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 Moscow Institute of Physics and Technology (MIPT), Moscow region, Russia and National Research Nuclear University "MEPhI", Moscow, Russia ‖ also at IPSA-DRII, France * * also at the Graduate University of Science and Technology, Vietnam Academy of Science and Technology † † 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. [10,11].In order to reach this unprecedented reduction of systematic uncertainties, our knowledge of neutrinonucleus interaction cross sections must be improved.In the range of energies used in current LBL neutrino oscillation experiments, a precise knowledge of neutrino interactions with nucleons is crucial for the extrapolation from the near to the far detector.Incorrect modeling of neutrino interactions can affect the reconstructed neutrino energy, which can introduce bias in the measurement of neutrino oscillation parameters.This reduction of interaction uncertainty in part will be accomplished through measurements using the planned near detectors for DUNE and Hyper-Kamiokande, but can also be achieved through performing measurements and exploring new techniques with current generation neutrino experiments.
Neutrino charged-current quasielastic (CCQE) interactions, also referred to as one-particle-one-hole (1p1h) excitations, can be written as: where ν ℓ is the incident neutrino of flavor ℓ, n and p are the struck neutron and outgoing proton respectively, and ℓ is the charged lepton.CCQE interactions are the dominant reaction at the T2K (Tōkai-to-Kamioka) neutrino beam energy (peaked at 0.6 GeV).Modelling these interactions with a bound nucleon inside a nucleus is complex, requiring treatment of the Fermi motion, removal energy, and nucleon-nucleon correlations.Various models exist for predicting the initial state nucleon momentum and removal energy, such as Fermi gas models or relativistic mean field models, and for modelling correlations between nucleons, for example the random phase approximation method [12][13][14][15][16]. Interactions with correlated pairs of nucleons, referred to as multi-nucleon or two-particle-two-hole (2p2h) excitations, are possible due to meson-exchange currents or short range correlations in the nucleus [16][17][18][19][20][21][22].These multi-nucleon interactions enhance the neutrino cross section in the energy range of T2K and can easily be confused for CCQE interactions, which can then bias the oscillation analysis if not considered.The global picture of neutrino cross-section data is still complicated as many results are in tension with each other, and the available models and Monte Carlo (MC) generators cannot accurately describe many different results across experiments [23,24].In recent T2K oscillation results [25], the dominant systematic uncertainty is from the nucleon removal energy on charged current quasielastic interactions, showcasing the need for further study of neutrino cross sections and cross-section modelling.
The near detectors (close to the neutrino source) used by T2K provide a unique opportunity to perform a combined measurement using the same neutrino beam with two detectors exposed to different but correlated spectra of incident neutrinos, and is the subject of the analysis presented in this paper.Neutrino detectors measure the rate of neutrino interactions, which is primarily a product of the neutrino flux and neutrino cross section.Changes in both the flux and cross-section models can cause the observed event rate to change, often in similar ways, and this degeneracy limits the ability to separate individual effects due to either the flux or cross section.The correlation between the different fluxes at the near detectors provides additional information that can be used to constrain the flux uncertainty and break some of the degeneracy between flux and cross section.The different neutrino energy spectra seen at each detector also presents an opportunity to study the energy dependence of neutrino interactions within the same analysis framework.An example is the energy dependence of multi-nucleon interactions, which comprise a non-negligible fraction of the samples used to measure the cross section presented in this paper (on average 10% across all samples).The multi-nucleon cross-section prediction as a function of energy for the Nieves et al. model [20], which is the default multi-nucleon model used in T2K, and the Martini et al. model [16] shows differences mainly related to normalization of by about a factor of two to three across the neutrino energy range used at T2K (shown in Fig. 1 from [26]).This variation between models motivated an additional systematic uncertainty for the T2K oscillation analysis [1].The analysis presented in this paper takes advantage of the T2K near detector setup to perform the first measurement using multiple detectors with different neutrino energy spectra.

II. THE T2K EXPERIMENT
T2K [27] is a second-generation long-baseline neutrino oscillation experiment based in Japan, which is able to measure neutrino oscillations with a ν µ (ν µ ) beam.The neutrino beam is produced at the Japan Proton Accelerator Research Complex (J-PARC) in Tōkai.It is first detected 280 m downstream from the source at the near detector complex, where the flavor composition of the incoming neutrino flux is not expected to be affected by oscillations, and then it travels 295 km to the Super-Kamiokande (SK) far detector [28,29], located in Hida, where oscillations significantly affect the flavor composition.The near detector complex houses the two detectors of primary interest for the analysis presented in this paper: a detector on the axis of the neutrino beam, called INGRID [30] (Interactive Neutrino GRID), and a detector located 2.5 degrees off-axis, called ND280 [31] (Near Detector at 280 meters).INGRID primarily serves as a neutrino beam and flux monitor, measuring the total rate of neutrino interactions and the beam direction.ND280 is dedicated to the study of the un-oscillated spectrum of neutrinos at 280 meters from the production point and neutrino interaction cross-section properties.
The Super-Kamiokande far detector is a deep underground 50 kton water Cherenkov detector.The SK detector, as with ND280, is situated at 2.5 degrees off-axis, meaning that it is exposed to the same relatively narrow energy band neutrino flux, peaked at the oscillation maximum, around 0.6 GeV.

A. Neutrino beam
T2K neutrinos come from in-flight decays of focused hadrons emitted from an extended, 91.4 cm long, monolithic graphite target.The target is bombarded with a 30 GeV proton beam produced at J-PARC.Interactions of beam protons inside the target initiate a chain of hadronic interactions, the charged products of which are focused upon exit from the target using a series of three magnetic horns.The polarity of the horn current determines whether a ν µ (neutrino mode) or νµ (antineutrino mode) enhanced beam is produced, by focusing predominantly positively or negatively charged pions and kaons, respectively.These mesons are then left to decay, e.g.via π ± → µ ± + ν µ (ν µ ), in a 96 m long decay volume, capped with a concrete beam dump at the downstream end.Behind the beam dump, a muon monitor [32,33] is used to measure the secondary beam stability.

B. INGRID
The INGRID detector is an on-axis neutrino detector located 280 m downstream of the proton target.It consists of 14 identical detector modules (referred to as standard modules) and an extra module called the Proton Module (PM).
The main purpose of the standard modules is to monitor the neutrino beam direction.The 14 identical standard modules are arranged in two identical groups along the horizontal and vertical axes, as shown in Fig. 3.Each of the modules consists of nine iron target plates and eleven tracking scintillator planes surrounded by veto scintillator planes to reject charged particles coming from outside the modules [34], as shown in Fig. 4. By contrast, the Proton Module was specifically developed for neutrino cross-section measurements.It is located at the beam center between the horizontal and vertical standard modules as shown in Fig 5 .It is a fullyactive tracking detector consisting of 36 tracking layers surrounded by veto planes (shown in Fig. 6), where each tracking layer is an array of two types of scintillator bars [35].Each scintillator plane covers an area of 120 × 120 cm 2 transverse to the beam direction.The tracking layers also serve as the neutrino interaction target, with the total target mass of the scintillator and fibers in the fiducial volume being 292.1 kg.

C. ND280
The off-axis near detector ND280 (Fig. 7), is a magnetized particle tracking device.It consists of a number of sub-detectors installed inside the refurbished UA1/NOMAD magnet, which provides a 0.2 T field used to measure the charge and momentum of particles passing through ND280.Inside the UA1 magnet, the neutrino beam first passes through the π 0 detector (P0D) [36] and then the inner tracker, both of which are surrounded by an electromagnetic calorimeter (ECal) [37].Moreover the UA1 magnet yoke is instrumented with plastic scintillator to perform as a muon range detector (SMRD) [38] in order to track high angle muons and "sand muons" coming from neutrino interactions in the rock upstream of the detector.The tracker region of ND280 consists of three Time Projection Chambers (TPC1, 2, 3) [39], interleaved with two Fine-Grained Detectors (FGD1, 2) [40].The upstream FGD1 detector is made of fifteen XY planes of polystyrene scintillator with each plane having 2 × 192 bars, while the downstream FGD2 contains seven polystyrene scintillator modules interleaved with six modules of water in between.The FGDs provide 1.1 tons target mass each for neutrino interactions and tracking of the charged particles coming from the interaction vertex, while the TPCs provide 3D tracking and determine the momentum and energy loss of each charged particle traversing them.The observed energy loss in the TPCs, combined with the measurement of the momentum, is used for particle identification (PID).The analysis presented here is focused on neutrino interactions on carbon, including only events occurring in FGD1.

III. EVENT SIMULATION AND SELECTION
The goal of this analysis was to perform a simultaneous fit to ND280 and INGRID data, extracting the muon neutrino flux-integrated differential cross section on hydrocarbon without pions in the final state as a function of the outgoing muon kinematics for both the off-and on-axis T2K flux.Signal events are defined by a neutrino interaction with an outgoing muon, zero pions, and any number of other hadrons in the final state and are referred to as CC-0π events (or topology).This signal definition is chosen because it is the most common interaction for the T2K oscillation analysis and to match what is accessible to the detectors: the outgoing finalstate particles that exit the nucleus.Particles produced in the neutrino interaction can re-interact as they leave the nucleus, potentially producing new particles or being absorbed, referred to as final-state interactions (FSI).Defining the signal in terms of the final-state particles reduces the model dependence of attempting to correct for FSI effects.Similarly, the cross section is measured as a function of the outgoing muon kinematics as opposed to using the reconstructed neutrino energy or momentum transfer to avoid as much model dependence as possible.

A. Event simulation
The T2K neutrino flux simulation [41] is based on the modeling of proton interactions with the graphite target and propagating the produced particles through the target station, allowing for further interactions.Interactions within the target are simulated using the fluka 2011 package [42,43] while out-of-target interactions and decays are handled by the geant3 [44] and gcalor [45] packages.Hadronic interactions and multiplicities are tuned using NA61/SHINE thin-target data [46][47][48] and data from other experiments [49][50][51].The proton beam conditions, horn current, and neutrino beam position are monitored and used as inputs to the flux simulation to provide additional constraints.Combined, this data-driven procedure gives an overall flux normalization prior uncertainty of about 8.5% at ND280 and 9.9% at INGRID for this analysis, which is dominated by hadron production and interaction uncertainties.The ND280 and INGRID flux predictions are produced simultaneously using the same input parameters, and this results in correlated uncertainties that are included in this analysis (and described further in Sec.IV B).
Neutrino interactions in the detectors and the outgoing kinematics of the produced final-state particles are simulated using the NEUT neutrino event generator [52,53].NEUT describes charged-current quasi-elastic (CCQE) neutrino-nucleon interactions using the spectral function (SF) approach from [54] with the quasi-elastic axial mass (M QE A ) set to 1.21 GeV/c 2 based on the K2K CCQE cross-section measurement in [55].Multi-nucleon correlations (also referred to as 2p2h interactions) are based on the model from Nieves et.al. [20].Resonant pion production (RES) is described by the Rein-Sehgal model [56] using updated nucleon form factors [57] and the resonant axial mass (M RES A ) set to 0.95 GeV/c 2 .Coherent pion production uses the updated Berger-Sehgal model [58].Deep inelastic scattering (DIS) interactions are modeled using the GRV98 parton distribution functions [59] with corrections from Bodek and Yang [60] to extend the validity of the treatment to lower fourmomentum transfer (Q 2 ≲ 1.5 GeV 2 ).NEUT begins modeling DIS processes for interactions with hadronic invariant mass W > 1.3 GeV/c 2 .For interactions with 1.3 < W < 2.0 GeV/c 2 a custom hadronization model [61] is used to interpolate between RES and DIS processes, while for W > 2.0 GeV/c 2 , pythia/jetset [62] is used for the hadronization model.Hadrons produced in the primary neutrino-nucleon interaction must propagate through the nuclear remnant before they can be detected.Interactions before the hadrons leave are referred to as final-state interactions (FSI), and are simulated using a semi-classical intra-nuclear cascade model [63,64].The MC productions for each detector use the same physics models in NEUT, but are based on slightly different versions, with ND280 using version 5.3.2 and INGRID using version 5.3.3,however this has a negligible impact on the analysis.
The propagation of the final-state particles through the detector medium after exiting the nucleus is performed using a geant4 [65] simulation.Both detector simulations use qgsp bert for the hadronic physics list [66].The detector readout simulation is handled by a custom electronics simulation separately for ND280 and INGRID [27].

B. Data samples
This analysis uses neutrino-mode data collected between 2010 and 2017 during T2K Runs 2 through 8.The ND280 sample corresponds to a total of 11.53 × 10 20 POT (protons on target), while the INGRID sample corresponds to a total of 6.04 × 10 20 POT.The breakdown of collected data by run period is listed in Tab.II.The INGRID detector configuration was changed after Run 4 where the Proton Module was moved to a different location in the detector hall, which limits the usable data for this analysis and is the main reason for the difference in total POT between the ND280 and INGRID samples.T2K Run 3b used a lower horn current (205 kA instead of 250 kA) during data-taking, and is included in the ND280 data set but excluded from the INGRID data set.This is what was done in previous ND280 [67][68][69] and INGRID analyses [70], and kept the same for this analysis for consistency and could be revisited for future versions.

C. Signal selection
The signal selection for this analysis is designed to select muon neutrino events with no detected pions in the final state and any number of visible protons, referred to as the CC-0π topology.The target material is the plastic scintillator in either FGD1 (for ND280) or the Proton Module (for INGRID).The individual selections for ND280 and INGRID were developed for previous analyses, described in Refs.[67][68][69] and Ref. [70] for ND280 and INGRID respectively, and minor updates necessary for the joint fit and the addition of new data were made for the analysis presented in this paper.

ND280
The ND280 selection first requires passing a set of data quality cuts, and then requires the interaction vertex to be within the FGD1 fiducial volume (FV).The FV is defined to include events with a vertex at least five scintillator bars from the edge in the X and Y directions, and excludes the first XY module as an upstream veto.Events with a single negatively charged muon candidate and any number of proton candidates sharing a common vertex are identified and classified into different samples based on the detectors (FGD1 or TPCs) used to measure the momentum of the muon candidate and the proton candidate(s), if any.This sample separation by detector and particle content allows for a more precise treatment of the detector systematics due to the different detector responses.Tracks are identified by their energy deposition and curvature compared to the expected distributions for each particle hypothesis.The momentum of each reconstructed track is measured either by curvature in the TPCs or by range in FGD1 (and ECAL).Events with a detected associated decay electron in FGD1 are treated as background as these events are likely to have produced an untracked stopped pion decaying into a muon followed by a Michel electron from the muon decay.The signal events are classified into the following samples: • Sample I (µTPC): defined by a single muon candidate in the TPCs and no other tracks; • Sample II (µTPC+pTPC): a muon candidate in the TPCs with one or more proton candidates in the TPC; • Sample III (µTPC+pFGD): a muon candidate in the TPCs and a proton candidate in FGD1; • Sample IV (µFGD+pTPC): a muon candidate in FGD1 (possibly reaching the ECAL) and a proton candidate in the TPC; • Sample V (µFGD): a muon candidate in FGD1 (possibly reaching the ECAL) and no other tracks.
The majority of events in the signal sample (∼ 62%) are events with a single reconstructed muon and no other tracks, with most muons reaching the TPCs.The kinematic distributions of each sample separated by true topology are shown in Figs. 8 and 9 with data plus statistical errors overlaid.Most samples achieve a high purity of CC-0π events (approximately 82% pure when integrated over all signal samples) with the main background coming from misidentified or unidentified pions from CC-1π + (events with a muon and single positive pion track) or CC-Other (events with a muon and multiple pion tracks) events.The normalization of the data and nominal MC is very similar when integrated across all the samples, but varies within 15% per sample.A noticeable feature is the slight deficit of data events compared to the nominal MC at very forward angles for the combined µTPC sample (representing ∼ 85% of the total event sample).The selected events are binned using the reconstructed muon momentum and cosine of the angle with respect to the beam direction (list of bin edges available in Appendix A).The binning scheme is designed such that bins are not finer than the detector resolution.
The cross section is extracted by adding the contributions from each sample, but the samples are kept separate in the analysis.This is important because events with and without protons and which subdetectors were used in the reconstruction are affected by different systematics and backgrounds.

INGRID
The INGRID selection first requires passing a set of data quality cuts, and then requires the interaction vertex to be within the Proton Module fiducial volume.The FV is defined to be the transverse central ±50 × ±50 cm 2 region of the Proton Module and excludes the first four scintillator layers as an upstream veto.Events with exactly one or two tracks sharing a common vertex are selected where one track must be minimum ionizingthe muon candidate -and the second track, if present, must be proton-like according to the PID.The INGRID PID algorithm is based on a boosted decision tree that uses the dE/dx along the track and the distribution of deposited energy with respect to distance from the end point of the track.Protons will tend to deposit more energy at the end of the track compared to muons or pions (referred to as a Bragg peak).The muon candidate track must either stop in the Proton Module or reach the standard module directly downstream, where it may also stop or traverse the entire module and escape.Events where the muon escapes out the side of the Proton Module are rejected.The momentum of a stopping muon candidate track is measured by calculating an equivalent distance traversed in iron, and muon candidate tracks that travel through the entire standard module and escape have a lower limit on their momentum.The momentum threshold for a muon to escape the standard module is approximately 1 GeV/c.The selected events are binned using the reconstructed momentum and angle with respect to the beam direction (list of bin edges available in Appendix A).Stopping and escaping events are considered together as a single CC-0π sample in this analysis (Sample IX).
The kinematic distribution of the signal sample separated by true topology is shown in Fig. 10 with data plus statistical errors overlaid.The INGRID sample is notably less pure than the ND280 sample, with a much higher background primarily coming from pions being misidentified as muons.For the cross-section extraction, the kinematic regions of p µ < 0.35 GeV/c and cos(θ µ ) < 0.50 are excluded from the analysis to remove regions of no acceptance due to the detector geometry.

D. Control regions
To provide a better constraint on the background contributions, a set of control samples are included in the analysis.As with the signal selections, the ND280 and INGRID control samples are designed to select similar types of events, however the additional capabilities of ND280 allow for more complicated event topologies.

ND280
The ND280 selection includes three control samples to select events with a pion, constraining the primary background contribution to the signal selection.These samples follow similar initial data quality cuts and criteria for identifying a muon candidate, but also require the identification of a pion candidate.A new addition for this analysis compared to previous ND280 analyses is the inclusion of a separate sample designed to identify low momentum pion events by detecting the presence of a Michel decay electron in FGD1.The control samples are categorized by the pion content as follows: • Sample VI (CC-1π + ): defined by a single muon candidate in the TPC and one π + candidate in the TPC; • Sample VII (CC-Other): a muon candidate, one π + candidate, and at least one additional track in the TPC; • Sample VIII (CC-Michel): a muon candidate in the TPC and a delayed Michel electron in FGD1 indicating the presence of a low momentum π + below tracking threshold.
The kinematic distributions of each sample separated by true topology are shown in Fig. 11.The data clearly shows a deficit compared to the nominal MC prediction for the CC-1π + sample while the opposite is seen in the CC-Other sample, highlighting the need to include the control samples for a data driven background constraint.This deficit of CC-1π + events has been observed in previous ND280 analyses [69,71].However in the CC-Michel sample, which contains mostly CC-1π + events, the data has a similar overall normalization compared to the nominal MC prediction, and also shows an excess at the peak of the distribution.This tension between the CC-1π + and CC-Michel samples and the impact on the analysis is discussed further in Sec.V.

INGRID
The INGRID selection includes a single control sample to select events with a single pion candidate track.Events must contain exactly two or three tracks that share a common vertex, with the highest-momentum minimum-ionizing track labelled as the muon candidate, the other minimum-ionizing track as the pion candidate, and a third track that if present must be proton-like.The PID cuts have been tuned for this sample to have a higher efficiency for selecting pion tracks compared to selecting proton tracks for the signal sample.The kinematic distribution of the control sample separated by true topology is shown in Fig. 12. Stopping and escaping events are considered together as a single CC-1π (events with a muon and a single charged pion track) sample in this analysis (Sample X).Similar to the ND280 control samples, the INGRID data shows a deficit of interactions producing a pion compared to the nominal MC prediction.

IV. ANALYSIS STRATEGY A. Binned likelihood fit
This analysis uses an unregularized binned maximum likelihood fit similar to the analyses in Refs.[67-69, 71, 72], to fit a set of signal and control samples to provide a data-driven background constraint, to unfold the detector effects, and to extract the number of selected signal events in the analysis bins.The ND280 and INGRID samples are fit simultaneously to extract the CC-0π cross section for each detector and produce a correlated result.For the purposes of the analysis, ND280 and INGRID events occupy different bins but are otherwise treated similarly.The analysis framework has been significantly improved compared to previous T2K CC-0π results (specifically Refs.[67][68][69]71]), for example including an improved treatment of the MC statistical uncertainty and principal component analysis to reduce the dimensionality of the fit.
This method varies the input MC using a set of fit parameters for both signal and background events to find the best agreement to the data, and the values and corresponding errors of these parameters at the best-fit point are then used for the cross-section extraction.The primary parameters of interest in the fit are the "template parameters" c i which scale the total number of signal events in each kinematic truth bin i (seventy in total for this analysis), and are completely free parameters with no prior constraint.The rest of the parameters are the systematic (or nuisance) parameters that describe variations to the flux, detector, and neutrino interaction model (described in Sec.IV B).Separate flux and detector parameters (including correlations when available) are included for ND280 and INGRID, while both detectors use the same neutrino interaction model parameters.
The best-fit parameters are found by minimizing the negative log-likelihood ratio (also approximated as the chi-square), and is split into a statistical and systematic contribution as follows: where and Equation 3 is the modified statistical log-likelihood ratio following the Barlow-Beeston method [73,74] for including the uncertainty of finite MC simulation.N MC j and N obs j are the number of simulated and observed events for each reconstructed bin j.The Barlow-Beeston scaling parameter for each bin β j is given by the following: (5) where σ 2 j is the relative variance of the number of MC events N MC j in the bin.In the limit of infinite MC simulation, σ 2 j → 0 and β j → 1 giving the standard Poisson log-likelihood ratio.Equation 4 is a Gaussian penalty term to account for the contribution from varying the systematic parameters ⃗ p during the fit compared to their fixed prior values ⃗ p prior and uncertainty.A covariance matrix V syst is used to describe the prior uncertainty and correlations between the parameters.
The input MC simulation for a reconstructed bin j is the sum of weighted signal and background events, and can be expressed as: where N sig ij and N bkg ij are the signal and background events for truth kinematic bin i and reconstructed bin j as predicted by the MC simulation, c i are the signal template parameters, and w ij are the weights as a function of the systematic parameters ⃗ p, and depend on the truth and reconstructed bins i and j.

B. Systematic uncertainties
There are three types of systematic uncertainties considered for this analysis and included in the fit as parameters: flux, detector, and neutrino interaction model uncertainties.
The neutrino flux uncertainty is parameterized as scale factors in forty total bins of true neutrino energy with separate flux parameters (or bins) for ND280 and IN-GRID.Only the ν µ flavor is considered for this analysis due to the small contribution of the νµ , ν e , and νe flavors.These parameters use the same energy binning scheme and can only affect events for their respective detector (the flux bin edges can be found in Appendix B).They have a prior constraint described by a covariance matrix, which includes the correlations between energy bins and between the fluxes at each detector.As shown in Fig. 13, the fluxes at ND280 and INGRID are highly correlated a priori.The high energy bin (10 to 30 GeV) for ND280 is less correlated than the lower energy bins due to an increase in the hadron uncertainties for that bin (mostly from kaon decay).Since the number of events in the analysis corresponding to this energy range is small, it has little effect on the analysis.For a given true neutrino energy bin, identical weights are given to signal and background events.The flux uncertainty is dominated by the hadronic multiplicity and decay modeling, along with other contributions, such as uncertainties in the horn current and alignment.The detector uncertainty is parameterized as scale factors on the event rate for each reconstructed bin with separate detector parameters for ND280 and INGRID.They have a prior constraint described by a covariance matrix including the correlations between the reconstructed bins for a given detector, however ND280 and INGRID use separate matrices and are completely uncorrelated in the fit.In principle, several detector uncertainties could be considered correlated between them, for example using the same pion secondary interaction modeling in geant4 for the detector simulation, but this was out of scope for this analysis.Independent and dedicated control samples for each detector are used to evaluate the detector uncertainties based on data-MC agreement.The largest contribution to the detector uncertainty is the pion secondary interaction modeling.Since the fit includes a detector parameter for each reconstructed bin, this adds up to many hundreds of parameters.Principal component analysis is used to reduce the total number of detector parameters by more than half by transforming the parameters to their eigenspace and removing the parameters that contribute less information according to their eigenvalues such that 99% of the total information in the covariance matrix is retained.
The neutrino interaction uncertainty is included in the fit using a set of twenty-one parameters that are designed to weight events based on aspects of the neutrino interaction model for both signal and background samples, including final state interactions (similar to Ref. [69]).A table listing the neutrino interaction parameter names and their priors can be found in Appendix C, and are described as follows.Variations to the signal model are included through changing the shape of the CCQE cross section by varying the axial mass M QE A , and with two parameters to alter the overall normalization of and shape of the 2p2h model.The resonant pion model has two shape parameters, the axial mass M RES A and the axial form factor at zero momentum transfer C A 5 , and a normalization for the non-resonant background I 1/2 .Additionally, two normalization parameters for CC-1π events are included to give the fit additional freedom to adjust the pion background and prevent over-fitting of the flux parameters.For deep inelastic scattering events, a custom shape parameter (DIS multi-pi shape) is used to give greater freedom at lower neutrino energy along with two normalization parameters for DIS and multi-pi events.The other major event topologies (coherent and neutral current) are each given a normalization parameter.Finally a set of six parameters are included to allow the pion FSI model to vary within the fit, separated by different reaction channels (absorption, production, charge exchange, and scattering) and pion momentum range.The prior uncertainty and correlations between parameters are encoded in a single covariance matrix.ND280 and INGRID share the same neutrino interaction parameters but the event weights are calculated for each detector separately.

C. Cross-section extraction
The flux-integrated differential cross section as a function of true muon kinematics x = p µ cos(θ µ ) for each detector is calculated using the following: where N sig i is the best-fit number of selected signal events in truth bin i summed across all samples, ϵ i is the bin-by-bin efficiency correction, Φ is the integral of the neutrino flux evaluated at the best-fit parameters, N nucleons is the number of target nucleons in the fiducial volume, and ∆x i is the bin width.The bin edges for the extracted cross section can be found in Appendix A.
The cross-section uncertainty is calculated by numerically propagating the post-fit uncertainty for the fit parameters assuming they follow a multivariate Gaussian distribution.The post-fit covariance matrix is Cholesky decomposed and used to create correlated random variations (or "throws") of the fit parameters that follow the same multivariate distribution as the covariance matrix.This procedure is repeated 10 4 times to sample the likelihood space encoded in the post-fit covariance matrix.For each thrown variation of the fit parameters, all the events are reweighted using the thrown parameter values and the cross section is recalculated as in Eq. 7. Additionally for each variation, the integrated flux is recalculated, the selection efficiency is allowed to vary based on the thrown parameters, and the number of targets is varied independently for each detector.The efficiency for each cross-section bin and its uncertainty for ND280 and INGRID are shown in Figs. 14 and 15 respectively.In general for both ND280 and INGRID the efficiency increases at forward angle and higher momentum as muons leave longer tracks.The number of target nucleons and its uncertainty for ND280 are 5.53 × 10 29 ± 0.67%, and INGRID are 1.76 × 10 29 ± 0.38%, and includes other elements in addition to carbon and hydrogen present in the fiducial volumes of each detector.The post-fit ν µ flux integral and uncertainty is 2.29 × 10 13 cm −2 and 6.0% for ND280 and 3.14 × 10 13 cm −2 and 6.1% for IN-GRID.The resulting distribution of cross-section values represent the plausible variations of the fit according to the post-fit uncertainties and correlations.Finally, the cross-section (dσ/dx) uncertainties are calculated using the ensemble of random throws and are parameterized using a covariance matrix, assuming the uncertainties are Gaussian distributed.
A set of fits were performed to estimate the total systematic uncertainty and the contribution from each systematic parameter class (flux, neutrino interaction, and detector) on the measured cross-section bins.A fit using only the template parameters is used as a baseline for the uncertainty (and corresponds approximately to the statistical uncertainty), and additional fits were performed that include each systematic parameter class to estimate its impact (in addition to the template parameters).Each fit used the nominal Monte Carlo prediction as the "data" so that the best-fit point is at the nominal value for every parameter; this guarantees each fit has the same best-fit point and allows for a more accurate comparison.The results are shown for the analysis crosssection bins for ND280 and INGRID in Figs.16 and 17 respectively, and show roughly equal contributions from the flux, neutrino interaction, and detector parameters to the total uncertainty.Additionally, the uncertainty for the low momentum bins in general is higher than mid to FIG. 14. Selection efficiency with post-fit uncertainty for the ND280 cross-section bins as function of true muon momentum in muon angle bins.Note that the final bin extending to 30 GeV/c has been omitted for clarity.
higher momentum bins for a given angle bin.This procedure cannot be applied in the same way to data as each fit will have a different best-fit point resulting in a different total uncertainty, preventing an equal comparison between the fits.

D. Validation
The cross-section extraction was validated using a series of mock data studies, where a known simulated data set was used as input to the fit and the performance of the extraction procedure was studied.These mock data sets cover a wide range of alterations, for example using data-driven modifications such as the deficit seen at low Q 2 by MINERνA [75], modifications to the resonant pion production model, or changes to the flux model.For each of these mock data sets, the cross-section extraction was able to recover the expected true cross section to within the 1σ uncertainties (often matching the true cross section nearly exactly), showing a robust procedure.The overall χ 2 agreement between the extracted and true cross section is also calculated as follows: where N is the number of cross-section bins, x i and x j are the ith and jth kinematic bin respectively, and V is the cross-section covariance matrix.
The post-fit p-value was calculated for the fit to data (or mock data) as another metric to check the validity of the result.First, numerous statistical and systematically varied samples of the MC prediction were produced and then fit to build a post-fit log-likelihood distribution from Eq. 2. The p-value is the fraction of the simulated likelihood distribution that is greater than the post-fit log-likelihood of the fit in question.A value of 0.05 for the p-value has been chosen as the threshold to require further investigation of the result.

V. RESULTS
The distributions of the reconstructed events used to calculate the cross section are shown in Figs.19 and 20 for ND280 and Fig. 21 for INGRID as a function of the muon kinematics compared to the nominal and post-fit MC predictions with data plus statistical errors overlaid.In general the fit is able to adjust the fit parameters as described in Sec.IV to match the observed data in the signal samples for both detectors within expected statistical fluctuations.
Nearly all systematic parameters had a post-fit value within their 1σ prior uncertainty, with the normalization parameters for DIS and multi-π events pulled further than 1σ to accommodate the difference between the nominal MC and the data in the control samples as discussed previously in Sec.III.Additionally, the 2p2h shape parameter was pulled to the boundary corresponding to pushing the distribution of 2p2h events toward smaller momentum and energy transfer.The p-value for the fit to data compared to the nominal input model was calculated to be approximately 0.01, and was investigated further to verify the robustness of the result.The extracted cross sections for ND280 and INGRID are shown in Figs.22 and 23 respectively, and includes an additional uncertainty to account for possible missing freedom in the fit indicated by the poor p-value.A discussion of the poor p-value is provided in the following section (Sec.V A).Additional plots of the extracted cross section including the last momentum bin extending to 30 GeV/c are shown in Appendix D.

A. Discussion on the poor p-value
As part of the data fit validation a p-value was calculated for the post-fit result as described in Sec.IV D to gauge the compatibility of the result with the nominal input model.The overall p-value for the original fit to data was calculated to be approximately 0.01, indicating that, relative to the input model, the observed data was an unlikely fluctuation.Based on the post-fit likelihood for each sample and the systematic contribution, the two main hypotheses were the pion modelling for background events (including the separation of samples for tracked versus Michel-tagged pions) and the high momentum bins considered in the fit.Several tests and different configurations of the fit were considered to see the impact on the result and to evaluate the need for an additional uncertainty.
To test the freedom of the pion model, two additional interaction model parameters were included to introduce more freedom at low pion and muon momentum where the majority of tension is present.One parameter was allowed to further alter the overall normalization of outof-fiducial volume events (which are more likely at lower muon momentum), and the other changes the kinematics of the intermediate ∆ decay in resonant pion interactions, which modifies the outgoing pion spectrum.In the fit to data, the rate of out-of-fiducial events was increased and the ∆-decay parameter was moved to increase the rate of low-momentum pions (which correlates with lowmomentum muons).However the total post-fit likelihood only improved by about 1% compared to the data fit without these extra parameters, showing little sensitivity to the change in pion kinematics.
Next an alternative configuration of the fit combines the ND280 CC-1π + and CC-Michel control samples (Samples VI and VIII respectively) in an attempt to quantify the role of the control samples on the background model, inspired by how these samples were used in previous T2K analyses [69,71].Overall the fit with the merged control samples gives similar results to the original fit, but with slightly better post-fit likelihood for the systematic contribution and a p-value of 0.024.The post-fit likelihood for the signal samples are nearly identical to the original showing little impact from the merged versus separated pion control samples.The merged sample performs slightly better compared to being separated, highlighting the tension between these samples as a driver of the poor p-value.
For the final test of the pion model, the fit was tested using only a single pion control sample (ND280 CC-1π + , ND280 CC-Michel, INGRID CC-1π) at a time.The post-fit likelihood for the other samples largely were the same and the systematic contribution to the likelihood was about 10% smaller (which is expected as fewer bins should in general easier to describe using the same systematic parameters) for each test case, indicating a robust result with respect to the pion control samples.The little to no effect on the signal samples across these tests is primarily due to the low contamination of pion events and the relative fraction of events in the signal versus control samples.
Another alternative configuration removes the high momentum bins extending to 30 GeV/c for each ND280 sample (except the backward bin) as these bins were contributing a disproportional amount to the total χ 2 relative to what would be expected from only statistical fluctuations.This was performed as an additional closure test to check that there was no pathological behavior regarding these bins that could be driving the p-value result.The fit performs moderately better as expected when removing problematic bins with a p-value of 0.065.Given the modest increase of the p-value and since the analysis is less sensitive to this kinematic region, these bins were considered not to be a problem.
The goal of these tests was to see if a large difference in p-value or post-fit likelihood was noticeable and if it required an additional uncertainty.The consistency of the post-fit results for the original fit to data and the variety of alternative fit configurations is evidence of a robust result for the analysis, and based on these additional investigations the small p-value is not an indication of a biased result and is driven by the measured data.Furthermore the slight variation in extracted cross section from the different tests were all well covered by the uncertainty on the original result.However, to be conservative, an additional uncorrelated uncertainty is added across all bins to cover the average difference between the original fit and the fit with merged ND280 control samples and the fit where the highest ND280 momentum bins were removed as these showed the largest difference in the extracted cross section.This additional uncertainty uses a similar method to the PDG to calculate how much the errors needed to be scaled to cover the different central values at 68% confidence level by construction [76].The percent error increase for each kinematic bin is shown in Fig. 18 where most bins had a 2% or less increase in error and the large increases generally correspond to the highest momentum bins.

B. Model comparisons
In this section, the measured cross sections are compared to a small selection of neutrino interaction models.The agreement between the measurement and models is quantified via the χ 2 (relative to the number of degrees of FIG.18. Percent error increase for each cross-section bin (flattened as a 1D array) from the the additional studies for the small p-value.freedom) as in Eq. 8 where σ true is replaced by the model prediction σ model .The model predictions were produced by generating a sufficiently large number of neutrino interactions on hydrocarbon using the T2K on-and offaxis flux for each model.Generated events that satisfy the CC-0π signal definition are selected to calculate the flux integrated cross section for each detector.The comparisons between this result and following models were facilitated by the nuisance software package [77].
(i) neut 5.5.0 [52,53] : several different NEUT configurations were used: two using the Benhar et al. [54] spectral function approach for the nuclear ground state with different values of the axial mass (M QE A = 1.03 and M QE A = 1.21 GeV/c 2 ), and one using the Nieves et al. model [20] for the nuclear ground state and quasi-elastic interactions (which is based on a local Fermi gas).All NEUT models use the same pion production model, multi-nucleon model from Nieves et al. [20], and FSI interaction model as described earlier in Sec.III A. This is a newer version of NEUT compared to the nominal MC production (NEUT 5.3.2) used in the analysis.
(ii) genie 3.0.6[78,79]: two different GENIE configurations were used: the "G18 01a" tune, which includes the Bodek-Ritchie modified relativistic Fermi gas (BRRFG) for the ground state nuclear model, the GENIE empirical multi-nucleon model [80], the Rein-Sehgal model for pion production, and the hA model for FSI; the "G18 10b" tune, which uses the local Fermi gas (LFG) for the nuclear model, the Figures 24 and 25 show the data compared to each generator's implementation of a LFG nuclear ground state plus the Nieves et al. multi-nucleon model.The pion production models will be roughly similar, however the FSI treatment is different between each prediction.The generators mostly differ in the normalization for the ND280 bins at the middle momentum range around the T2K flux peak energy of 0.6 GeV, however all show very similar INGRID predictions.For this particular set of models and generator versions, NEUT performs notably better than GENIE and NuWro.
Figures 26 and 27 show the data compared to several different multi-nucleon predictions using NuWro with a LFG ground state and the same parameters for all other aspects of the generation.The predictions are very similar between the different multi-nucleon models as implemented in NuWro, with a slight preference for SuSAv2.

C. Comparison to previous result
This analysis uses the same binning for the ND280 samples as the CC-0π analysis from Ref. [69], allowing for a direct comparison between the results.The main differences are the inclusion of more data for this result (T2K Run 8), increasing the neutrino-mode sample statistics by approximately a factor of two, and the configuration of the fits, where this analysis did a neutrinoonly joint fit of on-and off-axis data and Ref. [69] did a joint anti-neutrino and neutrino fit with only off-axis data.Both results are shown in Fig. 28 with the final bin extending to 30 GeV/c omitted for clarity.The majority of the bins agree within their 1σ error bars, and show a trend for this result to report a smaller cross section at medium to higher muon momentum (above 0.8 GeV/c) that is more pronounced at more forward-going angles.Additionally, the high fluctuation in the cross section seen in the 2.0 to 3.0 GeV/c momentum bin in the most forward angle bin (0.98 < cos θ µ < 1.00) is now smaller and closer in value to the neighboring bins compared to previous T2K CC-0π results.

VI. CONCLUSIONS
This paper presents the first measurement of neutrino interactions without pions in the final state using multiple energy spectra at T2K with the on-and off-axis near detectors.The analysis was performed using a joint maximum likelihood fit with signal and control samples from both detectors to minimize the background uncertainties and perform the unfolding from reconstructed to truth variables.The results include the cross-section measurement at each detector and the correlation between them, providing additional information compared to the individual measurements.Generator models continue to struggle to describe the data, and for the comparisons performed in this paper, the NEUT implementation of a LFG ground state plus the Nieves et al. multi-nucleon model has the smallest χ 2 /N ∼ 1.66, which is still very poor agreement.
This analysis is the next step in combined measurements at T2K and further opens up the possibility for more complex combinations of analyses.Only neutrinomode data was considered for this first analysis using multiple energy spectra, but future analyses will include the anti-neutrino data.Additionally, future versions of this analysis will include the T2K replica target measurements from NA61/SHINE [84] for the flux modeling, and updates of the neutrino interaction model.Since this analysis was finalized, the WAGASCI [85] and BabyMIND [86] detectors were added to the T2K near detector hall at an off-axis angle of 1.5 degrees and have started taking data.WAGASCI/BabyMIND data could be used to extend this analysis to use three different energy spectra, and provide additional statistics for interactions on both hydrocarbon and water.The upcoming J-PARC accelerator upgrade [87] will increase the beam power providing a higher rate of data taking.Finally, the ND280 upgrade [88] will increase the detector capabilities, providing increased angle coverage, better low momentum tracking, and additional target mass.
The data release for this analysis is hosted at Ref. [89].It contains the best-fit cross-section results, the nominal MC prediction, the associated covariance matrix, the analysis binning, and the flux histograms.

FIG. 1 .
FIG. 1. Multi-nucleon cross section on 12 C as a function of energy for the Nieves et al. and the Martini et al. models.
INGRID and ND280 are exposed to the same neutrino beam, but are placed at different angles relative to the beam center which gives a different integrated flux and energy spectrum for each detector.The neutrino flux peaks around 0.6 GeV at ND280 and around 1.1 GeV at INGRID, and the nominal ν µ fluxes are shown in Fig. 2. The beam composition at INGRID and ND280 when running in neutrino mode is shown in Tab.I.

FIG. 5 .
FIG. 5. A schematic view of the Proton Module and the standard modules.

FIG. 8 .
FIG.8.Event distribution for measured data and MC prediction in reconstructed muon momentum and angle for the ND280 signal samples stacked by true topology.The purity of each topology is listed in the legend, and the last bin for muon momentum contains all events with momentum greater than 5 GeV/c.

FIG. 11 .
FIG.11.Event distribution for measured data and MC prediction in reconstructed muon momentum and angle for the ND280 control samples stacked by true topology.The purity of each topology is listed in the legend, and the last bin for muon momentum contains all events with momentum greater than 5 GeV/c.

FIG. 12 .
FIG. 12. Event distribution for measured data and MC prediction in reconstructed equivalent distance in iron and angle for the INGRID control sample stacked by true topology.Through-going events are all placed in the final distance bin.The purity of each topology is listed in the legend.

1 FIG. 13 .
FIG.13.Input flux correlation matrix binned in neutrino energy for both ND280 and INGRID.The flux is highly correlated both across the energy spectrum and between the detectors.

FIG. 15 .
FIG. 15.Selection efficiency with post-fit uncertainty for the INGRID cross-section bins as function of true muon momentum in muon angle bins.Note that the final bin extending to 30 GeV/c has been omitted for clarity.

FIG. 19 .
FIG.19.Event distribution for measured data and the pre-/post-fit MC prediction in reconstructed muon momentum and angle for the ND280 signal samples.

FIG. 20 .
FIG. 20.Event distribution for measured data and the pre-/post-fit MC prediction in reconstructed muon momentum and angle for the ND280 control samples.

TABLE I .
Neutrino beam composition at INGRID and ND280.
[16]ear ground state with the Nieves multinucleon model, and multiple using a local Fermi gas (LFG) for the nuclear ground state with a different available multi-nucleon model: SuSAv2[83], Nieves et al., or Martini et al.[16].All configurations use the same models for pion production and FSI interactions and setM QE A = 1.03 GeV/c 2 .Tab.III contains the χ 2 for the joint result and the χ 2 considering a single detector.The joint χ 2 for a given model comparison will be slightly different from the sum of the individual χ 2 values due to the correlations between the detectors.Overall the generator predictions do not describe the data well according to χ 2 /N values ranging from approximately 1.5 to 3.0 for N = 70 degrees of freedom (measured bins) with 58 ND280 bins and 12 INGRID bins.The larger χ 2 values for ND280 compared to INGRID is primarily due to ND280 having a finer binning than INGRID.A unique aspect of this measurement is the ability to compare how a given model does for ND280 and INGRID individually and how the full result with correlations between ND280 and INGRID is better or worse than the naive sum.For example, the two GENIE models used in this paper show opposite behavior: one model describes ND280 better than the other but does worse describing INGRID and vice versa.

TABLE III .
Agreement between this result and the various model comparisons as measured by the χ 2 for both the joint result and when compared to each detector individually.ND280 has 58 cross-section bins and INGRID has 12 crosssection bins for a combined 70 total bins.