Simultaneous measurement of the muon neutrino charged-current cross section on oxygen and carbon without pions in the final state at T2K

This paper reports the first simultaneous measurement of the double differential muon neutrino charged-current cross section on oxygen and carbon without pions in the final state as a function of the outgoing muon kinematics, made at the ND280 off-axis near detector of the T2K experiment. The ratio of the oxygen and carbon cross sections is also provided to help validate various models' ability to extrapolate between carbon and oxygen nuclear targets, as is required in T2K oscillation analyses. The data are taken using a neutrino beam with an energy spectrum peaked at 0.6~GeV and comprises 57.34$\times$10$^{19}$ protons on target. The extracted measurement is compared with the prediction from different Monte Carlo neutrino-nucleus interaction event generators, showing particular model separation for very forward-going muons. Overall, of the models tested, the result is best described using Local Fermi Gas descriptions of the nuclear ground state with RPA suppression.


I. INTRODUCTION
The on-going long baseline (LBL) neutrino oscillation experiments, such as T2K and NOvA, are measuring the neutrino oscillation parameters with unprecedented precision and shedding light on the two known unknowns: neutrino Mass Hierarchy (MH) and Charge-Parity (CP) violation in the lepton sector [1][2][3][4][5]. A precise knowledge of neutrino interactions is a critical input for the study of neutrino oscillations not only for current LBL experiments but also for future experiments such as DUNE [6] and Hyper-Kamiokande [7]. Indeed, the precise determination of the MH and the measurement of the CPviolating phase in the PMNS mixing matrix [8,9] require the systematic error on predicted neutrino interaction event rates to be reduced to a few percent, of which the uncertainties related to neutrino interactions are currently the main contribution.
Although the presence of a near detector dramatically decreases uncertainties through constraints on the unoscillated neutrino flux, proper modelling of neutrino interactions is still critical for correct extrapolation of the expected event rate from the near to the far detector, which have different incoming neutrino energy spectra and may also have different acceptances and target materials. This is the case for T2K, where the near detector target regions are primarily composed of hydrocar- * 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.
bon, with only passive water sections, and have a limited acceptance to high-angle and backward-going particles, while the far detector, Super-Kamiokande [10], is a 4πacceptance Water Cherenkov detector. Beyond providing essential input for the prediction of the event rate at the far detector, the modelling of neutrino interactions is also important for estimating the bias and spread of any metric to determine the neutrino energy from its interaction products, which is a crucial input to neutrino oscillation analyses. Overall, the study of neutrino-nucleus interaction cross sections is imperative for current and future LBL experiments, with the study of interactions on water being particularly important for T2K. The neutrino-induced Charged Current Quasi Elastic (CCQE) interaction can be written as: where ν is the incoming neutrino, n and p represent the struck neutron and outgoing proton and is the charged lepton of the same flavour as the neutrino [11]. CCQE, also often referred to as '1p1h' (one-particle onehole), is the dominant reaction mode at T2K neutrino energies (peaked at 600 MeV) and therefore it is the interaction which is most important to characterise for T2K's neutrino oscillation measurements. While CCQE interactions with free nucleons are relatively simple to model [12], the situation becomes much more complex when the struck nucleon is bound inside a nucleus, that has an unknown initial momentum and binding energy. Moreover, the Final State Interactions (FSI) of outgoing hadrons inside the nuclear medium make CCQE interactions practically indistinguishable from meson-production interactions with subsequent mesonabsorption FSI. Interactions with multiple nucleons inside the nucleus can also leave a meson-less '2p2h' (two particle, two hole) [13] final state, which can also be confused with CCQE. Direct identification of solely CCQE interactions (or any specific interaction mode) is therefore difficult. In order to avoid highly model-dependent background subtractions, the experimental neutrino scattering community has developed the practice of publish-ing measurements of experimentally accessible final state topologies. In the case of T2K, the most relevant topology, accounting for the vast majority of events used by the far detector in oscillation analyses, are those with: one charged lepton; any number of nucleons; and nothing else (often called CC0π). Furthermore, the additional interaction modes and nuclear effects that contribute to a CC0π measurement are themselves important to understand for T2K neutrino oscillation measurements.
In this paper we present, for the first time, a combined measurement where the muon-neutrino-induced CC0π double differential cross sections on oxygen and carbon, as well as their ratio, are simultaneously extracted at the T2K off-axis near detector, ND280, as a function of the outgoing muon kinematics. This measurement provides a detailed probe of the CC0π interactions, which are of critical importance for T2K's oscillation analysis whilst also providing an additional exploration of nuclear-medium effects through the simultaneous analysis of neutrino interaction on two different nuclear targets. By measuring interactions on carbon and oxygen at the same time, and thereby providing a much improved understanding of how they may differ, this analysis complements other CC0π measurements on only carbon from T2K [14][15][16] in addition to those made by MINERvA [17][18][19][20][21] and Mini-BooNE [22] [23]. It also provides a validation and improvement on the first CC0π measurement on water for an incoming beam of muon (anti)neutrinos, published by T2K in Ref. ([24]) [25] using a different sub-detector at ND280 with different analysis techniques.
The paper is organised as follows: after a description of the T2K experiment in Sec. II, the data and Monte Carlo (MC) simulated data samples are outlined in Sec. III. The analysis strategy is then reported in Sec. IV A, including the description of the event selection, the cross section extraction procedure and the estimation of uncertainties. The paper ends with the presentation of the results, compared to a large number of models, in Sec. V, before conclusions are presented in Sec. VI.

II. THE T2K EXPERIMENT
The Tokai-to-Kamioka (T2K) experiment [26] is an accelerator-based long-baseline neutrino oscillation experiment located in Japan. Beams of predominantly muon neutrinos or anti-neutrinos are produced by directing a proton beam from the J-PARC accelerator complex in Tokai into a 90 cm long graphite target. The neutrinos then travel to the Super-Kamiokande far detector, 295 km from the neutrino production point [27]. The beam centre is directed 2.5 • away from the location of Super-Kamiokande, in order to achieve a narrowly distributed neutrino flux around the peak energy (∼ 600 MeV). The off-axis neutrino flux prediction, which will be discussed in more detail in Sec. III, is available in Ref. [28]. In order to characterise the unoscillated neutrino energy spectrum, to identify remaining intrin-sic backgrounds in the beam and to measure neutrino nucleus interactions, T2K also includes a near detector complex, located 280 m from the neutrino production point. It is the 2.5 • off-axis ND280 detector within this complex which is used for the analysis presented in this manuscript.
In this work, the two FGDs are used as the neutrino interaction targets whilst both the FGDs and TPCs are used as tracking detectors. The most upstream FGD (FGD1) primarily consists of polystyrene scintillator bars, with layers oriented alternately along the two detector coordinate axes transverse to the incoming neutrino beam, thus creating an 'XY module' and allowing 3D tracking of charged particles. The downstream FGD (FGD2) has a similar structure, but the polystyrene bars are interleaved with inactive water layers. The scintillator layers of both FGDs are made of 86.1% carbon, 7.4% hydrogen and 3.7% oxygen by mass, while the water modules are made of 73.7% oxygen, 15.0% carbon and 10.5% hydrogen; small fractions of Mg, Si and N are also present in both FGDs. A schematic of the two FGDs, as well as the chosen Fiducial Volume (FV) is shown in Fig. 2, illustrating that the FGD1 FV consists of 28 scintillator layers (i.e. 14 XY modules), while the FGD2 FV consists of 13 scintillator layers (i.e. 6 X modules and 7 Y modules) and 6 water modules. An XY module has a similar thickness to a water module. Overall, the considered total FV is made of ∼ 75% of hydrocarbon and ∼ 25% of water.

III. DATA AND MONTE CARLO SAMPLES
The analysis presented here uses T2K data spanning Runs 2 to 4, as reported in Tab. I, for a total of 57.34 × 10 19 Protons on Target (POT) taken with the beam mode producing predominantly muon-neutrinos (as opposed to anti-muon neutrinos).
The analysis of the neutrino data relies on the comparison of the measured quantities with simulation in order to correct for flux normalization, for detector effects and to estimate the systematic uncertainties.
The T2K flux simulation [27] is based on the modelling of interactions of protons with the fixed graphite target using the FLUKA 2011 package [34,35]. The modelling of hadron re-interactions and decays outside the target is performed using GEANT3 [36] and GCALOR [37] software packages. Multiplicities and differential cross sec-  tions of produced pions and kaons are tuned based on the NA61/SHINE hadron production data [38][39][40] and on data from other experiments [41][42][43], allowing the reduction of the overall flux normalisation uncertainty to 8.5%. The corresponding POT for simulated data is also reported in Tab. I. Neutrino interaction cross sections with nuclei in the detector and the kinematics of the outgoing particles are simulated by the neutrino event generator NEUT 5.3.2 [44,45]. The final state particles are then propagated through the detector material using Geant4 [46] before the readout is simulated with a custom electronics simulation.
NEUT version 5.3.2 describes CCQE neutrino-nucleon interactions according to the spectral function (SF) approach from Ref. [47] where the axial mass used for quasi-elastic processes (M QE A ) is set to 1.21 GeV, based on the Super-Kamiokande measurement of atmospheric neutrinos and the K2K measurement on the accelerator neutrino beam [48]. The resonant pion production process is described by the Rein-Sehgal model [49] with the axial mass M RES A set to 1.21 GeV. The modelling of 2p2h interactions is based on the model from Nieves et al. [50]. The deep inelastic scattering (DIS), relevant at neutrino energies above 1 GeV, is modeled using the parton distribution function GRV98 [51] with corrections by Bodek and Yang [52]. The FSI, describing the transport of the hadrons produced in the elementary neutrino interaction through the nucleus, are simulated using a semi-classical intranuclear cascade model [44,45].
As described in Sec. IV F and V, many other models and generators are considered for validations of the cross section analysis framework and the subsequent comparison with extracted results.

A. Goals and sample definition
The aim of this measurement is to extract the muon neutrino flux-integrated double-differential CC0π cross section simultaneously on oxygen and carbon nuclei as a function of the outgoing muon kinematics using the ND280 off-axis detector. For the first time the FGD1 and FGD2 detectors are used to simultaneously extract oxygen and carbon cross sections, thus accounting for correlations between them and also allowing a calculation of the cross section ratio. Since no single neutrino interaction target is completely dominated by oxygen, carbon interactions represent the main background for oxygen interactions. Both oxygen and carbon CC0π interactions are driven by the same physics and it would not be consistent to assume to know the latter to extract the former. A simultaneous measurement is therefore the best method to correctly disentangle the oxygen cross section from the carbon one in a Tracker based analysis.
In addition to using the two FGDs together to separate the two target nuclei, the reconstructed start point of the muon track in FGD2 is also employed to identify a sub-sample of events with a higher proportion of oxygen interactions. This technique is illustrated in Fig. 3, which demonstrates that interactions happening on water are mainly reconstructed in the X (Y) layers if the muon track is forward-(backward-) going. Overall, three categories of events are considered depending on the reconstructed starting position of the muon track: • samples with the muon track starting in FGD2X are oxygen-enhanced; • samples with the muon track starting in FGD1 and FGD2Y are carbon-enhanced.
This separation of carbon-and oxygen-enhanced event categories allows one to act as a control sample for the "background subtraction" of the other. Tab. II summarises the predicted sub-detector compositions for CC0π interactions.
A CC0π selection is applied in the FGD1 and FGD2 fiducial volumes and further split into FGD1, FGD2X and FGD2Y detector categories, depending on the starting position of the reconstructed muon track. In addition to the selection of CC0π events, this analysis also employs two control samples specifically designed to constrain and validate the modelling of the primary backgrounds to the main selection (these are also split FIG. 3. Schematic view of the FGD2 and of the technique employed to select oxygen-enhanced and carbon-enhanced samples based on the reconstructed muon track's start position. Yellow stars represent the true interaction position, while orange stars represent the reconstructed position. Interactions happening on water, are mainly reconstructed in the X (Y) layers if the muon track is forward-(backward-) going.
∼15% ∼60% into the three sub-detector categories). The details of the selection of signal and control samples are discussed in Sec. IV B. Following the identification of suitable signal and control samples, these are binned in terms of reconstructed muon kinematics and are used in a likelihood-fitter to subtract the background and unfold the detector response from the data (i.e. recover the number of selected signal events in 'true' muon kinematics). There is an unconstrained parameter controlling the scaling of the number of signal events in each bin of true muon kinematics for oxygen and carbon separately. Additionally, there are a variety of constrained (through a Gaussian penalty term) nuisance parameters allowing various background model variations and detector responses changes which are able to be constrained through dedicated control samples that are fit simultaneously with the signal samples. This fitting procedure is described in more detail in Sec. IV C. The results of the fit are then efficiency corrected and the flux and number for targets accounted for in order to extract the double differential cross section, as is detailed in Sec. IV D.
Systematic uncertainties are mainly evaluated by repeating the cross section extraction for a large ensemble of plausible variations to the input flux, detector and neutrino interaction models, whilst statistical uncertainties are evaluated using ensembles of data sets with Poissonian fluctuations of the number of real data events in each bin. This procedure, and the few exceptions to it, are discussed in Sec. IV E.

B. Event selections
The CC0π selection used in this analysis is the same as the one described for neutrino interactions in [16] and is summarised below. The selection achieves a wide acceptance in muon kinematic phase space by including high-angle and backward-going tracks in addition to the forward-going samples. As introduced in Sec IV A, this analysis uses FGD1 and FGD2 as a target for neutrino interactions whilst both the FGDs and the TPCs are used as tracking detectors. Additional information from the ECals and SMRD are also used in the case of characterising high-angle tracks.
After the first requirements on the data quality and the position of the vertex are fulfilled, the selection identifies interactions with only a single negatively charged minimally ionizing particle (the muon candidate) and any number of observed proton-like tracks (identified via the energy deposit of the track, its curvature in the TPC and/or its range in the FGD), which each must share a common vertex with the muon candidate. The particle type of each track is characterised by measuring its momentum (through its curvature if the track enters the TPC or its range if not) and energy loss. Interactions with an identified associated decay electron are also rejected, as these are likely to be from low momentum untracked pions decaying to muons and then to Michel electrons [53]. As introduced in Sec. IV A, each event is categorised based on whether it was observed to occur in FGD1, in an FGD2 X-layer or in an FGD2 Ylayer. For each sub-detector category (FGD1, FGD2X, FGD2Y), the selected events are then further divided into five exlusive signal samples depending on the detectors (FGD or TPC) used to measure the muon and proton (if there were any) kinematics and the observed proton multiplicity of the interaction (also shown in Fig. 4): sample I -µTPC: characterized by events with only one muon candidate in one of the TPCs; sample II -µTPC+pTPC: one muon and one proton candidate in one of the TPCs; sample III -µTPC+pFGD: one muon candidate in one of the TPCs and one or more proton candidates stopping in one of the FGDs; sample IV -µFGD+pTPC: one muon candidate tracked in one of the FGDs (and eventually the Ecal) and one or more proton candidates where one must enter one of the TPCs; sample V -µFGD: one muon candidate in one of the FGDs that reaches the ECal or SMRD and no identified proton candidate.
In Tab. III the number of selected events per signal sample and per sub-detector category is reported. Fig. 5 shows the event distribution per signal and control subsamples, compared with the T2K simulation predictions broken down per interaction and target nucleon type. The selection is highly dominated by events with one reconstructed muon and no other tracks. The predominance of CC0π interactions on carbon is evident in FGD1 and FGD2Y, while CC0π interactions on oxygen are dominant in FGD2X. It is also evident that the background comes principally from charged current events containing pions. These backgrounds primarily arise due to low momentum charged pions escaping identification. In order to constrain these backgrounds, two control samples are used in addition to the signal samples: sample VI -CC1π: characterized by events with one muon candidate and one π + candidate in the TPCs; sample VII -CC-others: one muon candidate + one π + candidate + an additional track in the TPCs; More details about the selection of these control samples can be found in [16]. In this analysis, the control samples are also divided into FGD1, FGD2X and FGD2Y categories, depending on the starting position of the muon track. The kinematics of the muon candidate in the first signal sample are shown in Fig. 6, where the predictions from the simulation are broken down by true interaction and target type. Similar plots for the other signal and control samples can be found in the supplementary material.
The ν µ CC0π cross section is extracted considering the contribution from all the samples, but it is important to keep the events with and without protons and with muon in different subdetectors separated in the analysis, as these are each affected by different systematic uncertainties, backgrounds and detector responses.
Following the selection, the events are binned according to the requirements of the cross section extraction. This involves ensuring the number of selected events in each bin is sufficient and that the binning is not finer than the detector resolution. For simplicity, the same binning is used for both carbon and oxygen cross sections and therefore the choice of the binning is driven by the oxygen events, since there are roughly three times more carbon events. The chosen binning is reported in Tab. IV. The corresponding efficiency for both oxygen and carbon events in the "truth" space (i.e. in the space    free from detector effects) is reported in Fig. 7. The slightly lower oxygen efficiency in the backward and high angle region is due to the difference between the FGD1 and FGD2 detector configurations, where in FGD2 there are the passive water layers interleaved with the active scintillator. The resultant loss in the efficiency mostly affects high-angle or backward tracks.

C. Fitting procedure
The analysis is performed using a binned likelihood fit with control samples to constrain the background, similarly to what is done in Ref. [14][15][16]54] in order to extract the selected number of signal events, unfolded from the detector response. This method is chosen as, in its unregularised form, it ensures no dependence on the signal model used in the simulation for the correction of detector smearing effects. Although model dependence can still enter through the efficiency correction, this is mitigated by choosing to extract a result as a function of observables which well characterise the detectors acceptance. Fitting-based unfolding methods, in contrast to commonly used iterative matrix-inversion methods (e.g. the commonly used method from [55] detector effects) space is allowed to float freely, whilst the background model predictions and the detector response are included as nuisance parameters with Gaussian penalty terms on the likelihood. In this analysis, a simultaneous fit is applied to all 21 of the signal and control sub-samples (s) described in Sec. IV B. For each of them, the predicted number of reconstructed events in the fit in the j th analysis bin, N j , can then be written as: Signal selection efficiency as a function of true muon kinematics using the binning adopted for the analysis (Tab. IV) for oxygen (black) and carbon (red) events. For readability purposes, the last momentum bin is cut at 5 GeV/c.
where i runs over the bins of the true muon kinematics, prior to detector smearing effects; are the numbers of signal (carbon and oxygen) and background events as predicted by the T2K Monte Carlo for the true bin i; w sig-C i , w sig-O i and w bkg i describe the alteration of the input simulation due to systematic parameters, described in Section IV E. The fit parameters of primary interest are the c i and o i : they are the factors that adjust the number of CC0π events on oxygen and carbon predicted by the MC to match the observed number of events in data. Finally, U ij is the detector smearing matrix that describes the probability to find an event of true bin i as reconstructed in bin j. This matrix is also altered by the detector systematic parameters, as described in Sec. IV E.
The best fit parameters are those that minimise the following likelihood: (2) or more explicitly: where N s j is the expected number of CC0π events in the sub-sample s and reconstructed bin j and N s,obs j is the observed number of events in each signal subsample s and reconstructed bin j. The second term (−2 ln(L syst )) is a Gaussian penalty term, where p are the nuisance parameters describing the effect of the systematics, 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. Finally, the two last terms (−2 ln(L reg p ) and −2 ln(L reg cosθ )) are additional and optional regularisation terms, similar to those used in Ref. [14,24].
Regularisation is the injection of some prior knowledge of the signal into the unfolding procedure in order to mitigate potential instability in the unfolded result, ensuring it is 'smooth' and physical. This can be required as, if the analysis binning is fine relative to the detector resolution, it is possible that many combinations of true bins lead to the same set of reconstructed bins [56]. With few exceptions, regularisation is routinely used in recent neutrino-nucleus cross-section measurements. In Eq. 3 a variant of Tikhonov regularisation is employed: the first regularisation term smooths the muon momentum bins within each cosθ µ bin, whilst the second one the cosθ µ bins, by using o k and c k , averaged values of the c i and o i over the considered angular bin. This is done separately for oxygen and carbon. Like all forms of regularisation, its presence introduces a bias in the extracted results, in this case to the shape of the input simulation, and potential underestimation of uncertainties. However, as detailed in Ref. [14,24], to reduce the risk of substantial bias towards the predicted shape, the 'L-curve' technique presented in Ref. [57] is used to choose the strength of the regularisation (p reg p and p reg θ ) directly from data. This technique is based on comparing the size of the regularisation term in the likelihood to the 'smoothness' obtained and balancing the two. This is discussed further in Appendix B. Particular care has also been given to verifying at each step of the analysis that the contribution from the two regularisation terms was minimal with respect to the dominant likelihood terms: −2 ln(L stat ) and −2 ln(L syst ). It was also always found that the regularisation on momentum bins accounts for a few percent of the total −2 ln(L), while the regularisation on angle bins accounts for some permille.
Despite the care taken to avoid bias, no regularisation method can be perfect and the application of any kind of regularisation will lead to at least some bias and underestimation of uncertainties, however small; therefore both regularised and unregularised results are reported. In general the regularised result is more stable with less strong off-diagonal covariances and so is better suited to 'by-eye' comparisons. Conversely, the unregularised result's large bin-to-bin variations and accompanying anticorrelations can cause misleading conclusions by-eye but is the result best suited to quantitative comparisons (e.g. the calculation of metrics for determining model agree-ment with the result). For this reason, χ 2 values from model comparisons are reported for both the regularized and unregularized results and show that any physical conclusions concerning data/model agreement are compatible with the two results, as is detailed in Sec. V and further discussed in Appendix B.

D. The extracted cross sections
The flux-integrated cross-sections and their ratio are evaluated in each bin i of muon momentum and angle (after the deconvolution of detector response): are the total number of signal events in bin i evaluated by the fit, O i and C i are the efficiencies, N FV O nucleons and N FV C nucleons are the number of nucleons in the fiducial volume, for oxygen and carbon respectively. Finally, Φ is the integrated flux for the T2K neutrino beam. In particular, the numbers of nucleons of the oxygen and carbon composing the fiducial volume of both FGD1 and FGD2 [58], have been estimated as:

E. Sources of uncertainties and their propagation
In order to produce meaningful results from the cross section extraction method presented in the previous sections, it is essential to evaluate and propagate potential sources of error. These include the statistical uncertainty on the data in addition to systematic uncertainties related to the modelling of the flux, of the detector response and of neutrino interaction cross sections.
Error Propagation. In order to propagate the impact of each systematic error source on the extracted cross section, elements of the cross section extraction procedure (the fit and the propagation to a cross section) are repeated for an ensemble of plausible variations ('toys') of the input MC. The way in which the ensembles of toys are built to characterise the uncertainty from each error source is detailed in the subsequent sub-sections. The sub-sections also detail the additional parameters that enter into the fits which, as discussed in Sec IV C, allow some of the sources of uncertainties to be constrained (mostly via the control samples). Statistical uncertainties are also calculated with toys in the same manner, but these are constructed by varying the number of entries in each reconstructed analysis bin according to a Poisson distribution centred around the number of events actually observed.
For the majority of the uncertainties, 1000 toys, in which each source of error is considered simultaneously, are used for propagation. For each toy a new cross section result is obtained following Eq. 4 and 5 where the impact of the uncertainties are included on all relevant parts of the cross-section extraction . The mean value of these results is taken as final cross section value and the spread is used to build a matrix of covariances to characterise the total uncertainty on the nominal extracted cross section (and, separately, on the extracted cross section ratio between oxygen and carbon). The covariances (V ij ) are constructed as: where t runs over the number of toys, the superscript FIT signifies the extracted results in the t th toy; dx j is the width of the j th bin in muon cosθ and momentum; and dσ FIT dx i,j are the mean differential cross-section values over 1000 toys in the j th or i th bin. V ij is therefore the total covariance matrix, including the statistical and systematic errors for the double differential cross sections. A 'shape only' matrix of covariances (W ij ) can also be calculated to be used to characterise the uncertainty on the result with normalisation information removed (this is useful for the model comparisons exhibited in Sec. V): where σ FIT,t indicates the integrated cross section over the full phase space as obtained in toy t. This method of error evaluation is used for all uncertainties other than those stemming from nucleon FSI and vertex migration, which are each discussed separately below. It should be noted that the method assumes that the distribution of toys within and between each extracted cross section bin is well approximated by a multi-variate Gaussian distribution. This was validated by analysing the ensembles of toys produced.
Flux uncertainty. The T2K flux prediction and uncertainties have previously been described in [27]. In each toy of the error propagation, the T2K flux covariance matrix is used to draw a random variation of the flux. The main impact of the flux is a larger overall normalisation uncertainty on the extracted cross section which enters through variations of the denominator in Eq. 4. The flux is not constrained in the cross section extraction procedure and so the resultant normalisation systematic uncertainty on the extracted cross section is, as in other T2K analyses (e.g. Ref. [15]), approximately 8.5%.
Detector response uncertainties. The detector response uncertainties considered are largely the same as described in Ref. [16] and are correlated between FGD1 and FGD2. The dominant systematics come from the uncertainties on the amount of background from the modelling of the pion secondary interactions and the TPC particle identification accuracy. To propagate the impact of the detector systematics, 500 toys of detector response variations are produced as variations to the input MC, considering the effect of all the detector systematics together. From this, a covariance matrix is built to characterise the uncertainties on the total number of reconstructed events in each bin of each sample used in the fit, for a total of 609 bins. This covariance matrix is then used to produce toys in the error propagation procedure described at the start of this section. Nuisance parameters are also added to the fit to constrain the impact of the detector uncertainties through the control samples. The number of nuisance parameters corresponds to the total number of reconstructed bins (609). Therefore, in order to reduce the number of fit parameters (which is essential for both the fit stability and to allow a reasonable computation time), a coarser reconstructed binning is used for these. Thus a second covariance matrix in this coarser binning is also produced to allow a calculation of the penalty arising from modifications of these parameters in Eq. 3.
Vertex migration uncertainty. Misreconstruction can lead to the reconstructed vertex position 'migrating' forwards when the first reconstructed hit is a layer downstream of the true one, or backwards when the reconstructed vertex is a layer upstream of the true vertex. The forward migrations come from a hit reconstruction inefficiency and constitute a small error that is treated as part of the other detector systematics. Backward migrations can come from low energy backward going particles whose energy deposits are mistakenly associated with the reconstructed muon track and therefore move the vertex one or more layers upstream. This latter uncertainty is particularly important to this analysis, in which samples in the FGD2 detector are divided depending on the position of the first reconstructed hit to attempt to isolate an oxygen enhanced sample of interactions, as described in Sec. IV B. The nominal simulations predict that about 14% of selected CC0π events in FGD2 are backward migrated and the uncertainty related to the estimation of this number has been evaluated in detail for this analysis.
In the case of a backward migrated event, the charge of the first hit (i.e. the starting point of the reconstructed track) is usually deposited by the stopping hadrons and not by the muon. Also, when a backward going hadron track is incorporated within the forward going muon track, the position of the first hits (hadron hits) are not expected to perfectly match with the rest of the track. Therefore the deviation (with respect to the rest of the track) and charge of the first few hits of a track can be used to estimate the backward migration rate. A fit of these variables, in which the backward migration rate was a free parameter, allowed a conservative estimation of the mismodelling of the backward migration rate to be around 30%. To estimate the impact of this uncertainty, an alternative input MC is produced where the reconstructed vertex of 30% of backward migrated tracks is artificially moved to the position of the true vertex (i.e. 30% of backward migrated events are moved to the category of non-migrated events). This alternative MC is used to fit the data and the difference between the cross section result obtained in this case and the nominal result is taken as uncertainty within each bin of muon kinematics. The backward migration uncertainty is considered as uncorrelated and so is added in quadrature to the diagonal elements of the covariance matrix as obtained in Eq. 6. As can be seen in Figs. 9 -11, the backward migration uncertainty affects mainly the oxygen cross section in the backward and high angle regions.
Number of target nucleons uncertainty. As discussed in Sec. IV D, the uncertainty on the number of nucleon targets for oxygen and carbon is 0.7% and 0.5% respectively. This uncertainty is propagated to the final results by varying the number of oxygen and carbon targets for each toy, taking into account the correlations. The uncertainty on the other materials, estimated to be at level of 10%, is also taken into account when producing the toys.
Modelling of signal and background interactions. The extraction of a cross section requires an estimation of the signal efficiency. Ideally the former should be a property of the detector but, without a very fine binning in as many observables that fully characterise the acceptance of the detector, there will always be some impact of the signal model on the detector efficiency. For example, the presence and multiplicity of additional nucleons can cause an event to be vetoed by the selection more or less often. In this analysis the signal is almost entirely made up of interactions from CCQE, 2p2h and resonant pion production with a subsequent pion absorp-tion FSI. The uncertainty on the neutrino-nucleon aspect of CCQE interactions is considered through variations of the nucleon axial mass M QE A (±0.41 GeV), that is fully correlated between oxygen and carbon. The uncertainty on the nuclear ground state model is controlled through variations of the Fermi motion and removal energy, very similarly to what is described in [59] but to be more conservative no correlations are assumed between oxygen and carbon nuclei. The uncertainty on 2p2h interactions includes a normalisation and a shape term. The former is taken to have a 100% uncertainty and the latter is treated as described in [2]. The 2p2h parameters are partially (30%) correlated between oxygen and carbon. Finally, pion absorption FSI and proton ejection FSI probabilities are also varied, details of the former can be found in [2] whilst the latter is described in more detail below. All the signal model variations are used, together with all the other systematics parameters, to create alternative input MC samples, but are not constrained in the fitter. It is clearly critical for a measurement's usefulness that the extracted cross section should not depend strongly on the modelling of it and indeed in this measurement these signal modelling uncertainties make up only a small portion of the overall error budget (and generally less than a 5% error) across almost all bins of the measured muon kinematics. The only exception is the backward going angular bin and the highest momentum bin of the high angle slice (0 < cos θ µ < 0.6), where the error can reach 10%. Beyond this, further tests to expose any significant model dependence in the cross section extraction are described in Sec. IV F.
The cross section extraction also relies on a prediction of the background event rate in each bin, which ideally should be well constrained by control samples. Although this analysis is high in signal purity (87% for FGD1 and 82% for FGD2), the backgrounds still require careful treatment. The dominant background is from resonant pion production in which neither the pion nor any associated Michel electron is observed directly. The variations of pion production processes are detailed in [59]. The same reference also details how pion FSI (in addition to the absorption process described above) are treated through parameters that alter different process interaction probabilities within the FSI cascade of the nominal MC.
These model uncertainties are propagated like the others, where many toys of plausible model variations are created by varying a set of underlying model parameters and modifying the input MC accordingly. Many of these parameters (all of those associated with the background processes other than pion FSI) are also allowed to float in the fit with a prior uncertainty (entering via the penalty term discussed in Sec. IV C and shown in Eq. 3) which is of the same size as the variation of the parameters used to build the toys.
Although the majority of the model uncertainties have been treated in similar ways for other T2K analyses, the analysis of nucleon FSI requires the same special treatment as detailed in [16]. Using current software tools, nucleon FSI cannot easily be varied using the input MC for the analysis, an uncertainty is built using two specially built samples of the NuWro event generator [60] (version 11q) with and without FSI. As discussed above, the primary way in which nucleon FSI enters into the uncertainty on the extracted cross section is through alterations to the efficiency. The difference in the efficiency of the two NuWro samples is therefore taken as a conservative additional uncertainty. As demonstrated in Figs. 9 and 10, this is generally small: less than 5% (and generally closer to 2%) for all bins other than for the very highest momentum bin at high or backward angles (where it can reach up to 15%).

F. Cross-section extraction validations
In order to validate the cross-section extraction procedure and diagnose any significant model dependence within it, a large number of 'mock data' studies were performed. It was validated that the procedure was able to accurately extract the true signal cross section from alternative simulations that were treated as data. These mock data sets include two different neutrino interaction generators as input data (GENIE 2.8.0 and NuWro 11q) in addition to large ad-hoc modifications of the input MC to simulate extreme variations in the signal. Importantly the modifications was calculated and applied in much finer binning than the analysis bins of Tab. IV in order to allow an alteration of events within each bin and therefore a representative variation of the signal efficiency. For each mock data sample, the regularisation strength, p reg p and p reg θ , was also re-estimated and their values were found to be fairly stable, p reg p being 4 or 5 and p reg θ between 4 and 7. The cross-section extraction was validated using a χ 2 test performed as: where σ meas and σ truth are the extracted and true cross sections (i.e. the cross section predicted by the MC acting as mock data) respectively. The values of the χ 2 were found, in all cases, to be lower than the number of analysis bins, indicating compatibility between the extracted cross section and the truth. The χ 2 were also calculated for different numbers of toys used in the uncertainty propagation method to calculate the covariance matrix (V ij ) in order to find the number of toys required to achieve a good statistical precision of the matrix elements (this was found to be 800 toys). Importantly, for each mock data set, the χ 2 was found to be very similar for regularised and unregularised results, showing that very little bias is introduced when the regularisation is applied for each of these mock data sets. The impact of regularisation was also evaluated on the real data and is discussed in Sec. V and Appendix B.

V. RESULTS AND COMPARISONS WITH MODELS
The event selection and cross section extraction procedure detailed in Sec. IV A is applied to the data samples introduced in Sec. III. Using the L-curve method discussed in Sec. IV C regularisation strengths are chosen as p reg p = 4 and p reg θ = 7 (see Appendix B for more details), similar to what was found in the mock data studies detailed in Sec. IV F. In this section the regularised results are shown, but for completeness the unregularised results are available as supplementary Material and a comparison between regularised and unregularised results is presented in Appendix B. As is detailed below, the use of regularisation has very little impact on model discrimination (as is shown in Tab. V).
The uncertainties on the extracted result and on the corresponding covariance matrix are calculated as detailed in Sec. IV E. 1000 toy fits were performed on the data, a number that was found to be sufficient to accurately calculate covariances. In Fig. 8, the distribution of the reconstructed events in the analysis binning for all the signal samples summed together is shown, as well as the comparison with the nominal MC and the mean of the fitted MC (over the many toys). Overall the fit is able to well reproduce the observed distributions. Similar plots for the control samples are available in Appendix A, showing these to also be accurately reproduced by the fit.
The final errors in each bin of the extracted cross section and cross section ratio are summarised in Figs. 9 -11, showing an approximate breakdown by error source. This breakdown is made by first running 1000 toys from only statistical fluctuations of the data before adding the systematic fluctuations and then each of the additional uncertainties described in Sec. IV E (the vertex migration and nucleon FSI). As expected, the statistical uncertainty on the oxygen cross section is higher than the one for carbon, since the number of oxygen events is roughly 1/3 of the number of carbon events. It can also be seen that the systematic uncertainties affecting the O/C ratio are reduced, since many of them (e.g. flux systematics) are fully correlated between oxygen and carbon. However, the ratio suffers from a higher statistical uncertainty, due to the intrinsic anti-correlation existing between the oxygen and carbon template parameters in each bin. The final correlation and covariance matrices (as calculated using Eq. 6) are shown in Fig. 12. From the correlation matrix it can be seen that the analysis binning choice, relative to the available statistics, and the application of a data-driven regularisation had mitigated the impact of anti-correlations between adjacent bins in the unfolding. However, it can also be seen that important correlations still remain, especially in the less statistically limited carbon cross section, demonstrating the importance of quantitative comparisons of the data to models which consider all elements of the data covariance (such as the χ 2 comparison shown in Eq. 9) The extracted double differential cross sections per nucleon are shown for oxygen and carbon together in Fig. 13. Generally, a slightly higher oxygen cross section is observed, with the exception of the most forward going angular bin, in which the carbon cross section is a little larger.

A. Comparisons to models
In the following, the measured cross sections, and their ratio, are compared to different neutrino-interaction models and the level of agreement is quantified by the χ 2 statistics, as follows: It should be noted that, apart from when considering the ratio measurement, the overall normalization uncertainty (fully correlated between bins) constitutes a relatively large fraction of the uncertainty, between 20% and 60% depending on the bin. Therefore the χ 2 statistics may suffer from 'Peelle's Pertinent Puzzle' (PPP) [61,62], which describes how the implicit assumption in Eq. 9 that the variance is distributed as a multivariate Gaussian may not be well suited to highly correlated results. Therefore, to mitigate this problem the shape only χ 2 is also provided in Tab. V. This is estimated as follows: , (10) where σ model int.

int.
are the total integrated cross sections per nucleon estimated from the model and from the data, respectively.
The comparison of the measurements presented in this paper to the various models is performed in the framework of NUISANCE [63]. A sufficiently large number of events are generated on carbon and oxygen from each model using the T2K flux. From each model the events corresponding to this analysis' signal definition (CC0π) are then selected and used to calculate a cross-section per target nucleon.
In all the LFG models other than the one used by GiBUU, the Random Phase Approximations (RPA) corrections are applied, as computed in Ref. [79].
In Figs. 14-16, the result is compared to generators using differing models for the CCQE contribution and for the corresponding nuclear ground state: LFG (NEUT), SF (NuWro), SuSAv2 (GENIE) and GiBUU, while in Figs. 17-19, data is compared to: NEUT with SF, NuWro with LFG, GENIE with LFG and RMF(1p1h) + SuSAv2(2p2h). Finally Fig. 20 shows the breakdown by neutrino true interactions contributing to the CC0π channel for the NEUT 5.4.1 predictions.
The values shown in brackets in the legend of each figure represent the χ 2 as obtained from Eq. 9 for the entire measurement (oxygen and carbon, 58 bins). The χ 2 (full and shape-only) for all models are summarised in Tab. V. The oxygen-only and carbon-only χ 2 are also reported in the same table. These χ 2 have been obtained considering only the 29 oxygen or carbon bins and neglecting the correlations between the two measurements; although they thus neglect some information with respect to the full results, it remains interesting to consider them to quantify model agreement with each individual target. In addition to the χ 2 and χ 2 shape metrics, a partial χ 2 excluding the last cosθ µ bins is also shown in order to isolate the impact of this very forward bin where models seem to struggle the most (as is evident from Fig. 14). As can be  11. Summary of the uncertainties for the oxygen over carbon cross section ratio as obtained over 1000 toys with regularisation. The statistical error is in blue. Systematic errors are then sequentially added in quadrature starting with all of those addressed via the prior variation propagation method (light blue), followed backward migration (green). As described in the text, proton FSI errors are considered to be fully correlated between oxygen and carbon and thus canceled out in the ratio. seen from the table, the last cosθ µ bin is often responsible for a large portion of the χ 2 .
Finally, in Tab. VI, the values of the integrated cross sections per nucleon for carbon and oxygen are reported alongside the ratio for the integrated regularised and unregularised results. This is then compared with the expectations from all the tested models.

B. Discussion
Overall, from the models shown, the Valencia (LFG) model predictions for 1p1h and 2p2h (i.e. NEUT 5.4.1 LFG, Genie 3 LFG hN and Genie 3 LFG hA) show the lowest χ 2 in comparison with our data. This is evident from the full and shape-only χ 2 and also from the comparison plots themselves, indicating a genuine agreement considering all correlations and accounting for possible misleading full χ 2 from PPP. The agreement between the GENIE and NEUT implementations of the model is not surprising since where they differ is predominantly in the extrapolation of the Valencia inclusive model to exclusive predictions, which has a small impact when measuring only muon kinematics. It is also important to note that a large portion of the disagreement of other models stems from the most forward bin, where the role of RPA suppression is most important (without it the agreement would be very poor here); anyway, a slightly lower χ 2 for the Valencia model remains when considering only more intermediate kinematics, particularly when considering the shape-only χ 2 (as indicated in Tab. V). More generally, it can also be seen from the plots that, without the most forward angular bin, models that use dramatically different nuclear physics assumptions give similar predictions, all of which are generally in agreement with the result. This is mainly because model differences in this region of lepton kinematics are largely just normalisation changes which are not easily resolvable within current flux uncertainties. Separating these models is more possible by additionally measuring hadron kinematics (for example as T2K has measured in [14]), although when this is done none of these models is capable of describing all the data. It is interesting to note that the GiBUU prediction (also based on a LFG nuclear model) shows a large χ 2 in the most forward bin, where RPA effects are most important. GiBUU does not include an RPA suppression as it is suggested that its more sophisticated nuclear ground state description accounts for a large por- tion of the role of RPA [80]. GiBUU's transport approach to modelling FSIs (a more complete approach than commonly used cascades) also predicts a significantly larger pion absorption in the forward region [81], which could also contribute to the over prediction. Similarly to GiBUU, the SF models also show that a more sophisticated model of the nuclear ground state does not mean better agreement with the data. The SF predictions in NuWro and NEUT are very similar and, like GiBUU, struggle to describe the most forward bin. It can be seen from Fig. 20 that this is the region where 2p2h contributes most strongly and it may be that the addition of the Valencia 2p2h (based on a Fermi gas model) is too strong when applied on top of a SF prediction.
It can be seen that the SuSAv2 model (as implemented in GENIE) is also unable to describe the most forward bin, but this should not be surprising. SuSAv2 is based on extracting scaling functions from RMF and assuming super-scaling, however it is well known that at low momentum transfer (likely to be at forward angles) this is not so well satisfied [68]. As can be seen from Figs. 17-19, RMF is much more able to describe the forward bin for carbon (although struggles for oxygen).
Considering again Tab. V, it is clear that, in general, the χ 2 shape values show the same trend as the total χ 2 . The oxygen-only and carbon-only χ 2 show, in general, that all generators tend to slightly better agree with the carbon measurement than with oxygen measurement, other than the RMF(1p1h)+SuSAv2(2p2h) model that seems to slightly better reproduce the oxygen cross sections. Concerning the ratio, it can clearly be seen that model predictions of the differences between carbon and oxygen are so small that the data has very little power to offer any particular conclusion other than all tested models can describe the ratio reasonably well. Since the uncertainties in the ratio measurement are dominated by statistics of the data samples, more data in future T2K analyses will allow a greater precision. The integrated result on carbon and oxygen can also be considered, which has a much smaller statistical uncertainty and shows that all the generators predict a lower integrated cross section for both oxygen and carbon with respect to what is measured. Whilst the carbon disagreement is usually within one standard deviation, this is not true for oxygen.

VI. CONCLUSIONS
In this paper, carbon and oxygen CC0π muon neutrino double differential cross section measurements, as well as their ratio, as a function of muon kinematics has been presented as obtained from the ND280 tracker. The anal-ysis is performed with a joint fit on carbon-and oxygenenhanced selected samples of events, thus allowing a simultaneous extraction of the oxygen and carbon cross sections with proper correlations. The measurements have been done with and without the use of a data-driven Tikhonov regularisation; comparisons of the results show excellent compatibility and therefore demonstrate the ab-sence of significant model bias in the unfolding of detector smearing effects from the data. Overall, the accurate and simultaneous extraction of differential cross sections on the two nuclear targets that T2K utilises in its oscillation analysis allows to characterise the aspects of neutrino interaction cross sections most pertinent to this experiment. More generally the results also allow to probe nuclear medium effects that are expected to be the cause of dominant uncertainties in next generation LBL neutrino oscillation experiments. In this way it is hoped that measurements presented will be used to assist in the validation of input models to oscillation analyses whilst also providing new data for theorists and model builders to improve or tune their predictions.
An extensive comparison of the extracted results to some of the most commonly used and sophisticated neutrino interaction models available today shows a preference for CCQE models based on a relatively simple Local Fermi-gas nuclear ground state, as opposed to more involved spectral function or mean-field predictions. With current statistical uncertainties, the strength of this preference is currently dominated by the most forward angular slice where the nuclear physics governing low energy and momentum transfer interactions becomes most important. This is also where relatively poorly understood 2p2h and FSI effects are largest relative to the CCQE prediction. It therefore remains possible that the more sophisticated CCQE models are correct but are undermined by the more simple FSI models or the 2p2h predictions based on a Fermi-gas ground state that currently need to be added on top. Outside of this forward slice all tested models give predictions compatible with the results, despite containing very different nuclear physics making further model discrimination difficult.
Future analyses will aim to improve model separation through both the simultaneous measurement of hadron and lepton kinematics in addition to combing the current joint analysis on oxygen and carbon with the analysis on neutrinos and anti-neutrinos recently published in [16] whilst benefiting from improved constraints on the flux model.
The data release for the results presented in this analysis is posted at the link in Ref. [82]. It contains the analysis binning, the oxygen and carbon ν µ double-differential cross sections central values, their ratio and associated covariance and correlation matrices. Fig. 24 shows a comparison of the extracted correlation matrices for regularised and unregularised results, whilst Fig. 25 shows a comparison of the two extracted results. As expected, errors bin per bin are in general larger than for regularised results with stronger anti-correlation between nearby bins. It should be noted that the only reason why adjacent bins are more strongly correlated than others is due to regularisation. It can also be seen that there is more bin-to-bin variation in the unregularised result, particularly for the lower statistics oxygen measurement, which is compensated by the larger anticorrelations between them. Despite these differences, the regularised and unregularised results remain absolutely compatible and, as can be seen in Tab. V, the χ 2 values obtained from model comparisons are very similar for the two results. Critically, physics conclusions drawn from the regularised and unregularised results are the same.
We strongly encourage future users of the regularised results to validate any quantitative statistic from a model comparison against what is obtained from the unregularised results.