Measurement of neutrino and antineutrino oscillations by the T2K experiment including a new additional sample of ν e interactions at the far detector

The T2K experiment reports an updated analysis of neutrino and antineutrino oscillations in appearance and disappearance channels. A sample of electron neutrino candidates at Super-Kamiokande in which a pion decay has been tagged is added to the four single-ring samples used in previous T2K oscillation analyses. Through combined analyses of these five samples, simultaneous measurements of four oscillation parameters, j Δ m 232 j , sin 2 θ 23 , sin 2 θ 13 , and δ CP and of the mass ordering are made. A set of studies of simulated data indicates that the sensitivity to the oscillation parameters is not limited by neutrino interaction model uncertainty. Multiple oscillation analyses are performed, and frequentist and Bayesian intervals are presented for combinations of the oscillation parameters with and without the inclusion of reactor constraints on sin 2 θ 13 . When combined with reactor measurements, the hypothesis of CP conservation ( δ CP ¼ 0 or π ) is excluded at 90% confidence level. The 90% confidence region for δ CP is ½ − 2 . 95 ; − 0 . 44 (cid:2) ( ½ − 1 . 47 ; − 1 . 27 (cid:2) ) for normal (inverted) ordering. The central values and 68% confidence intervals for the other oscillation parameters for normal (inverted) ordering are Δ m 232 ¼ 2 . 54 (cid:3) 0 . 08 ð 2 . 51 (cid:3) 0 . 08 Þ × 10 − 3 eV 2 =c 4 and sin 2 θ 23 ¼ 0 . 55 þ 0 . 05 − 0 . 09 ( 0 . 55 þ 0 . 05 − 0 . 08 ), compatible with maximal mixing. In the Bayesian analysis, the data weakly prefer normal ordering (Bayes factor 3.7) and the upper octant for sin 2 θ 23 (Bayes factor 2.4). DOI: 10.1103/PhysRevD.96.092006


I. INTRODUCTION
Neutrino oscillations have been firmly established by multiple experiments.Super-Kamiokande (SK) observed an energy and pathlength dependent deficit in the atmospheric muon neutrino flux [1]; and Sudbury Neutrino Observatory (SNO) resolved the long-standing solar neutrino problem by demonstrating that the previously observed deficit of electron neutrinos from the Sun was due to flavor transitions [2].These two experiments, together with accelerator-based (K2K [3], MINOS [4]) and reactor-based (KamLAND [5]) long-baseline experiments measured the two mass-squared differences between mass eigenstates and two of the three mixing angles in the PMNS matrix.
The mixing angle, θ 13 , has been measured as nonzero by T2K [6,7], by reactor experiments [8][9][10], and more recently by NOνA [11].Establishing that all three mixing angles are nonzero opens a way to study CP violation in the leptonic sector through neutrino oscillations.CP violation in neutrino oscillations arises from δ CP , an irreducible CP-odd phase in the PMNS matrix.This phase introduces a difference in the appearance probability between neutrinos and antineutrinos.To investigate this phenomenon, after taking data with a beam predominantly composed of muon neutrinos in order to observe the appearance of electron neutrinos at the far detector, T2K has switched to taking data with a beam predominantly composed of muon antineutrinos.A direct measurement of CP violation can then be obtained by comparing ν µ → ν e and ν µ → ν e channels.
To produce neutrinos, protons extracted from the J-PARC main ring strike a target producing hadrons which are then focused and selected by charge with a system of magnetic horns.The hadrons decay in flight, producing an intense neutrino beam.A beam predominantly composed of neutrinos or antineutrinos can be produced by choosing the direction of the current in the magnetic horn.T2K uses the so-called off-axis technique with the beam axis directed 2.5 • away from SK in order to produce a narrow-band neutrino beam, peaked at an energy of 600 MeV, where the effect of neutrino oscillations is maximum for a baseline of 295 km.Neutrinos are also observed at a near detector complex, installed 280 m from the target, comprising an on-axis detector (INGRID) which provides day-to-day monitoring of the beam profile and direction, and a magnetized off-axis detector (ND280), at the same off-axis angle as SK, which measures neutrino interaction rates before oscillation.
The analyses described in this paper are based on an exposure of 7.482×10 20 protons on target (POT) in neutrino mode (ν-mode) and 7.471×10 20 POT in antineutrino mode (ν-mode) collected at SK during seven physics runs as detailed in Tab.I.The neutrino oscillation parameters are measured by combining ν µ and ν µ disappearance channels with ν e and ν e appearance channels, using the same analysis techniques described in Ref. [12].The analyzed data set is the same as in Ref. [12], but an additional SK sample is included in the oscillation analysis.Previously, for the appearance channel, only the SK single-ring e-like interactions without additional activity in the detector were used for the oscillation analysis.The analysis presented in this paper includes an additional sample enriched in ν e interactions in which the e-like ring is accompanied by delayed Michel electron due to the decay chain π + → µ + → e + of π + 's produced in the neutrino interactions.This sample is currently only used in ν-mode, and increases the statistics of the ν e sample in SK by roughly 10%.The paper is organized as follows: the neutrino beam and the modeling of the neutrino fluxes are described in Sec.II.The neutrino interaction model developed for this analysis will then be described in Sec.III, followed by the selection of neutrinos in the near detector complex in Sec.IV.The neutrino flux and neutrino interaction inputs, and near detector selections are combined to reduce flux and cross-section uncertainties at the far detector as will be shown in Sec.V.The far detector selections are described in Sec.VI.The neutrino oscillations and the T2K oscillation analyses frameworks are then described in Sec.VII and in Sec.VIII respectively.Sec.IX is dedicated to a description of the impact of the uncertainties of the neutrino interaction model on the T2K oscillation analyses.Finally, the results of the oscillation analyses are presented in Sec.X and in Sec.XI and some concluding remarks are given in Sec.XII.

II. THE T2K BEAM
The neutrino beam is produced by the interaction of 30 GeV protons from the J-PARC main ring accelerator on a 1.9 interaction-length graphite target.Secondary hadrons, mainly pions and kaons, leaving the target pass through three electromagnetic horns [13], which are operated at a current of either +250 kA or −250 kA to focus positively or negatively charged particles respectively.The outgoing hadrons decay in a 96-m long decay volume, where a relatively pure beam of muon neutrinos is produced by the decay of positively charged hadrons in positive focusing mode (ν-mode), and a beam mostly composed of muon antineutrinos is produced in negative focusing mode (ν-mode).Protons and undecayed hadrons are stopped in a beam dump, while muons above 5 GeV pass through and are detected in a muon monitor (MUMON [14]), and are used to monitor the secondary beam stability.The T2K beamline hardware has been described in detail elsewhere [15].
The T2K neutrino flux at the near and far detectors in case of no neutrino oscillation is predicted by a simulation which has been described in detail in Ref. [16].Interactions of the primary proton beam, whose profile is measured for each run period by a suite of proton beam monitors, as well as subsequently-produced pions and kaons, are simulated within the graphite target by the FLUKA 2011 package [17,18].The predicted hadron production rates inside and outside the target are then adjusted based on the results from the latest analysis of the full 2009 thin-target dataset by the NA61/SHINE experiment [19][20][21], as well as other hadron production experiments [22][23][24].Particles which exit the target and subsequently decay are tracked through the horns and decay volume by a GEANT3 [25] simulation using the GCALOR [26] package.The predicted (anti)neutrino fluxes at the far detector for T2K Run 1-7 is shown for both ν-and ν-modes in Fig. 1.
Most of the "right-sign" ν µ flux (i.e.ν's in ν-mode and ν's in ν-mode) comes from mesons produced inside the target, dominantly "right-sign" (focused) pion and kaon production, which subsequently decay to produce muons.Interactions producing "right-sign" ν e 's also predominantly come from interactions in the target, with a larger fraction of ν e produced by kaon (rather than pion) decays.Interactions producing the "wrong-sign" ν µ flux have a higher fractional rate of out-of-target interactions, which are dominated by protons, neutrons, and pions scattering in the horns and decay volume walls.Interactions producing the "wrong-sign" ν e flux have a significant fraction of K 0 production from proton or neutron interactions, as well as charged kaon production.
In general, the ν-and ν-mode fluxes are similar at low energy, although the "right-sign" ν µ (and ν e ) flux in ν-mode is ∼15% higher around the flux peak than the "right-sign" νµ (and νe ) flux in ν-mode.The "wrongsign" background ν flux is also lower in ν-mode compared to the ν flux in ν-mode, especially at high energy.These differences are due to the higher production multiplicities of positive, rather than negative, parent particles.Uncertainties in the neutrino flux prediction arise from the hadron production model, proton beam profile, offaxis angle, horn current, horn alignment, and other factors.For each source of error, the underlying parameters in the model are varied to evaluate the effect on the flux prediction in bins of neutrino energy for each neutrino flavor as described in detail elsewhere [16].The uncertainties on the unoscillated ν µ and νµ beam fluxes at the far detector are shown in Fig. 2 and are currently dominated by uncertainties on hadron production.The uncertainties on the background ν e and νe fluxes from the beam are 7-10% in the relevant region.

III. NEUTRINO INTERACTION MODEL
The neutrino interaction model used in this analysis is based on NEUT [28] version 5.3.2, which includes many significant improvements over the old version 5.1.4.2, used in previous T2K oscillation analyses (described in detail in Ref. [27]).This model is constrained where possible by external experiments that are used to provide initial cross section parameter uncertainties.Such uncertainties are reduced using ND280 data, as explained in Sec.V. Alternative models are used to build simulated datasets to test the robustness of the T2K analysis against model-dependent assumptions, as explained in Sec.IX.This section describes the updated NEUT interaction model and alternative models used for the oscillation analyses.

A. Neutrino interaction model used in the oscillation analyses
The interaction rate at T2K energies is dominated by charged current quasi-elastic (CCQE) events, ν l n → l − p (ν l p → l + n).Because CCQE is a two-body process and the neutrino direction is known, the neutrino energy can be reconstructed from the outgoing lepton kinematics alone.However, nuclear effects and other processes, which have the same experimental signature of a single muon and no final state pions (CC0π, or CCQE-like), are indistinguishable from CCQE and can affect the reconstructed neutrino energy and thus the oscillation result if not accounted for [29][30][31][32][33][34].The T2K cross-section modeling has been updated to include recent theoretical models of these processes (full details can be found in Ref. [35]).In previous analyses, the CCQE model was based on the Llewellyn-Smith neutrino-nucleon scattering model [36] with a dipole axial form factor and BBBA05 vector form FIG. 2. The T2K fractional systematic uncertainties on the SK flux arising from the beamline configuration and hadron production prior to constraints from near detector data.Uncertainties are given for ν's in a ν-mode beam (top left), ν's in a ν-mode beam (top right), ν's in an ν-mode beam (bottom left), and ν's in an ν-mode beam (bottom right).For the ν-mode plots, the total current uncertainties (NA61 2009 data) are compared to the total uncertainties estimated for the previous T2K results (NA61 2007 data) [27].
factors [37], and used the Smith-Moniz Relativistic Fermi Gas (RFG) model [38] to account for the fact that the nucleons are bound in a nucleus.The main improvements available in NEUT v5.3.2 are the inclusion of the Spectral Function (SF) model from Ref. [39], which provides a more sophisticated description of the initial state of the nucleus than the RFG: the inclusion of the multinucleon interaction (2p2h) model from Refs.[40,41]; and the implementation of the Random Phase Approximation (RPA) correction from Ref. [40].The 2p2h model includes interactions with more than one nucleon bound within the nucleus, which contribute considerable strength to the CCQE-like cross section, and add significant smearing to the reconstructed neutrino energy distribution (as it is not a two-body process).RPA is a nuclear screening effect due to long-range nucleon-nucleon correlations which modifies the interaction strength as a function of four-momentum transfer, Q 2 .The models make different physical assumptions, so they cannot be combined arbitrarily.Two candidate models were considered for the default model: SF and RFG+RPA+2p2h.RFG+RPA+2p2h was selected as the default because it was most consistently able to describe the available Mini-BooNE [42,43] and MINERvA [44,45] CCQE-like data (see Ref. [35] for details).
Various parameters have been introduced to describe the theoretical uncertainties and approximations in the RFG+RPA+2p2h model, which are constrained using the near detector data.The variable parameters are the axial-mass, M A , the Fermi momentum, p F , the binding energy, E b , and the 2p2h cross-section normalization.Given the poor agreement between MINERvA and Mini-BooNE datasets [35], and slight inconsistencies between the signal definitions from the two experiments, no external constraints are applied on the variable parameters prior to the ND280 fit.Given the absence of firm predictions on the scaling of various nuclear effects with the nucleus mass number, the Fermi momentum, binding energy and 2p2h normalization are treated as uncorrelated between 12 C and 16 O.The 2p2h normalization is also considered to be uncorrelated between neutrino and antineutrino interactions.All of these parameters can be separately constrained by the T2K near detector data with the inclusion of samples with interactions on both 12 C and on 16 O, with ν-and ν-mode data.
The NEUT model for resonant pion production is based on the Rein-Sehgal model [46] with updated nucleon form factors [47], and with the invariant hadronic mass restricted to be W ≤ 2 GeV to avoid doublecounting pions produced through Deep Inelastic Scattering (DIS).Three variable model parameters are considered: the resonant axial mass, M RES A ; the value of the axial form factor at zero transferred 4-momentum, C A 5 ; and the normalization of the isospin non-resonant component predicted in the Rein-Sehgal model, I 1 2 .Initial central values and uncertainties for these parameters are obtained in a fit to low energy neutrino-deuterium single pion production data from ANL [48] and BNL [49] for the resonant pion production channels ν µ p → µ − pπ + , ν µ n → µ − pπ 0 and ν µ n → µ − nπ + .For the dominant production channel, ν µ p → µ − pπ + , the reanalyzed dataset from Ref. [50] was used.Resonant kaon, photon and eta production is also modeled using the Rein-Sehgal resonance production amplitudes, with modified branching ratios to account for the decay of the resonances to kaons, photons, or eta, rather than to pions.External neutrino-nucleus and antineutrino-nucleus pion production data from MiniBooNE [51][52][53], MINERvA [54], Sci-BooNE [55] and K2K [56] were used as a cross check to ensure that the broad features of all datasets were consistent with the uncertainties on the interaction level parameters (M RES A , C A 5 and I 1 2 component) or the uncertainties on Final State Interactions (FSI, which will be described shortly).A full fit to the external data is difficult due to strong correlations between FSI parameters and the neutrino-nucleus interaction model parameters, and a lack of information on the correlations between external data points.
The coherent pion production model used is the Rein-Sehgal model described in Refs.[57,58].However, recent results from MINERvA [59] are better described by the Berger-Sehgal model [60], so a rough reweighting of the coherent events as a function of the outgoing pion energy, E π , is applied to approximate the Berger-Sehgal model using the weights and binning given in Tab.II.Normalization uncertainties of 30% are introduced separately for CC-and NC-coherent events, based on comparisons to the MINERvA data (after the weights in Table II are applied), which are fully correlated between 12 C and 16 O.
The Deep Inelastic Scattering (DIS) model is unchanged from previous analyses (described in Ref. [27]).The DIS cross section is calculated for W ≥ 1.3 GeV, using GRV98 structure functions [61] with Bodek-Yang corrections [62].Single pion production through DIS is suppressed for W ≤ 2 GeV to avoid double counting with the resonant pion production contributions, and uses a custom hadronization model described in Ref. [28].For W > 2 GeV, PYTHIA/Jetset [63] is used for hadronization.A CC-other shape parameter, x CC−Other was introduced to give flexibility to the CC-DIS contribution.This parameter applies to CC resonant kaon, photon and eta production, as well as CC-DIS events, and it scales the cross section by (1 + x CC−Other /E ν ).It was designed to give greater flexibility at low E ν , and the initial uncertainty was set by NEUT comparisons with MINOS CC-inclusive data [64].
In addition to the previously described NC-coherent parameter, two other NC-specific parameters have been introduced in this analysis.A study in Ref. [65] showed that the NEUT neutral current single photon production (NC1γ) cross section prediction was approximately a factor of two smaller than a recent theoretical model [66].Because of this, the NC1γ cross section has been set to be 200% of the NEUT nominal prediction, with an uncertainty of 100%.An NC-Other normalization parameter is applied to neutral current elastic, NC resonant kaon and eta production, as well as NC-DIS events, with an initial uncertainty set at 30%.
Hadrons produced inside the nucleus may undergo FSI before leaving the nuclear environment, which changes the outgoing particle content and kinematics in the final state.NEUT models FSI for pions, kaons, etas and nucleons using a cascade model described in Ref. [28].Interactions are generated inside the nucleus according to a Woods-Saxon density distribution [67], and all outgoing hadrons are stepped through the nucleus with interaction probabilities calculated at each step until they leave the nucleus.Particles produced in DIS interactions are propagated some distance without interacting to allow for a formation zone, where the initial step size is based on results from the SKAT experiment [68].The allowed pion interactions in the nucleus are: charge exchange, where the charge of the pion changes; absorption, where the pion is absorbed through two or three body processes; elastic scattering, where the pion only exchanges momentum and energy; and inelastic scattering, where additional pions are produced.If an interaction occurs, new and modified particles are also added to the cascade.For pion momenta p π ≥ 500 MeV, nucleons are treated as free particles, and separate high (p π ≥ 500 MeV) and low (p π < 500 MeV) energy scattering parameters are introduced for charge exchange and elastic scattering.Initial interaction uncertainties are obtained from fits to a large body of pion-nucleon and pion-nucleus scattering data for nuclei ranging from carbon to lead, as described in Ref. [27].The variable parameters included to vary the pion FSI cross section are summarized in Tab.IX.Pion FSI parameters are assumed to be fully correlated between 12 C and 16 O.Uncertainties on nucleon, kaon and eta FSI interaction probabilities are not considered in the current analysis.
To account for effects which may potentially affect ν ( -) e but not ν ( -) µ cross sections, such as radiative corrections or second class currents (see, for example, Ref. [69]), which are not included in the NEUT cross section model, additional uncertainties which affect ν ( -) e have been introduced.These include an uncorrelated 2% uncertainty on the ν e /ν µ and νe /ν µ cross section ratios to account for radiative corrections, and an additional 2% uncertainty which is fully anticorrelated between ν e and νe to allow for second class currents.
The full list of cross-section uncertainties and their values before and after the ND280 data constraints is pro-vided in Tab.IX.

B. Alternative neutrino interaction models for studies of simulated data
Neutrino interactions with 12 C and 16 O nuclear targets at the near and far detectors may be affected by important nuclear effects which are not well understood.Various theoretical models are available to describe such effects, which are based on different approximations and with different ranges of validity.None of the available models are capable of describing all the available measurements of neutrino-nucleus cross sections from T2K and from other experiments.It is therefore crucial to test that the T2K oscillation analysis is insensitive to reasonable modifications of the neutrino interaction model described in Sec.III A, which will now be refered to as the "reference model".With this aim, various simulated datasets have been built based on alternative models.The following effects have been considered: variations of the distribution of the momentum of the initial nucleons in the nucleus and of the energy needed to extract the nucleons from the nucleus (the binding energy); uncertainties on the long-range nuclear correlations modifying the cross section as a function of Q 2 ; and modifications of the modeling of multi-nucleon interaction model (2p2h), including short-range nuclear correlations and meson exchange currents.
To test the nuclear effects in the initial state two alternative models have been considered beyond the RFG simulation used as reference: the SF developed in Ref. [39] and the Local Fermi Gas (LFG) model from Ref. [40].The LFG model also differs from the reference model in the implementation of the binding energy.In the latter, an effective value is considered, based on the average momentum of nucleons within the nucleus, while the LFG model considers the different state of the initial and final nucleus after the nucleon ejection, naturally including a different binding energy for neutrino and antineutrino interactions.The simulated datasets built with this alternative model will be referred to as the "alternative 1p1h model".
The correction to the CCQE cross section due to longrange nuclear correlations, described by RPA in the reference model, has been parametrized as a function of Q 2 in terms of five free parameters.A joint fit to the MiniBooNE [42,43] and MINERvA [44,45] ν µ and νµ datasets has been performed to extract an alternative, data-driven RPA correction, labeled "effective-RPA" in the following.The effective-RPA correction deviates from the reference RPA at high Q 2 , as can be seen in Fig. 3.
The model in Ref. [70] has been considered as an alternative 2p2h model, which differs from the reference model in many respects.The alternative 2p2h cross section is twice as large for neutrino interactions but has a similar strength for antineutrino interactions, except at FIG. 3. Effective-RPA with error band from the fit to external data compared with RPA corrections computed in Ref. [40].

high neutrino energies (E ν
1 GeV) where it is about 30% larger, as can be seen in Fig. 4. The difference between the 2p2h normalization for neutrino and antineutrino interactions is constrained with the ND280 data in order to avoid biases in the CP asymmetry measurement in the oscillation analysis.The alternative model has also been used for one of the studies of simulated data.Another important difference between the two models consists in the relative proportion of nucleonnucleon correlations, meson exchange currents and their interference, the first being strongly enhanced in the alternative model.This difference affects the estimation of the neutrino energy from the outgoing lepton kinematics.This estimation assumes the CCQE hypothesis, and it is well known that the 2p2h contribution biases the neutrino energy reconstruction [30,31] if not properly taken into account in the simulation.The reference model includes 2p2h events, and so this effect is included in the T2K neutrino oscillation analysis.Nevertheless, the different 2p2h components produce different biases in the neutrino energy estimation, as shown in Fig. 5, therefore incorrectly estimating the relative proportions of nucleon-nucleon correlations and meson exchange current can cause a residual bias in the neutrino energy estimation.To address this, three simulated datasets have been built.In the first, the multinucleon interactions have been reweighted as a function of neutrino energy, separately for neutrino and antineutrino, to reproduce the alternative model (referred to as the "alternative 2p2h model" in the following).In the other two simulated datasets, the full 2p2h cross section has been assigned either to meson exchange currents ("Delta-enhanced 2p2h") or nucleon-nucleon correlations only ("not-Delta 2p2h") by reweighting the muon kinematics as a function of muon angle, muon momentum and neutrino energy.

IV. THE ND280 COMPLEX
The precise measurement of neutrino oscillations in T2K requires a good understanding of the neutrino beam properties and of neutrino interactions.The two previous sections have described the neutrino flux model and neutrino-nucleus interaction model, and constraints on those models based on external measurements.As we will show in Sec.VIII D, with only that information, the precision on the measurements of neutrino oscillations parameters would be limited.In order to reduce the model uncertainties a near detector complex has been built 280 m downstream of the production target.The goal of the near detectors is to directly measure the neutrino beam properties and the neutrino interaction rate.The near detector complex comprises an on-axis detector (INGRID) and an off-axis detector (ND280).INGRID is composed of a set of modules with sufficient target mass and transverse extent to monitor the beam direction and profile on a day-to-day basis.The ND280 is composed of a set of subdetectors, installed inside a magnet, and is able to measure the products of neutrino interactions in detail.
In this section, the methods used to select high purity samples of neutrino and antineutrino interactions in IN-GRID and ND280 will be described, and the results are compared with the predictions obtained from the beam line simulation and the interaction models.The use of the ND280 data to reduce the systematic uncertainties in the T2K oscillation analysis will be described in Sec.V.

A. On-axis near detector
The INGRID detector is used to monitor the neutrino beam rate, profile and center.Those parameters are used to determine the off-axis angle at SK. INGRID is centered on the neutrino beam axis and samples the neutrino beam with a transverse cross section of 10 m × 10 m using 14 modules positioned in the shape of a cross.Each INGRID module holds 11 tracking segments built from pairs of orthogonally oriented scintillator planes interleaved with nine iron planes.There are also three veto planes, located on the top, bottom and one side of each module.The most upstream tracking plane is used as a front veto plane.The scintillator planes are built from 24 plastic scintillator bars instrumented with fibers connected to multi-pixel photon counters (MPPCs) to detect scintillation light.More details can be found in Ref. [71].

Event selection and corrections
Neutrino and antineutrino interactions within INGRID modules are first reconstructed independently in the horizontal and vertical layers of scintillators.Pairs of tracks in the two different orientations are then matched by comparing the most upstream point to form 3D tracks.The upstream edges of the different 3D tracks are then compared in the longitudinal and transverse direction with respect to the beam direction in order to construct a common vertex.The subsequent reconstructed event is rejected if the vertex is reconstructed out of the fiducial volume, if the external veto planes have hits within 8 cm from the upstream extrapolated position of a reconstructed track, or if the event timing deviates from more than 100 ns to the expected event timing.
In order to reduce the systematic uncertainty on the track reconstruction in ν-mode, the selection has been improved from the one used in Ref. [27].To reduce the impact of MPPC dark noise, the reconstruction is only applied to events where two consecutive tracking planes have a hit coincidence on their horizontal and vertical planes.This condition was not used in Ref. [27] and has been applied only to the ν-mode in the analyses presented here.A total of 12.8 × 10 6 and 4.1 × 10 6 neutrino events are reconstructed respectively in ν-and ν-mode, with estimated purities of 99.6% and 98.0% respectively.
The selected number of events in each module is corrected to take into account the impact of the detector dead channels, the event loss due to non-reconstructed neutrino interactions caused by pile-up, the variation of the iron mass between the modules, the time variation of the MPPC noise during the data taking and the contamination from external background as the previous IN-GRID analysis [27,71].

Systematic uncertainties
The systematic uncertainties on the event selection are estimated using the simulation and control samples.The sources of error are the same as those identified in Ref. [27] and include the neutrino target mass, the accidental coincidence with MPPC dark noise, the hit efficiency, the event pile-up, the cosmic and beam-induced backgrounds along with errors associated to event selection cuts.The method for estimating the uncertainty has not been changed since Ref. [27] for ν-mode, and is also applied here for ν-mode data.The uncertainties are evaluated to be 0.9% and 1.7% for neutrino and antineutrino data respectively.The larger uncertainty in ν-mode mainly arises from a discrepancy between data and simulation for interactions producing a track that cross less than four tracking planes.

Results of neutrino beam measurement
The stability of the neutrino flux is monitored by measuring the event rate, that is the total number of selected events per protons on target.Fig. 6 shows the intensity stability as a function of time for both ν-and ν-modes.Most of the data have been taken with the horn currents set to an absolute value of 250 kA, except for a small fraction of ν-mode data taken during T2K run 3 in which horns were operated at 205 kA.The average event rates are compared with the simulation, and the ratios = 0.984 ± 0.001(stat.)± 0.017(syst.). ( The quoted systematic uncertainties do not include the uncertainties on flux and cross section model, but they only include INGRID detector systematic uncertainties.The numbers of expected events in the Monte Carlo are obtained with the cross section models described in Sect.III A. The spatial spread of the neutrino beam is measured using the number of reconstructed events in each INGRID module.This produces a measurement of the number of events as function of the distance from the center in both vertical and horizontal directions.The two distributions are fit with a Gaussian, and the neutrino beam center and width are given by the mean and the sigma of the fit.The measurement of the position of the beam center is crucial to determine the off-axis angle, and therefore, the neutrino beam energy at SK.A deviation of 1 mrad of the beam direction would shift the peak neutrino energy by 2%.Fig. 6 shows the beam direction stability for all data taking periods.The variations are well within the design goal of 1 mrad.The average angles are θ beam, ν X = 0.027 ± 0.010(stat.)± 0.095(syst.)mrad θ beam, ν Y = 0.036 ± 0.011(stat.)± 0.105(syst.)mrad (2) for ν-mode, and θ beam, ν X = −0.032± 0.012(stat.)± 0.121(syst.)mrad θ beam, ν Y = 0.137 ± 0.020(stat.)± 0.140(syst.)mrad (3) for ν-mode.All values are compatible with the expected beam direction.

B. The off-axis ND280 detector
The off-axis near detector ND280 measures the neutrino energy spectrum, flavor content, and interaction rates of the unoscillated beam.These measurements are crucial to reduce the uncertainties on neutrino flux and interaction models which affect the prediction on the number of expected events at the far detector.
The ND280 detector consists of a set of subdetectors installed inside the refurbished UA1/NOMAD magnet, which provides a 0.2 T field, used to measure the charge and the momentum of particles passing through ND280.For the analyses described in this paper, ν µ and ν µ charged current interactions are selected in the tracker region of ND280, which consists of three Time Projection Chambers (TPC1, 2, 3) [72], interleaved with two Fine-Grained Detectors (FGD1, 2) [73].
The upstream FGD1 detector consists of fifteen polystyrene scintillator modules, while the downstream FGD2 contains seven polystyrene scintillator modules interleaved with six water modules.The FGDs provide target mass for neutrino interactions and track 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 of the charged tracks produced in neutrino interactions in order to measure exclusive CC event rates.The major updates in the near detector analysis with respect to Ref. [27] are the use of interactions in FGD2 and the inclusion of data taken with the ν-mode beam.
The charge and particle identification ability of the tracker is important because it provides separation between µ + (produced by νµ CC interactions) and µ − (produced by ν µ CC interactions) when T2K runs in ν-mode.Moreover, by including both FGD1 and FGD2 samples, the properties of neutrino interactions on water can be effectively isolated from those on carbon, reducing the uncertainties related to extrapolating across differing nuclear targets in the near and far detectors.The near detector analysis described here uses a reduced dataset comprising 5.81 × 10 20 POT in ν-mode and 2.84 × 10 20 POT in ν-mode, as shown in Tab.III.The event selection in ν-mode beam is unchanged since the previous analysis described in Ref. [27].Muonneutrino-induced CC interactions are selected by identifying the µ − produced in the final state as the highestmomentum, negative-curvature track in each event with a vertex in FGD1 (FGD2) Fiducial Volume (FV) and crossing the middle (last) TPC.The energy lost by the selected track in the TPC must be consistent with a muon.
All the events generated upstream of FGD1 are rejected by excluding events with a track in the first TPC.The selected ν µ CC candidates are then divided into three subsamples, according to the number of identified pions in the event: CC-0π, CC-1π + and CC-Other, which are dominated by "quasi-elastic" (CCQE), CC resonant pion production, and DIS interactions, respectively.Pions are selected in different ways according to their charge.A π + can be identified in three ways: an FGD+TPC track with positive-curvature and an energy loss in the TPC consistent with a pion; an FGDcontained track with a charge deposition consistent with a pion; or a delayed energy deposit in the FGD due to a decay electron from stopped π + → µ + .In this analysis, π − 's are only identified by selecting negative-curvature FGD+TPC tracks, while π 0 's are identified by looking for tracks in the TPC with charge depositions consistent with an electron from a γ conversion.The output of the ν-mode tracker selection are six samples, three per FGD.The selected CC-0π and CC-1π + samples in both FGDs before the ND280 fit are shown in Fig. 7.For each of the selected samples, the number of observed and predicted events are shown in Tab.IV.The main difference between ν-and ν-modes is the increase in the number of interactions produced by "wrongsign" neutrinos.Once differences in the flux and the cross section are taken into account, the wrong-sign contamination in ν-mode is expected to be approximately 30%, while the wrong-sign contamination in ν-mode is approximately 4%.
The lepton selection criteria of νµ ( ν µ ) CC interactions is similar to the one used in the neutrino beam mode, except for the condition that the highest-momentum, positively (negatively) charged particles must also be the highest momentum track in the event.This additional cut is essential to reduce the background due to π + (π − ) generated in neutrino (antineutrino) interactions that can be misidentified as the muon candidate.The selected νµ CC (ν µ CC) candidate events are divided in two subsamples: CC-1-Track, dominated by CCQE-like interactions; and CC-N-Tracks (N>1), a mixture of resonant production and DIS.These two subsamples are  defined by the number of reconstructed tracks crossing the TPC.For these selections the CC candidates are not divided into three subsamples as in Sec.IV B 1, according to the number of identified pions in the event in order to avoid samples with low statistics.The output of the ν-mode tracker selection are eight samples, four per FGD.For each of the selected sample, the number of predicted events and the ones observed in data are shown in Tab.V.The four selected samples in FGD1, before the ND280 fit, are shown in Fig. 8.

C. ND280 detector systematic uncertainties
In order to assess systematic uncertainties related to the ND280 detector modeling, various different control samples are used, as described in Ref. [27].The control samples include muons produced in neutrino interactions outside ND280, cosmic muons, interactions upstream of TPC1, and stopping muons.All control samples are independent of the samples used for the ND280 analyses described earlier.The method to propagate the systematic uncertainties in the near detector analysis is also un-TABLE V. Observed and predicted event rates for different ND280 samples collected in ν-mode beam.Before the ND280 fit that will be described in Sec.V uncertainties of 20% on the event rates are expected.new methods used to evaluate charge misidentification and FGD tracking efficiency uncertainties and uncertainties due to interactions outside the fiducial volume.The ToF between FGD1 and FGD2 is used to select events with a backward muon candidate in the FGD2 samples.The ToF systematic uncertainty is ∼0.1% for the ν-mode FGD2 samples, and ∼0.01%for the ν-mode FGD2 samples.The ToF uncertainty is smaller in ν-mode because fewer backward-going µ + are produced in ν µ interactions than backward-going µ − in ν µ interactions.The charge mis-identification uncertainty is parametrized as a function of the momentum resolution in the TPCs.
The FGD tracking efficiency for CC-events, where either a short pion or proton track is also produced, is estimated using a hybrid data-MC sample.This sample uses events with a long FGD-TPC matched muon candidate track with the addition of an FGD-isolated track generated via particle gun with a common vertex.
The method used to estimate the number of out-offiducial volume (OOFV) events has been refined by estimating the number of events, and the error, separately for each detector in which the OOFV events occur, rather than averaging over the number of OOFV events produced in all of the detectors outside the tracker as previously.
Most sources of systematic error are common between ν-and ν-modes because the selection criteria are similar, as described in Sec.IV B 2. However, as ν-mode data are divided into CC 1-Track and CC N-Track samples, only based on the number of reconstructed FGD-TPC matched tracks, most uncertainties relating to the FGD reconstruction are not relevant.The exceptions are the FGD-TPC matching and ToF uncertainties, which apply to both modes.Other differences between modes arise because some errors change with the beam conditions (sand muons, pile-up and OOFV) and are evaluated independently for each run period.
The total systematic uncertainties are shown in Tab.VI.
The dominant source of uncertainty for all ND280 samples comes from the pion re-interaction model, used to estimate the rate of pion interactions in the FGDs.This is due to differences between the GEANT4 model, used to simulate pion re-interactions outside the nucleus, and the available experimental data.For example, the system- atic uncertainty related to pion interactions affecting the FGD1 ν µ CC-0π (ν µ CC-1-Track) sample is 1.4% (4.9%), with a total error of 1.7% (5.4%).The pion re-interaction uncertainty is larger for ν-mode samples than for ν-mode samples because π − interactions on carbon and water are less well understood than π + interactions at the relevant energies, and because the fraction of π + from wrong-sign contamination in ν-mode misidentified as µ + candidate is larger than the fraction of π − misidentified as µ − in the ν-mode.

V. NEAR DETECTOR DATA ANALYSIS
The predicted event rates at both the ND280 and SK are based on parametrized neutrino flux and interaction models, described in Sec.II and Sec.III A. These models are fit to the precisely measured, high statistics data at the ND280, producing both a better central prediction of the SK event rate and reducing the systematic uncertainties associated with the flux and interaction models.The near detector analysis uses event samples from both FGD1 and FGD2, and from the ν-mode and νmode data, giving 14 samples in total.These, along with their associated systematic uncertainty, were described in Sec.IV B.

A. Near detector likelihood and fitting methods
The form of the ND280 likelihood and the fitting method are the same as described in Ref. [27].The 14 event samples are binned in p µ and cos θ µ , giving 1062 fit bins in total, though only the p µ projection is shown for clarity.The likelihood assumes that the observed number of events in each bin follows a Poisson distribution, with an expectation calculated according to the flux, cross-section and detector systematic parameters discussed above.A multivariate Gaussian likelihood function is used to constrain these parameters in the fit, with the initial constraints that are described in Secs.II, III A and IV B. The near detector systematic and near detector flux parameters are treated as nuisance parameters, as are the cross-section systematic parameters governing neutral current and electron neutrino interactions.The fitted neutrino cross-section and unoscillated SK flux parameters are passed to the oscillation analysis, using a covariance matrix to describe their uncertainties.
One significant difference with respect to Ref. [27] is that, as discussed in Sec.III A, the CCQE cross-section parameters (except the nucleon binding energy, E b ) have no external constraint.These parameters are constrained solely by the ND280 data.In addition, in order to alleviate possible biases on the estimation of the oscillation parameters (see Sec. IX for more details), the differences between the reference model and the alternative model for the 1p1h component of the neutrino-nucleus interaction cross section described in Sec.III are taken into account in the likelihood.This is done by adding the difference in the expected number of events between the two models in each p µ and cos θ µ bin to the diagonal of the ND280 detector covariance matrix V d .Finally, another significant difference is the inclusion of event samples from FGD2, which contains a water target, and ν-mode data samples.

B. Fit results
The fit produces central values for the flux, crosssection and detector systematic parameters along with a covariance.Fig. 9 shows the values of the unoscillated SK flux parameters and Fig. 10 the cross-section parameters before and after the fit as a fraction of the nominal value, along with their prior constraints.These parameter values are listed in Tables VII, VIII and IX, showing the best-fit point for each along with its uncertainty, calculated as the square root of the diagonal of the covariance.
Most noticeable in these results is the 10-15 % increase in the neutrino flux, seen across all species and energies in both ν-and ν-modes.
Small changes are seen in the central values of the CCQE cross-section parameters, with the fit increasing the Fermi momentum parameter while reducing the nucleon binding energy and axial mass parameters.More interestingly, the 2p2h normalization is increased to approximately 1.5 times its nominal value, indicating that the fit is sensitive to differences in lepton kinematics between CCQE and 2p2h interactions.The antineutrino 2p2h normalization is reduced compared to the neutrino parameter, highlighting a difference in the neutrino and antineutrino CC-0π event rates that cannot be explained by flux or detector systematics.The fit also reduces the value of the charged current single pion parameters, as seen in the previous analysis [27].This accounts for the relative deficit observed in the CC-1π sample compared to the CC-0π sample.

Goodness of fit and fit validation
The goodness of fit for the near detector analysis was estimated using mock datasets including statistical uncertainties.Mock datasets are generated by simultaneously varying the systematic parameters in the fit according to their prior covariance then applying these to the nominal MC.These were then fit and the minimum negative log-likelihood value found.The distribution of the minimum negative log-likelihood values is shown by the histogram in Fig. 11, with the value from the data fit indicated with a red line.The overall p-value for the fit is 8.6%.Figure 12 shows the same distribution for the flux and cross-section parameter priors, demonstrating that the fitted parameter values propagated to the oscillation analysis are reasonable.
In addition, the Bayesian analysis which simultaneously fits both near and far detector samples, that will be described in Sec.VIII B, was used to cross-check the primary result by only fitting near detector data.The results of this fit are compared to the best-fit parameters from the near detector analysis in Fig. 13, showing excellent agreement between the two.

C. ND280 postfit distributions
The expected muon momentum spectrum after the ND280 fit for the CC-0π and CC-1π + samples in ν-mode and the FGD1 samples in ν-mode are shown in Fig. 14 and Fig. 15 respectively.After the ND280 fit, the expected distributions show in general a better agreement with the data.The numbers of postfit predicted events for all the 14 samples are shown in Tab.X.The effects of the ND280 fit on the different neutrino interactions contributing to each ND280 sample are detailed in Tab.XI and Tab.XII.

VI. FAR DETECTOR EVENT SELECTION AND SYSTEMATICS
The T2K far detector is Super-Kamiokande (SK) [74], a 50 kton water Cherenkov detector located in the Kamioka Observatory, Gifu, Japan.SK is divided into two concentric cylinders, defining an inner detector (ID) instrumented with 11,129 20-inch photomultiplier tubes (PMT) and an outer detector (OD) instrumented with 1885 PMTs.The outer detector is mainly used as a veto for entering backgrounds while neutrino interactions are selected in a fiducial volume inside the ID.
In order to precisely measure neutrino oscillations parameters, together with the large target volume, high acceptance and efficient discrimination is necessary to distinguish the leptons produced in ν µ and ν e interactions.Vertex, momentum reconstruction, and Particle identifi-  cation (PID) in SK are done by observing the Cherenkov radiation produced by charged particles traversing the detector.These particles produce ring patterns that are recorded by the PMTs and are the primary tool used for the PID.Muons produced by ν ( -) µ CC interactions are usually unscattered thanks to their large mass and produce a clear ring pattern.In contrast, electrons from ν ( -) e CC interactions produce electromagnetic showers resulting in diffuse ring edges.In addition to the shape of the Cherenkov ring, the opening angle also helps to distinguish between electrons and muons.At the typical energies of leptons produced by neutrinos from the T2K beam, the probability to misidentify a single electron (muon) as a muon (electron) is 0.7% (0.8%).In SK it is not possible to distinguish neutrinos from antineutrinos event-by-event since the charge of the outgoing leptons cannot be reconstructed.For this reason the selection that will be described in this section is identical for data taken in ν-mode and in ν-mode.
T2K data are extracted from the incoming stream in ±500 µs windows centered on the beam trigger.A scan of each window recovers individual events which are then classified.
Events in which Cherenkov light is deposited exclu-sively in the ID comprise the fully contained (FC) sample.PMTs in the OD that register light are grouped into clusters.If the largest such cluster contains more than 15 PMTs the event is moved to an OD sample.Low energy (LE) events are separated into a dedicated sample by requiring that the total charge from the ID PMT hits in a 300 ns window be greater than 200 photoelectrons (p.e.), which corresponds to the charge observed from a 20 MeV electromagnetic shower.Events are also designated as LE if a single ID PMT hit constitutes more than half of the total p.e. observed to reject events due to noise.LE events are not included in the analysis presented here.Events at the far detector are timed with respect to the leading edge of the beam spill, taking into account the neutrino time of flight, the Cherenkov photon propagation time and delays in the DAQ electronics.Fig. 16 shows the event timing (∆T 0 ) distribution for all OD, LE, and FC events within ±500 µs of the beam arrival time.For FC events a visible energy (E vis ) greater than 30 MeV is also requested.A clear peak is observed around ∆T 0 = 0 in the FC sample.For an event to be incorporated into the analysis, ∆T 0 must lie in the interval [-2,10] µs.A total of 4 events have been observed outside this range.The expected number of such events is 4.17, estimated by using data taken with no beam.FC events within the spill window can be seen in Fig. 17 where the beam structure with eight bunches is clearly visible.The dotted lines represent the fitted bunch center times with a fixed bunch interval of 581 ns.A fiducial volume is defined within the ID, 2 m away from the detector wall, with a fiducial mass of 22.5 kt.Events whose vertex is reconstructed within this volume are selected into the fully contained fiducial volume sample (FCFV).Visible energy is defined as the energy of an electromagnetic shower that produces the observed amount of Cherenkov light.In total, 608 events are classified as FCFV.The expected number of background events from non-beam related sources in accidental coincidence is estimated to be 0.0145.

A. SK charged-current quasi-elastic selection
Charged-current interactions ( ν ( -) + N → l ± + X) in the narrow energy range of the T2K beam most commonly produce single-ring events at SK because most of the resulting particles, except for the primary lepton, do not escape the nucleus, or are below detection threshold.The energy of the incoming neutrino can be calculated assuming the kinematics of a CCQE interaction and neglecting Fermi motion: where E rec ν is the reconstructed neutrino energy, m i and m f are the initial and final nucleon masses respectively, and m i = m i − E b , where E b = 27 MeV is the binding energy of a nucleon inside 16 O nuclei.E l , p l and θ l are the reconstructed lepton energy, momentum, and angle with respect to the beam, respectively.The selection criteria for both ν ( -) e CC and ν ( -) µ CC events were fixed using MC studies before being applied to data.Events are determined to be e-like or µ-like based on the PID of the brightest Cherenkov ring.The PID of each ring is determined by a likelihood incorporating information on the charge distribution and the opening angle of the Cherenkov cone.The PID likelihood distribution for νmode FCFV single-ring events is shown in Fig. 18.The same criteria are applied to events observed for both ν-  and ν-mode data taking.ν ( -) e CC candidate events are selected using the criteria listed in Tab.XIV.The E vis requirement removes low energy NC interactions and electrons from the decay of unseen parents that are below Cherenkov threshold or fall outside the beam time window.The π 0 -like event rejection uses an independent reconstruction algorithm which was introduced in previous analyses [27].The cut 25 GeV is required, as above this energy the intrinsic beam ν ( -) e background is dominant.The numbers of events remaining in the neutrino and antineutrino beam data after successive selection criteria for a simulation sample produced with the oscillation parameters of Tab.XIII are shown in Tab.XIV.
After all cuts, 32 events remain in the ν e CC candidate sample and 4 in the ν e CC candidate sample, as shown in Fig. 19.Kolmogorov-Smirnov (KS) tests of the accumulated events as a function of POT are compatible with a constant rate, with p-values of 0.99 and 0.25 respectively.
The vertex distributions of the candidate event samples are checked for signs of bias that might suggest background contamination.Fig. 20 shows the vertex distribution of the ν e CC candidate events in the SK tank coordinate system.Combined KS tests for uniformity in r 2 and z yield p-values of 0.1 and 0.6 for the ν e and ν e samples respectively.ν ( -) µ CC candidate events are selected using the criteria shown in Tab.XV.The momentum cut rejects charged pions and misidentified electrons from the decay of unobserved muons and pions.Fewer than two Michel electrons are required to reject events with additional unseen muons or pions.After all cuts are applied, 135 events remain in the ν µ CC candidate sample and 66 in the ν µ CC candidate sample as shown in Fig. 21.

B. SK charged-current single pion selection
A new far detector event sample has been included in the oscillation analysis described here.As mentioned previously, single-ring events produced by quasi-elastic interactions are the most common at T2K neutrino energies.By modifying the event selection criteria to include an additional Michel electron candidate, it is possible to select events with a pion produced below Cherenkov threshold.Michel electrons are tagged at SK by searching for secondary hit clusters within the same trigger window and in subsequent trigger windows that pass time, charge and vertex goodness criteria.Note that this selection tags only positive π + 's as π − 's are absorbed by nuclei before they decay when stopped in water.Fig. 22 shows the true Michel electron momentum distribution for true CC1π + simulated events reconstructed inside the fiducial volume.Fig. 23 shows the true pion momentum distribution for selected signal events and the selection efficiency.The Cherenkov threshold for charged pions is also shown.The efficiency falls above this threshold as the pion produces more light and the event-fitting algorithm is increasingly likely to find a second ring, thus disqualifying events from the sample.The event reduction for the full selection is shown in Tab.XVI.There is a larger background from misidentified muons in this sample compared with the single-ring selection as such events are more likely to contain a Michel electron candidate.It can be seen in Tab.XVI that there is an apparent difference between the data and the MC expectation after cut number five.Applying an equivalent cut sequence to atmospheric neutrino data and MC (in this case the neutrino direction is not known so a cut on E vis is used instead of E rec ν ) yielded no such discrepancy.This single-ring CC1π + selection has been implemented as a new sample for ν e appearance with neutrino beam data.The reconstructed energy equation is modified from the CCQE case by recognizing that the outgoing baryon is a ∆ ++ instead of a proton and neglecting nuclear effects: where m ∆ ++ is the mass of the ∆ ++ (1232.0MeV/c 2 ).constructed in the data while 3.1 events are expected for the oscillation parameters of Tab.XIII.Fig. 26 shows the vertex distribution of the ν e CC1π + candidate events in the SK tank coordinate system.

C. SK detector systematic uncertainties
This section discusses the estimation of the uncertainty in the selection efficiency and background for the oscil-lation samples that result from modeling of the SK detector.This topic has been covered in detail in previous publications [27] but there have been a number of updates, particularly related to the addition of the ν e CC1π + sample.
Control samples unrelated to the T2K beam are used to assess the uncertainties.Cosmic-ray muon samples are used to estimate uncertainties related to the FC, fiducial-volume and decay-electron requirements, for the selections of both ν ( -) e and ν ( -) µ CC candidates.The error  from the initial FC event selection is negligible.The uncertainty in the fiducial volume is estimated to be 1% using the vertex distribution of cosmic ray muons which have been independently determined to have stopped inside the ID.The uncertainty due to the Michel electron tagging efficiency is estimated by comparing cosmic-ray stopped muon data with MC.The rate of falsely identi-fied Michel electrons is estimated from MC and a 100% uncertainty on that rate is assumed.Other studies of systematic uncertainty in SK modeling divide simulated events into categories according to their final state topologies, with the criteria shown in Tab.XVII.These topologies do not correspond exactly with true interaction modes due to subsequent in- teractions within the nucleus or with neighboring nuclei or because one or more particles are produced below Cherenkov threshold.The dominant uncertainties are described in the following paragraphs.Atmospheric neutrino data are used to assess possible mismodeling of the ring counting (RC), particle identification (PID), and π 0 rejection for the first four topologies shown in Tab.XVII.In addition to flux and cross-section uncertainties from the atmospheric neutrino analysis, additional parameters are included in the fit to alter the cut values applied to the  The energy for the CCQE-like sample is reconstructed using Eq. 4. For the CC1π + sample the modified interaction kinematics in Eq. 5 are used.The expectation is based on the parameters of Tab.XIII.
MC for the three classifiers above, thereby allowing for possible mismodeling of the events.Separate parameters are used for each of the four final state topologies and for different ranges of visible energy.A likelihood is defined comparing the data and MC which is then marginalized over flux and cross section parameters using a Markov chain Monte Carlo to estimate corrected efficiencies for the four final state topologies in bins of E vis and their covariance.The measured shifts between the nominal and fitted efficiencies are included in the final uncertainty assuming full correlation between samples.
To estimate the uncertainty in modeling π 0 's we con-struct a set of hybrid data-MC control samples.Events in these samples are created by overlaying a single electronlike ring from the SK data with a simulated photon ring.The kinematics of the photon ring are chosen such that when combined, the two rings follow the decay kinematics of π 0 events from the T2K MC.Hybrid samples with both rings from the SK MC are also produced for comparison with the hybrid data.The difference in the selection efficiency when using the oscillation analysis sample candidate criteria is used to determine the systematic error.This procedure is performed twice, with the data electron as the higher energy and lower energy ring, and the In previous analyses, the CC1µ events in the ν e ap-pearance sample that did not involve muon decay-inflight (81% of CC1µ) were assigned a conservative 150% error due to the ring-counting, PID and π 0 rejection cuts.However while these events form only ∼1% of the ν e CCQE-like candidate sample, they are ∼10% of the ν e CC1π + candidate sample and thus the uncertainty was re-evaluated with a dedicated analysis using stopped cosmic-ray muons and atmospheric neutrinos.Based upon this work, this uncertainty has been reduced to 63.2%, where the error on the particle identification is the dominant source of uncertainty.For CC1µ events involving a muon decay-in-flight, a 16% uncertainty is assigned, unchanged with respect to the previous analysis.
All aspects of the SK detector simulation that can affect the modeling of the candidate event selection described above are propagated using a vector of systematic uncertainty parameters, s, which scale the nominal expected number of events in bins of the observable kinematic variables E rec ν or p l for the true neutrino interaction mode categories.

VII. NEUTRINO OSCILLATION FRAMEWORK
The previous sections have described how the near detector samples are used to reduce the uncertainties on the neutrino fluxes before oscillation, and on the neutrino interaction model parameters.The best-fit values of the neutrino and antineutrino flux and cross-section parameters together with their correlations are extrapolated to SK, where the oscillation probabilities are measured.In this section, the general model that describes the oscillations between the three standard neutrino flavors, ν e , ν µ and ν τ , are described, with a focus on the relevant oscillation channels for the T2K experiment.Finally, external measurements of neutrino oscillation parameters, which are used as prior constraints in some fits, are summarized.

A. Oscillation probabilities in the PMNS framework
As anticipated by Pontecorvo, the neutrino flavor eigenstates do not correspond to mass eigenstates, but are linear superpositions of them: where only the three active neutrinos, k = 1, 2, 3 are considered.U αk is an element of the 3 × 3 unitary matrix, called the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix [76,77], which can be parametrized as: Vertex X (cm) where s ij ≡ sin θ ij and c ij ≡ cos θ ij , θ ij are the three mixing angles and δ CP is the CP-violating phase.The Majorana phases are neglected here as the three-flavor oscillation probability is invariant under their rotation.
The neutrino oscillation probability is determined by six parameters: the mixing angles θ 12 , θ 13 , θ 23 , which define the amplitude of the oscillation probability; the differences in the squared masses of the eigenstates, ∆m , that define the oscillation frequency and the position of the oscilla-tion maxima as a function of L/E; and the CP-violating phase, δ CP .While it is known that ∆m 2  21 > 0, the neutrino mass ordering (MO) and hence the sign of ∆m 2   32   has not yet been determined.
Effects due to charged-current coherent scattering of ν e with the electrons in the Earth matter can affect the neutrino event rate.The matter effect potential is proportional to +N e for ν e and −N e for νe , where N e is the matter electron density.
Neglecting matter effects, the ν ( -) µ survival probability can be written as but the θ 23 octant, i.e., whether θ 23 > π/4 or θ 23 < π/4, cannot be distinguished.At the T2K baseline and peak neutrino energy, the impact of the matter effect is small, and a full calculation including them would only modify the P (ν µ → ν µ ) and P (ν µ → νµ ) probabilities by ∼0.1%.
The ν ( -) e appearance probability, approximated to first order in the matter effect (see [78]), is where "(−)" corresponds to the change of sign required for antineutrino oscillations.ρ is the average density of the Earth matter through which the neutrinos travel and α eV 2 /c 4 = 7.56 × 10 −5 ρ g/cm 3 E ν [GeV].In all oscillation analyses described here, the exact three-flavor oscillation probabilities are used, rather than the approximations shown in this section.
As well as being proportional to sin 2 θ 13 , the first, leading, term of Eq. 10 is proportional to sin 2 θ 23 , which makes the appearance probability (in contrast to the disappearance probability) sensitive to the octant of θ 23 .The precise determination of θ 23 and the solution of the octant degeneracy are important for the determination of δ CP .Because of the precise measurement of θ 13 from reactor experiments (see Tab. XVIII), the largest uncertainty in the oscillation model which affects the determination of δ CP is due to the uncertainties on the value of θ 23 .
By looking at the other terms in Eq. 10, it is easy to understand why the ν ( -) e appearance transition is the golden channel in the search for CP violation.The term containing sin δ CP in Eq. 10 has an opposite sign for neutrinos and antineutrinos, and is the only term that can violate CP symmetry.Since all the mixing angles have been measured as nonzero, the requirement for violating CP symmetry in neutrino oscillations is sin δ CP = 0, i.e. δ CP = 0, π.
In T2K, the term proportional to sin δ CP can change the appearance probability by as much as ±30%.This means that an extreme value of δ CP such as δ CP = −π/2 would increase (decrease) the ν e (ν e ) by 30% with respect to δ CP = 0.In Eq. 10 there are also other terms containing cos δ CP .While these terms do not violate CP symmetry, they change the shape of the appearance spectrum and are important for a precise determination of the value of δ CP .
In summary, at first order, the term in sin δ CP defines the normalization of the ν ( -) µ → ν ( -) e spectrum, while the term in cos δ CP provides information on the shape and allows the discrimination of the CP-conserving phases δ CP = 0 and π.
As mentioned above the appearance probability is also affected by the matter effects.Matter effects can enhance either ν µ → ν e , or νµ → νe depending on the sign of ∆m 2 32 .For instance if the mass ordering is normal (inverted) ν µ → ν e oscillations are enhanced (suppressed) while νµ → νe oscillations are suppressed (enhanced).
With its low neutrino energy and baseline of 295 km, the T2K experiment is not sensitive to the matter effect.The two possible solutions for the mass ordering (MO) change the appearance probability by approximately ±10% and a definitive measurement of the MO cannot be made by exploiting the matter effects in T2K beam.The interplay between effects due to δ CP and effects due to the MO on the ν ( -) e -appearance probabilities are shown in Fig. 27.Each ellipse shows the effect due to δ CP on the two appearance probabilities and the two ellipses represent the two possible solution for the mass ordering.The two regions overlap for more than half of the δ CP range, indicating the large degeneracy among the MO and δ CP at the T2K baseline.
The goal of T2K is the precise measurement of θ 13 , θ 23 , |∆m 2  32 | and to search for evidence of CP violation.Since T2K has excellent sensitivity to the absolute value of ∆m 2  32 but not to its sign, the absolute value of ∆m 2   32   and the MO are considered as two separate parameters.As mentioned above, the violation of CP symmetry would manifest itself as a difference between neutrino and antineutrino oscillation probabilities.In order to find a hint of this violation it is necessary to directly compare ν µ → ν e and νµ → νe oscillation probabilities, as also shown in Fig. 27, where the anti-correlation as a function of both δ CP and MO is visible.The measurement of the ν ( -) µ survival probabilities is important for determining |∆m 2  32 |, i.e. the position of the first oscillation maximum, and sin 2 2θ 23 , drastically reducing the uncertainty on δ CP .The degeneracy between δ CP and the θ 23 octant can be solved by comparing ν e and νe appearance probabilities.Indeed, in contrast to δ CP , sin 2 θ 23 has the same effect on both the ν µ → ν e and the νµ → νe oscillation probabilities.For this reason, T2K measures the ν µ → ν µ (ν µ → νµ ) disappearance and the ν µ → ν e (ν µ → νe ) appearance probabilities in ν-and ν-mode simultaneously in the analyses described in this paper.The joint analysis also allows all of the correlations between the oscillation parameters to be properly taken into account.

B. External constraints on the oscillation parameters
For some of the oscillation parameters, Gaussian priors, based on external experiments, are used in the oscillation analyses.These parameters are the solar oscillation parameters, to which T2K is insensitive, and, in some cases, the value of sin 2 θ 13 as measured by reactor experiments.In the latter case, the weighted average of Daya Bay, RENO and Double Chooz as reported in Ref. [75] is used.The Gaussian priors used in the analysis are shown in Tab.XVIII.For all the results shown in this paper, it will be specified whether reactor constraints on sin 2 θ 13 are used.
All the parameters that are measured by T2K are included in the oscillation analyses with flat priors which extend beyond the allowed regions for those parameters.The boundaries on the flat priors used in the analysis are 0.3 < sin 2 θ 23 < 0.7, 2<∆m 2  32 (×10 −3 eV 2 /c 4 )<3, −π < δ CP < π, and 0 < sin 2 (2θ 13 ) < 0.4 (when no reactor measurements are used).In this section, the oscillation analyses which have been developed to provide oscillation parameter estimations are described: three different analyses have been developed within frequentist and Bayesian frameworks.
For all analyses, the appearance and disappearance channels in both ν-and ν-modes are analysed simultaneously, incorporating constraints from several T2K samples and external experiments.The oscillation analyses compare the rate and distribution of events in analysis bins defined by combinations of the reconstructed neutrino energy, E rec ν (see Equation 4), the reconstructed lepton momentum, and the reconstructed angle of the outgoing lepton with model predictions.Systematic uncertainties are characterized by systematic model parameters, which tune the model prediction, and by a covariance matrix, which determines the sizes and correlations of their uncertainties.The fit to the near detector data is applied as a multivariate Gaussian penalty term to constrain flux and cross-section uncertainties common to the near and far detectors.
For a given set of systematic and oscillation parameters, the expected distribution is computed in the following way.The expected SK MC samples are produced assuming nominal systematic parameters without oscillation.The effect of oscillations are then included by applying the neutrino oscillation probability to events defined by neutrino energy and interacting flavor.In order to do this all the analyses use the full three-flavor oscillation probabilities.The number of events is weighted by the oscillation probability and for the systematic parameters, as will be described in Sec.VIII D. Then the expected number of events in each bin of the far detector samples are extracted from a Poisson distribution, leading to a likelihood of the form where o and f are particular choices of the vector of the neutrino oscillation parameters (free in the fit), and of the nuisance parameters.n obs i is the observed number of events in the ith analysis template bin, and ) is the expected number of events in the ith bin.
The method described in Ref. [79] is used to project out the nuisance parameters.The high dimensionality of the likelihood function is solved by computing the marginal likelihood.This marginal likelihood for the parameters of interest is found by integrating the product of the likelihood function and the priors over all the parameters that are not of interest, and defined as: This defines the likelihood which is a function only of the parameters of interest o, given a dataset, x, and a model prediction defined by a set of systematic and oscillation parameters.The marginal likelihood of Eq. 12 can be used to perform both a hybrid frequentist-Bayesian analysis as well as a full Bayesian analysis.In both cases the prior p.d.f. for flux, cross section and SK detector parameters is a multivariate Gaussian, with mean µ and covariance V .All the analyses have been validated with simulated datasets using different neutrino interaction models and different sets of oscillation parameters for the generation of the simulated samples.Examples of the validations will be shown in Sect.IX.

A. Frequentist analyses
Two analyses have been developed in the frequentist framework.The main difference between the two analysis is that one uses 2D templates in {E rec ν , θ lep } for the e-like samples and the other uses {p lep , θ lep } templates.The reasons for these choices will be described in Sec.VIII C. Both analyses construct confidence intervals or regions by fitting the SK dataset, using the likelihood ratio −2∆ ln L( o) = −2 ln[L( o; x)/L max (x)] as the test statistic, where L( o; x) is as defined in Eq. 12.If the Gaussian approximation holds, Eq. 12 converges to a χ 2 distribution, and the negative log-likelihood ratio is sometimes denoted as ∆χ 2 ( o).
Due to the difficulties in providing confidence intervals with high dimensionality, a subset of the oscillation parameters is projected out by using the marginalization method described above, with the priors shown in Tab.XVIII.
Two approaches can be used in order to calculate the confidence intervals.In the Gaussian approximation, where the parameter dependence is linear, the socalled "constant ∆χ 2 method" provides the correct coverage: the confidence region consists of the set of oscillation parameter values that satisfies the condition −2∆ ln L( o) ≤ X crit ( o).The precomputed values of X crit ( o) [75], that define the α% confidence level for different dimensionality, do not depend on the value of o.
Since in the neutrino oscillation analysis the gaussian approximation condition does not always hold, a toy MC method, recommended by Feldman and Cousins [80], should be used to define the values of X crit ( o) as a function of o, in order to provide the correct coverage.An ensemble of 10,000 toy datasets was produced for each point of the oscillation parameter grid, fine enough to describe the variation of X crit ( o) as a function of the free parameters.The flux, cross-section and detector systematic parameters were sampled from the same multivariate Gaussian with mean µ and covariance V used as prior p.d.f. in the marginalization method, as described above.The oscillation parameters constrained by external measurements, such as sin 2 θ 12 , ∆m 2  21 and sin 2 θ 13 , were sampled by using the p.d.f.'s with parameters given in Tab.XVIII.The remaining oscillation parameters were sampled using the distribution of −2∆ ln L( o) obtained from the fit of the T2K data.This method will be referred as "posterior-predictive".
Since the Feldman-Cousins method is computationally intensive, it is used only for the most important result of the oscillation analysis, i.e. the confidence intervals as a function of δ CP and MO.For all the other results the constant ∆χ 2 method, which is more practical and provides good coverage to a first order approximation, is used.

B. Bayesian analysis
The Bayesian analysis performs the integration of Eq. 12 using a Markov Chain Monte Carlo (MCMC) method.The Metropolis-Hastings algorithm [81] is used to populate the space of the oscillation and nuisance parameters, distributed according to the posterior probability density function.
This algorithm uses a weighted random walk to explore the parameter space with a chain of many points, X i .A new proposed step to X i+1 is accepted with a probability equal to the minimum between 1 and P (X i+1 )/P (X i ), i.e., it is accepted with 100% probability when P (X i+1 ) > P (X i ) and with a probability equal to P (X i+1 )/P (X i ) when P (X i+1 ) < P (X i ).When a step to X i+1 is rejected, the point X i is repeated in the chain, and a a new step is proposed again from the previous step X i .Typically, a chain consists of 10 6 − 10 7 steps.The population of accepted points in parameter space thus constructed is distributed according to the posterior probability density.
Credible intervals are built by selecting the highest density region in the space of the parameters of interest, typically one or two, that contains α% of the points.A kernel density estimator (KDE) [82,83], is used to extract the most probable values of the oscillation parameters of interest, A KDE PDF is built by smearing the step points with a Gaussian function which has a width inversely proportional to the number of points.The most probable values of the parameters of interest were found by maximizing the KDE PDF using MINUIT [84].

C. Analysis templates
The three analyses use different far detector templates.All of them use the reconstructed neutrino energy, E rec ν , for the ν µ and νµ candidate samples, while different templates are used for the ν e and νe candidate samples.One of the two frequentist analyses and the MCMC Bayesian analysis use E rec ν and the angle between the lepton and the neutrino beam direction, θ lep .The choice of these two parameters is due to the fact that E rec ν is important to infer the oscillation probability function, which depends on the true neutrino energy.θ lep is used because it provides additional separation between ν and ν candidate events.The ν and ν cross sections differ as a function of Q 2 , and leptons produced in ν interactions tend to be emitted at larger angles than those produced in ν interactions.
As was shown in Tab.XIV, while the number of νµ → νe in ν-mode is negligible, in ν-mode, the number of oscillated ν e is roughly 30% of the number of oscillated ν e and thus it is crucial to use templates which are able to distinguish between the two cases.In Fig. 28, the expected {E rec ν , θ lep } distribution in the single ring e-like sample in ν-mode for the signal (ν µ → νe ) and for the oscillated ν e background is shown.
In Fig. 29 the 2D templates for the three ν e candidate samples are shown.The ν e CC1π + sample has no events with reconstructed energy below 400 MeV due to the ∆ ++ production requirement in Eq. 5 used to reconstruct the neutrino energy.
The second frequentist analysis also uses 2D templates but instead of using E rec ν , it uses lepton momentum, p lep , and θ lep .This choice provides similar ν-ν separation as {p lep , θ lep } templates are directly related to E rec ν .Indeed, a particular value of E rec ν corresponds to a slice in the {p lep , θ lep } phase space.For this reason, {E rec ν , θ lep } and {p lep , θ lep } analysis templates provide very similar sensitivities to the oscillation parameters.
Finally, another difference among the three analyses is that the Bayesian analysis performs a simultaneous fit of the ND280 and SK datasets in order to validate the extrapolation of the flux, cross-section and detector systematic parameters from the near to far detector.The other two analyses perform a fit of the far detector data using the priors obtained by the near detector data fit for the flux and cross section systematic parameters.The two methods show very good agreement in the estimation of flux and cross-section parameters as was shown in Sec.V.

D. Systematic uncertainties on the oscillation analyses
The systematic uncertainties of the oscillation analysis are separated into four different categories: flux, cross section, and SK detector and momentum scale parameters.The flux and many of the cross-section parame- ters (mentioned in Tab.IX) are correlated by the ND280 measurements as shown in Sec.V, effectively reducing the systematic uncertainties of the measurements at SK.
The flux parameters assign the systematic uncertainties on the neutrino and antineutrino flux.These parameters are applied as normalization factors for different neutrino flavor subcategories and true energy bins.The flux uncertainties at SK are constrained by the ND280 measurements.
The cross-section parameters are applied based on the true neutrino interaction category of each event.Most of these parameters are normalization factors, but some also contain shape information, modifying the kinematic distributions in p lep − θ lep and E rec ν probability density functions for corresponding oscillation analysis samples.
There are several categories of SK uncertainty: detector efficiency, final state interaction (FSI), secondary interaction (SI), and photo-nuclear (PN) effect uncertainties.The update on the SK detector efficiency uncertainty related to the addition of the CC1π e-like sample is described in Sec.VI C. In addition, the FSI, SI, and PN effect uncertainties on the CC1π + e-like sample are evaluated using the same method introduced in Ref. [27].The increased uncertainties for this sample is mainly due to the larger backgrounds affecting it, and to larger uncertainties on the pion FSI and SI as shown in Tab.XIX, where the SK systematic uncertainties are separated into the SK detector and the FSI+SI+PN contributions.
The SK energy scale uncertainty is applied independently from other parameters.The energy scale uncertainty is applied as a normalization of E rec ν for each event, which may vary the total event rate by shifting the events into the cut regions of the visible energy and reconstructed neutrino energy selection criteria.The SK energy scale uncertainty is estimated to be 2.4%.
The effect of the systematic uncertainties on the predicted event rates of the ν-and ν-mode samples are summarized in Tab.XIX and Tab.XX respectively.The 1σ uncertainties are obtained by throwing large sets of toy experiments, varying only the selected systematic parameters for each experiment, and calculating the relative uncertainties from the distributions of the event rates.The anti-correlations between flux and cross section parameters reduce the systematic uncertainties when both these sources are taken into account.

IX. SENSITIVITY OF OSCILLATION PARAMETERS TO NEUTRINO INTERACTION MODELING
The neutrino interaction uncertainties are one of the main contributors to the systematic uncertainty on all oscillation measurements and there is a global effort underway to improve the understanding of neutrino cross sections.This has lead to the creation of a number of interaction models which can, at least partially, describe existing cross-section data but which also produce different predictions for the oscillated event rates and spectra at SK.This analysis uses the high statistics near detector data to constrain both the flux and cross-section model uncertainties, improving the prediction of the far detector event rate, and reducing the uncertainty on that prediction.However, the near and far detectors observe different neutrino energy spectra (mainly due to oscillations) so the underlying neutrino interactions they sample will be different.The different design of the two detectors also means that they are sensitive to different regions of the lepton kinematic phase space.The near detector fit therefore tunes the flux and cross-section models to a set of neutrino interactions and a lepton kinematic phase space that is not the same as that observed at the far detector.These differences could then be incorrectly attributed to neutrino oscillation effects in the oscillation analysis, which would result in biased oscillation measurements.
This has been studied in a phenomenological context using a simplified oscillation analysis, summarized in Ref. [29].The results of the studies show that longbaseline experiments may be biased by cross-section model choices, as might be expected from the qualitative arguments above.Given the potential impact of these model choices it is important to investigate these issues using the full T2K oscillation analysis framework.This has been performed for a range of neutrino interaction model variations, which were discussed in Sec.III B.

Production and analysis of simulated datasets
The analysis is performed using replica datasets created by changing one part of the nominal MC.Simulated data are created by applying the event selection described in Sec.IV B and Sec.VI to the nominal T2K MC.A weight is then applied to each selected event, calculated as the ratio of the altered interaction model to the nominal cross section model.For the SK simulated data, the relevant oscillation probability is also applied.This produces event samples corresponding to the alternative interaction model.An example of this is shown in Fig. 30 for the ND280 samples, where the left plot shows the ratio of the alternative 2p2h simulated data to the nominal MC for the FGD1 CC-0π sample while the right shows the ratio for the alternative 1p1h model.
For each interaction model variation, simulated data are generated at both the ND280 and SK.The nominal flux and interaction models are then fit to the ND280 simulated data using the procedure described in Sec.V.The result is then used as the ND280 input to the oscillation analysis, described in Sec.VIII A. The SK simulated data are fit to produce a set of likelihood contours for the oscillation parameters.These likelihood contours are then compared to the expected likelihood contour for fits where the nominal MC is used as data.For each oscillation parameter a bias is calculated, defined as the difference in the parameter best-fit point between the simulated data fit and the nominal data fit, divided by the 1σ uncertainty on that parameter from the nominal fit.

Results
Results from the ND280 fit are used for central values and uncertainties for the neutrino flux and cross-section model parameters.When fitting the nominal MC, the parameter estimators are found to be unbiased while their uncertainties are reduced.For the fits to simulated data generated from alternative models the values of both the flux and cross-section parameters change, as shown in Fig. 31.These show the results obtained by fitting the data generated using the alternative 2p2h model, rather than the nominal model.Figure 31 demonstrates how differences between the true and model interaction cross sections can be absorbed by both the flux and crosssection parameters.The ND280 fit result is used to predict the unoscillated event spectra at SK, producing the ν-mode ν µ and ν e event spectra shown in Fig. 32.The predicted spectra falls between the nominal and 2p2h models and the uncertainty typically covers much of the difference between the models.
To evaluate possible bias in the oscillation parameter estimators, a fit is performed with SK simulated data.The likelihood surfaces for the oscillation parameters of interest obtained by this fit are shown in Fig. 33.These results use the true parameter values listed in the "maximal oscillation" column of Tab.XXI.Fig. 33 shows that if the alternative 2p2h model is correct, using the Nieves model in the nominal MC produces small (∼2%) biases in the measured values of sin 2 θ 23 and ∆m 2  32 (< 1%).The study also shows a change in the δ CP likelihood contour, increasing the exclusion around δ CP = π/2 by 1.2 units of ∆χ 2 .Overall, as seen in the bottom right of Fig. 33, the effect on the appearance analysis is small, compared to the statistical uncertainty.
Figure 34 shows fits to the same alternative 2p2h simulated data as in Fig. 33, but using the non-maximal  oscillation parameters from Tab. XXI.For non-maximal disappearance, a larger change is observed in the sin 2 θ 23 likelihood contour.The effect on the point estimate and 68% CL interval is relatively small.Fig. 35 shows the results of the fits to alternative 1p1h simulated dataset described in Sec.III B. The alternative 1p1h model changes the T2K MC as a function of both neutrino energy and the angle at which the lepton is produced relative to the neutrino direction.Similarly to the alternative 2p2h case, this model variation introduces a small bias in the measurement of ∆m 2  32 and has a small effect on the δ CP likelihood contour.

Summary
Detailed studies have been performed to evaluate the sensitivity of the T2K oscillation analysis to neutrino interaction modeling.By using the near detector data to estimate flux and interaction model parameters, the oscillation parameter estimates are found to be essentially unaffected by changes to the interaction model.The largest observed effect on the oscillation parameter likelihood contours is shown in Fig. 34, while Fig. 35 is more representative of the majority of the model variations discussed in Sec.III B.
A summary of the maximum bias observed for all of the alternative models studied is shown in Tab.XXII.The bias is presented as a fraction of the expected 1σ uncertainty on each oscillation parameter, and is the maximum bias seen from all true oscillation parameter values tested.T2K does not expect a significant constraint on δ CP given the integrated POT available for this analysis.As a result, it is difficult to quantify the effect on δ CP in a single number, and so this parameter is not included in the summary table.For reference, the largest effect observed for δ CP is shown in Fig. 34.

X. νe APPEARANCE ANALYSIS
T2K has already observed ν e -appearance in a ν µ beam [7], but no direct evidence of ν e -appearance has been reported so far.A search for this process, with the data collected by T2K in ν-mode, has been performed using the analysis tools described in Sec.VIII.
In order to look for ν e -appearance independently from ν e -appearance, a parameter β is introduced which multiplies the ν e -appearance probability In this analysis, β can have two values: β = 0, which corresponds to no ν µ → ν e oscillations; and β = 1, which corresponds to the appearance probability as predicted by the PMNS framework.
The two hypotheses are tested using rate-only and rate-plus-shape analyses.The test statistic is different in the two cases.The test statistic used in the rateonly analysis is the number of e-like candidates observed at SK in ν-mode, while for the rate-plus-shape analysis, the likelihood ratio for the two hypotheses, assumining PMNS oscillations and no ν µ → ν e is used: A toy ensemble is generated according to the prior knowledge on the oscillation parameters defined in Sec.VII B for both values of β.The T2K information from the ν µ -disappearance, ν µ -disappearance and ν e -appearance1 channels are taken into account by using a posterior predictive method [85].In this method, the data from these channels are used to constrain the prior parameter space in generating the toy ensemble for the ν e sample while preserving their correlations.
The sensitivity of the analysis is computed by producing a simulated sample without statistical fluctuations using the values of the oscillation parameters defined in Tab.XIII.The test statistics for this particular dataset are 6.28 expected ν e candidates for the rate-only analysis, corresponding to the total number of expected events in ν-mode in Tab.XIV, and ∆χ 2 = 2.54 for the rateplus-shape analysis.When compared to the generated toy ensembles for the no-ν e appearance hypothesis the p-values are 0.047 for both, the rate-only and rate-plusshape analyses.When compared with the toy ensemble generated with β = 1, p-values of 0.52 and 0.41 are found for the rate-only and rate-plus-shape analyses respectively.
The same analysis is then applied to the e-like data selected in ν-mode.The test statistics in data for the rateonly and rate-plus-shape analyses are N obs = 4 ν e candidates and ∆χ 2 = −2.51,obtained from the {E rec ν , θ lep } distribution of the e-like candidates in ν-mode shown in Fig. 28.As shown in Tab.XXIII this test statistics return a p-value of 0.41 (0.46) for the β = 0 hypothesis for the rate-only (rate-plus-shape) analysis, while the pvalue for β = 1 is 0.21 (0.07).This analysis shows that, with the available data, the rate-plus-shape analysis ex- The observed p-values for the β = 0 hypothesis are larger than expectation due to the lower than expected number of observed events in the e-like ν-mode sample.The lower p-value for the rate β = 1 hypothesis is driven by the larger discrepancy between the selected ν e and ν e rate (32 and 4 events respectively) than predicted by the simulated dataset generated, even under the assumption of a maximal CP-violation δ = −π/2 (28.55 and 6.28 events respectively).The p-value in the rate-plus-shape analysis for the β = 1 hypothesis is reduced further, due to the distribution of the ν e -appearance candidates in electron momentum and angle.Their distribution is more compatible with the background expected distribution, as shown in Fig. 19.

XI. JOINT NEUTRINO AND ANTINEUTRINO OSCILLATION ANALYSIS RESULTS
In this section, joint three-flavor oscillation analyses performed with both the frequentist and the Bayesian approaches are presented.The five SK samples introduced in Sec.VI are used, which allows the simultaneous study of the ν e and ν e appearance channels and the ν µ and ν µ disappearance channels.The oscillation parame-ters |∆m 2 32 |, sin 2 θ 23 , sin 2 θ 13 , δ CP , and the mass ordering are determined with and without using the measurement of sin 2 θ 13 from reactor experiments as a constraint.

A. Frequentist analysis
Although two frequentist analyses were introduced in Sec.VIII A, the results are similar, so detailed results are only presented for the analysis which uses {E rec ν , θ lep } templates in this section.Comparisons between both frequentist, and the Bayesian analysis are shown in Sec.XI C.

Results without reactor constraints
This section describes the results obtained by the frequentist analysis when only T2K data are used to estimate the oscillation parameters.The main parameters of interest in this case are δ CP and sin 2 θ 13 that can be directly compared to the reactor measurements.The point estimates for these oscillation parameters and the constant ∆χ 2 = 1 intervals are given for normal and inverted ordering in Tab.XXIV.The ∆χ 2 surfaces are shown in Fig. 36.These intervals have been produced via marginalization of all nuisance and oscillation parameters which are not of interest, as described in Sec.VIII A.
Two-dimensional contours of constant ∆χ 2 for the parameters δ CP and sin 2 θ 13 , along with a comparison with the constraint on sin 2 θ 13 from reactor experiments, are shown in Fig. 37.The point estimate and constant ∆χ 2 confidence intervals of sin 2 θ 13 from T2K data only are slightly larger than what is found by the reactor experiments.However, the T2K-only measurement of sin 2 θ 13 is still consistent with the reactor measurement at the 68% confidence level.
As mentioned above, in this analysis a fifth sample selecting ν e candidates at SK with one delayed Michel electron in the final state has been added for the first time.A comparison of the ∆χ 2 surface for δ CP only including the four single-ring samples used in previous analyses and the results obtained with the inclusion of the fifth sample is shown in Fig. 38.As expected, a small improvement is observed when the new sample is included.Fig. 39 shows a comparison of the constraints in the δ CP -sin 2 θ 13 plane when appearance channels taken in νmode and in ν-mode are considered independently.Both ν-and ν-mode disappearance channels are used in both fits.The ν-and ν-mode datasets alone prefer different values of sin 2 θ 13 , which is driven by the absolute appearance rate.It is clear that the ν-mode appearance sample does not have the power to exclude a zero value of sin 2 θ 13 by itself.In either case, the reactor measurement of sin 2 θ 13 sits near the upper and lower 68% confidence contours for the ν-mode and ν-mode samples respectively.

Results with reactor constraints
Here, the oscillation parameters obtained by the T2K data fit where sin 2 θ 13 is marginalized using the reactor constraint given in Tab.XVIII.The best fit values of the T2K data with the reactor constraint are summarized in Tab.XXV.Fig. 40 shows the 90% constant ∆χ 2 surface in the sin 2 θ 23 -∆m 2  32 plane, assuming normal mass ordering.The interval is compared with other experiments, showing good agreement with IceCube and SK and some tension with MINOS and NOνA.
The NOνA collaboration published ν µ -disappearance results disfavoring maximal mixing for sin 2 θ 23 at 2.6σ [86].The T2K data in the ν µ -and ν µ -disappearance channels, together with the T2K best fit and the expected spectrum produced using the NOνA best fit value for sin 2 θ 23 (higher octant) and ∆m 2  32 , are shown in Fig. 41.Fig. 42 shows the 2D sin 2 θ 13 -δ CP confidence level contours for the data fit including the reactor constraint.The comparisons with the four-sample joint fit are also shown to demonstrate the effect of the inclusion of the CC1π + e-like sample in the appearance analysis.Compared to the best-fit results obtained with the T2K-only data fit in Sec.XI A 1, the inclusion of the CC1π + e-like sample in the data fit with the reactor constraint results in a shift of best-fit value for the δ CP phase towards the maximally violating phase of −π/2.
Since there is a physical boundary at δ CP = ± π 2 , calculating the coverage near the boundary using a Gaussian approximation may not produce accurate results.To solve this problem, the coverage of the 1D ∆χ 2 distribution as a function of δ CP is computed using the Feldman-Cousins approach, described in Sec.VIII.In order to perform the study, 10,000 toy MC experiments were generated for different values of δ CP and the mass ordering.The 1D ∆χ 2 surface obtained with the Feldman-Cousins approach is used to evaluate the 90% confidence intervals for δ CP in both ordering cases, as shown in Fig. 43 A comparison of one-dimensional constant ∆χ 2 contours for normal ordering for δCP using T2K-only data for the four-and five-sample fits.are allowed at 90% confidence for normal (inverted) ordering.
A useful way to visualize the results is to compare the observed number of events in the ν-mode (in both CCQE-like and CC1π + -like samples) and ν-mode e-like samples with the expected events for different values of δ CP , sin 2 θ 23 , and mass ordering.As it is shown in Fig. 44 the T2K data falls outside the physically allowed region.
In order to quantify whether the T2K dataset is consistent with the PMNS framework in terms of significance, an additional toy MC study was performed.An ensemble of 10,000 simulated datasets was obtained in the same way as described in Sec.VIII for the Feldman-Cousins method, with δ CP = −π/2 and normal mass ordering.The values of −2∆ ln L that contain 68.3% and 95.5% of the MC toys were computed and compared to the distribution obtained with the fit of the T2K dataset.As shown in Fig. 45, the T2K data −2∆ ln L distribu- FIG.39.Contours in the sin 2 θ13-δCP plane using T2K-only data, obtained by analysing either the ν-or ν-mode appearance datasets are compared for both orderings.Both ν-and ν-mode disappearance datasets were used in all fits.The yellow band corresponds to the reactor value on sin 2 θ13 from the PDG 2015 [75].40.Allowed region at 90% confidence level for oscillation parameters sin 2 θ23 and ∆m 2 32 using T2K data with the reactor constraint (sin 2 (2θ13) = 0.085 ± 0.005).The normal mass ordering is assumed and the T2K results are compared with NOνA [86], MINOS [87], Super-K [88], and IceCube [89].
tion shows a less extreme fluctuation than at least 5% of the toys MC for all the values of δ CP and mass ordering, i.e. if the experiment is repeated many times and the true value is δ CP = −π/2 with normal ordering, more than 5% of the experiments are expected to show a more extreme statistical fluctuation than the current T2K dataset over the whole range of δ CP and mass ordering.From Fig. 45, the fraction of experiments that would exclude δ CP = 0, π at 90% or 2σ confidence level can be estimated.Assuming a true value of δ CP of -π/2 and normal ordering, 24.3% (21.3%) of toy MC experiments exclude δ CP = 0 (π) at 90% CL.The same can be repeated for different values of δ CP and mass ordering as shown in Tab.XXVI.This section describes the results obtained by the Bayesian analysis when using only T2K data to estimate the parameters sin 2 θ 23 , ∆m 2  32 , sin 2 θ 13 and δ CP with the MCMC method described in Sec.VIII B. In contrast with the frequentist analysis presented in Sec.XI A, the Markov chain walks in a parameter space where the sign of ∆m  41.Comparison of the T2K data in νµ (left) and νµ (right) disappearance channels with the expected spectra obtained with the T2K most probable values of the oscillation parameters and using the NOνA most probable values for sin 2 θ23 (higher octant) and ∆m 2  32 taken from Ref. [86].The ±1σ credible intervals, which have a 68.3% probability of containing the true value, are computed, for each parameter, from the posterior probability density marginalized over all the other parameters as shown in Fig. 46.Fig. 46 also shows the correlations between the oscillation parameters with the map of the marginal posterior density probability and the credible intervals in the space formed by two parameters.
The proportion of the MCMC points with sin 2 θ 23 > 0.5 or < 0.5 gives the posterior probability of the octant.Similarly, the relative proportion of steps with ∆m 2  32 > or < 0 gives the posterior probability of each mass ordering.They are shown in Tab.XXVIII.A Bayes factor can be computed as a ratio of the posterior probabilities [90].The Bayes factor for normal ordering is  e -appearance event rates in the ν-mode samples and in the ν-mode sample as a function of δCP for different values of sin 2 θ23 and both mass orderings, compared to T2K data.The dashed line distinguishes the two solutions for the octant of θ23.

Results with reactor constraints
This section presents the results obtained with the MCMC analysis when adding a Gaussian prior on sin 2 θ 13 with the value given in Tab.XVIII.The posterior mode marginalized over the nuisance parameters is given in Tab.XXIX.Including the reactor prior on sin 2 θ 13 , the best-fit is closer to that obtained by the reactor experiments compared to the T2K-only results.The δ CP best-  fit is closer to the maximum violating value of −π/2 due to the correlations with sin 2 θ 13 shown in Fig. 46.
The MCMC algorithm uses a flat prior on δ CP , but its dependence on this choice of prior has been tested by computing the credible intervals with a flat prior on sin δ CP .The two sets of intervals are in reasonable agreement as shown in Fig. 47.
The Bayes factor for the mass ordering and the θ 23 octant can be computed with the method described in FIG.46.The two-dimensional histograms represent the marginal posterior probability in the two-parameter space as a blue gradient.The white solid (dashed) line is the 90% (1σ) credible interval.The one-dimensional histograms represent the posterior probability density (post.proba.) of the oscillation parameter in the x-axis of the column marginalized over all other parameters.The blue areas are respectively the 1σ (dark), 90% (medium), and 95% (light) credible interval.
Sec. XI B 1. Using the values from Tab. XXX, they are found to be B(NH/IH) = 3.71 and B(sin 2 θ 23 > 0.5/ sin 2 θ 23 < 0.5) = 2.39 respectively.Also in this case, these cannot be considered decisive.

C. Comparison among the oscillation analyses
The frequentist likelihood multiplied by Gaussian penalty terms for the nuisance parameters and uniform priors for the oscillation parameters is equivalent to the Bayesian posterior density, for the same priors.In order to compare the analyses, the posterior probability densities sampled by the Bayesian analyses are converted into a ∆χ 2 function and the intervals are recalculated to extract confidence intervals that are compared with the frequentist analyses.Fig.48 shows the constant ∆χ 2 68% and 90% intervals for all three oscillation analyses in the sin 2 θ 13 -δ CP plane, assuming normal ordering and only using T2K data.Differences exist among the three methods as the 2D templates fitted in the appearance samples are different and the Bayesian analyses does a combined fit of near and far detector samples but no major differences are found between the contours.Fig. 49 shows the constant ∆χ 2 68% and 90% intervals in the sin 2 θ 23 -∆m 2  32 plane for both frequentist and Bayesian fits.Both distributions and intervals agree between fitters, validating the extrapolation of the constraints on the nuisance parameters obtained in the near detector fit to SK.

D. Best fit spectra
An estimate of the oscillation parameters ∆m 2  32 , sin 2 θ 13 and δ CP have been obtained with both frequentist and Bayesian analyses.The agreement between the fit results and the data has been evaluated by compar-   ing the expected spectra after the oscillation fit with the data points as shown in Fig. 50.The best-fit spectra are obtained by sampling 2000 points from the MCMC including the reactor constraint and fitting a Gaussian distribution to calculate the most probable value for the predicted number of events in each energy bin.
In order to extract the ratio of oscillated to unoscillated spectra, the expected spectra are also tuned to the no oscillation case.A coarser binning than the one used in the fit has been used for readability.and ν e and ν e appearance channels by using five samples selected at the far detector, thus, including the new additional CC1π + sample.The data related to this measurement can be found in [91].A comprehensive study has been performed to evaluate the sensitivity of oscillation parameter estimates to neutrino interaction modeling, showing that the impact of these uncertainties is small compared to the total uncertainties on the measurement of all the oscillation parameters.
The general approach followed in this paper, which combines separate analyses of beamline, neutrino interactions, near and far detector selections through sets of systematic parameters and their covariances will be extended with additional data which will be collected by T2K in the coming years in both ν-and ν-modes, and improved near and far detector samples.This is expected to greatly enhance the sensitivity of the T2K experiment to the measurement of the CP-violation phase δ CP as well as more precise measurements of the atmospheric parameters |∆m 2 32 | and sin 2 θ 23 .

FIG. 1 .
FIG. 1.The T2K unoscillated neutrino flux prediction at SK for ν-(left) and ν-(right) modes.The binning used for the flux systematic parameters is also shown.

ν
FIG.4.Multi-nucleon interactions (2p2h) cross section on12 C as a function of energy from the models of Nieves (reference model in the text)[40] and Martini (alternative model in the text)[70].

FIG. 6 .= 1 .N MC, ν 205kA = 1 .
FIG. 6.The event rate at INGRID as a function of time is shown in the top panel.The horizontal dashed redline corresponds to the best fit for the different horn currents with a constant function.The central and bottom panels show the neutrino beam direction measured at INGRID and at MUMON along the horizontal and vertical transverse directions with respect to the beam as a function of time.The dashed vertical lines separate the seven T2K physics runs.

FIG. 9 .
FIG. 9.The SK flux parameters for the ν ( -) µ (left) and ν ( -) e (right) neutrino species in ν-(top) and ν-modes (bottom), shown as a fraction of the nominal value.The bands indicate the 1σ uncertainty on the parameters before (solid, red) and after (hatched, blue) the near detector fit.

FSIFIG. 10 .
FIG.10.Cross-section parameters before (solid, red) and after (hatched, blue) the near detector fit, shown as a fraction of the nominal value (given in Tab.IX).The extent of the colored band shows the 1σ uncertainty.

FIG. 11 .FIG. 12 .
FIG. 11.Distribution of the minimum negative log-likelihood values from fits to the mock datasets (black), with the value from the fit to the data superimposed in red.

2 FIG. 15 .
FIG. 15.Top: data-MC comparison of νµ (right-sign) CC-1-Track (left) and CC-N-Tracks (right) samples after the ND280 fit.The simulation is broken-down by neutrino reaction type.Bottom: same muon momentum distributions for νµ background samples.

FIG. 17 .
FIG.17.∆T0 distribution all FC events observed during T2K Run 1-7 zoomed in on the expected beam arrival time.

FIG. 18 .
FIG.18.The PID likelihood distribution of the observed νmode CC event samples after FCFV and single-ring cuts have been applied.The data are shown as points with statistical error bars and the shaded, stacked histograms are the MC predictions.The expectation is based on the parameters of Tab.XIII.

Fig. 24 FIG. 19 .
FIG.19.The reconstructed energy spectra of the observed νe in ν-mode (left) and νe in ν-mode (right) CC candidate event samples assuming CCQE interaction kinematics.The data are shown as points with statistical error bars and the shaded, stacked histograms are the MC predictions.The expectation is based on the parameters of Tab.XIII.

FIG. 20 .
FIG. 20.Two-dimensional vertex distributions of the observed νe CC candidate events in (X, Y ) and (R 2 , Z) for ν-mode (top) and ν-mode (bottom).The arrow indicates the neutrino beam direction and the dashed line indicates the fiducial volume boundary.Events indicated by open square markers passed all of the νe selection cuts except for the fiducial volume cut.

FIG. 21 .
FIG. 21.The reconstructed energy spectra of the observed νµ in ν-mode (left) and νµ in ν-mode (right) CC candidate event samples assuming CCQE interaction kinematics.The data are shown as points with statistical error bars and the shaded, stacked histograms are the MC predictions, and the rightmost bin includes overflow.The expectation is based on the parameters of Tab.XIII.

1 FIG. 22 . 8 FIG. 23 .
FIG.22.The true momentum distribution for tagged Michel electrons from true CC1π + simulated events reconstructed in the fiducial volume (left) and the tagging efficiency for these events (right).The expectation is based on the parameters of Tab.XIII.

FIG. 24 .
FIG.24.Difference between true and reconstructed energy of the νe CCQE-like (left) and CC1π + (right) simulated samples.The energy for the CCQE-like sample is reconstructed using Eq. 4. For the CC1π + sample the modified interaction kinematics in Eq. 5 are used.The expectation is based on the parameters of Tab.XIII.

FIG. 25 .
FIG. 25.The reconstructed energy spectrum of the observed νe CC1π + candidate events assuming the modified interaction kinematics in equation 5.The data are shown as points with statistical error bars and the shaded, stacked histograms are the MC predictions.The expectation is based on the parameters of Tab.XIII.

FIG. 26 .
FIG. 26.Two-dimensional vertex distributions of the observed νe CC1π + candidate events in (X, Y ) and (R 2 , Z).The arrow indicates the neutrino beam direction and the dashed line indicates the fiducial volume boundary.Events indicated by open square markers passed all of the νe CC1π + selection cuts except for the fiducial volume cut.

FIG. 27 .
FIG. 27.The oscillation probability for νµ → νe and νµ → νe transitions as a function of δCP and MO for Eν = 0.6 GeV and a baseline of 295 km.The other oscillation parameters are fixed to the values shown in Tab.XIII.The dashed line indicates P (νµ → νe) = P (νµ → νe).

FIG. 28 .
FIG. 28.Expected {E rec ν , θ lep } distributions for the signal νµ → νe (top) and wrong-sign background νµ → νe (bottom) in the e-like CCQE-like selection in ν-mode.The superimposed grey dots correspond to the data.The expectation is based on the parameters of Tab.XIII.

2 FIG. 29 .
FIG. 29.{E rec ν ,θ lep } templates for the νe CCQE-like (top), νe CC1π + (middle) in ν-mode and νe CCQE-like (bottom) in ν-mode candidate samples.Both signal and background events are included in the expected distributions based on the oscillation parameters of Tab.XIII.The superimposed grey dots correspond to the data.

FIG. 30 .FIG. 31 .FIG. 32 .
FIG. 30.Ratio of the alternative 2p2h simulated data (left) and alternative 1p1h simulated data (right) to the nominal MC, shown as a function of reconstructed lepton kinematics for the FGD1 CC-0π sample.

FIG. 33 .
FIG.33.The likelihood surfaces (red) from the oscillation fit of the alternative 2p2h simulated data at SK, using the result of the fit to the ND280 simulated data as an input.The oscillation parameters of Tab.XXI (maximal mixing) were used.The likelihood surfaces from the nominal fit are shown by the black line.The contours for sin 2 θ23 (top left), ∆m2  32 (top right), δCP (bottom left) and δCP vs sin 2 θ13 (bottom right) are shown.For the one-dimensional likelihoods the point of minimum ∆χ 2 is given in each legend.

FIG. 34 .
FIG.34.The likelihood surfaces (red) from the oscillation fit of the alternative 2p2h simulated data at SK, using the result of the fit to the ND280 simulated data as an input.The oscillation parameters of Tab.XXI (non-maximal mixing) were used.The likelihood surfaces from the nominal fit are shown by the black line.The contours for sin 2 θ23 (top left), ∆m2  32 (top right), δCP (bottom left) and δCP vs sin 2 θ13 (bottom right) are shown.For the one-dimensional likelihoods the point of minimum ∆χ 2 is given in each legend.

FIG. 35 .
FIG.35.The likelihood surfaces (red) from the oscillation fit of the alternative 1p1h simulated data at SK, using the result of the fit to the ND280 simulated data as an input.The oscillation parameters of Tab.XXI (maximal mixing) were used.The likelihood surfaces from the nominal fit are shown by the black line.The contours for sin 2 θ23 (top left), ∆m2  32 (top right), δCP (bottom left) and δCP vs sin 2 θ13 (bottom right) are shown.For the one-dimensional likelihoods the point of minimum ∆χ 2 is given in each legend.
FIG. 38.A comparison of one-dimensional constant ∆χ 2 contours for normal ordering for δCP using T2K-only data for the four-and five-sample fits.

FIG. 42 .
FIG. 42.A comparison of two-dimensional constant ∆χ 2 contours in the δCP -sin 2 θ13 plane using T2K data with the reactor constraint, for both four-sample (red) and five-sample (black) analyses with normal (left) and inverted (right) mass ordering hypotheses.The contours are produced by marginalizing the likelihood with respect to all parameters other than the parameters of interest.

FIG. 43 .
FIG.43.One dimensional ∆χ 2 surfaces for oscillation parameter δCP using T2K data with the reactor constraint.The critical ∆χ 2 values obtained with the Feldman-Cousins method are used to evaluate the 90% confidence level with the proper coverage.

FIG. 44 .
FIG. 44.Total predicted ν( -) e -appearance event rates in the ν-mode samples and in the ν-mode sample as a function of δCP for different values of sin 2 θ23 and both mass orderings, compared to T2K data.The dashed line distinguishes the two solutions for the octant of θ23.

FIG. 45 .
FIG.45.One-dimensional marginal ∆χ 2 surfaces for oscillation parameters δCP and sin 2 θ13 using T2K data with the reactor constraint.The contour is produced by marginalizing the likelihood with respect to all parameters other than the parameter of interest.The red line shows the critical ∆χ 2 values obtained with the Feldman-Cousins method, used to evaluate the 90% confidence level with the proper coverage.The green line show the ∆χ 2 obtained with the fit to the T2K data.

32 2 . 4 FIG. 47 .
FIG.47.δCP marginal posterior probability as obtained with the MCMC method.The credible intervals for ±1σ, ±90% and ±95% are shown when using a flat prior in δCP (usual fit), and when converting to a flat prior in sin δCP.

15 FIG. 50 .
FIG.50.Comparison of the best-fit oscillated MC energy spectra, unoscillated spectra and T2K data for the five samples used in the fit: µ-like sample in ν-mode and ν-mode (top left and right), single ring e-like appearance sample in ν-mode and ν-mode (middle left and right), CC1π + e-like appearance sample in ν-mode (bottom).The larger unoscillated spectra in the CC1π + e-like sample compared to the single ring sample is due to the relatively large background of νµ in the CC1π + sample, which does not disappear in the no-oscillation case.The ratio of the best fit to unoscillated spectra are also shown.

TABLE I .
T2K data-taking periods and collected POT used in the analyses presented in this paper.

TABLE II .
Weights applied to coherent pion interactions as a function of the pion energy, Eπ.

TABLE VI .
Systematic uncertainty on the total event rate affecting the near detector samples.

TABLE VII .
Prefit and postfit values for the SK ν-mode flux parameters

TABLE IX .
Prefit and postfit values for the cross-section parameters used in the oscillation fits.If no prefit uncertainty is shown then the parameter had a flat prior assigned.If a parameter was not constrained by the ND280 fit this is noted in the postfit column.

TABLE XI .
Prefit and postfit expected fraction of events for each neutrino interactions for samples taken in ν-mode.∆T0 distribution of all FC, OD and LE events within ±500 µs of the expected beam arrival time observed during T2K Run 1-7.The histograms are stacked in that order.

TABLE XIII .
[27]es of the oscillation parameters used for the Monte Carlo simulation at SK.The values of sin 2 θ12, ∆m221 and sin 2 θ13 are taken from Ref.[75], while all the other oscillation parameters correspond to the most probable values obtained by the Bayesian analysis in Ref.[27].

TABLE XIV .
Event reduction for the νe CC selection at the far detector.The numbers of expected MC events divided into five categories are shown after each selection criterion is applied.The MC expectation is based upon three-neutrino oscillations with the parameters as shown in Tab.XIII.The visible energy, E vis , is greater than 100 MeV d There is no reconstructed Michel electron e The reconstructed energy, E rec ν , is less than 1.25 GeV f The event is not consistent with a π 0 hypothesis neutrino energy for the final CC1π a There is only one reconstructed Cherenkov ring b The ring is e-like c

TABLE XV .
Event reduction for the νµ CC selection at the far detector.The numbers of expected MC events divided into four categories are shown after each selection criterion is applied.The MC expectation is based upon three-neutrino oscillations with the parameters as shown in Tab.XIII.
g There is only one reconstructed Cherenkov ring h The ring is µ-like i The reconstructed momentum, pµ, is greater than 200 MeV/c j There are less than two reconstructed Michel electrons

TABLE XVI .
Event reduction for the νe CC1π + selection at the far detector.The numbers of expected MC events divided into five categories are shown after each selection criterion is applied.The MC expectation is based on the parameters of Tab.XIII.There is only one reconstructed Cherenkov ring b The ring is e-like c The visible energy, E vis , is greater than 100 MeV d There is one reconstructed Michel electron e The reconstructed energy, E rec ν , is less than 1.25 GeV f The event is not consistent with a π 0 hypothesis a 12 s 13 s 23 (c 12 c 23 cos δ CP − s 12 s 13 s 23 ) cos φ 23 sin φ 31 sin φ 21 12 c 23 s 12 s 13 s 23 sin δ CP sin φ 32 sin φ 31 sin φ 21 + 4s 2 2c 12 c 23 s 12 s 23 s 13 cos δ CP sin 2 φ 21

TABLE XIX .
Effect of 1σ variation of the systematic uncertainties on the predicted event rates of the ν-mode samples.

TABLE XX .
Effect of 1σ variation of the systematic uncertainties on the predicted event rates of the ν-mode samples.

TABLE XXI .
Oscillation parameter values used as inputs for the studies of simulated data.

TABLE XXII .
Summary of the maximum bias seen in oscillation parameters when fitting simulated datasets (defined in Sect.III B), presented as a fraction of the expected 1σ uncertainty on each parameter.Fits were performed for a number of true oscillation parameter assumptions.The numbers shown here are the maximum bias found for a given parameter and alternative model across all true oscillation parameter values tested.

TABLE XXV .
Best-fit results and the 1σ confidence interval of the T2K data fit with the reactor constraint with normal and inverted hypotheses.

TABLE XXVI .
The fraction of toy experiments for which δCP = 0, π and normal and inverted ordering are excluded at 90% and 2σ confidence is shown for different true values of δCP and mass ordering.10,000 toy experiments are used for each set of values.
2  32can flip, and results are presented for both mass orderings.The best-fit point and ±1σ credible interval for each parameter, obtained with the KDE method, are

TABLE XXVIII .
Posterior probabilities for the mass orderings and sin 2 θ23 when fitting T2K data only with an MCMC method.sin 2 θ23 < 0.5 sin 2 θ23 > 0.5 Line Total

TABLE XXIX .
Best-fit results and the 1σ credible interval of the T2K data fit with the reactor constraint with the MCMC analyses including both mass orderings.

TABLE XXX .
Posterior probabilities for the mass orderings and sin 2 θ23 with an MCMC method when fitting T2K data with the reactor constraint.sin 2 θ23 < 0.5 sin 2 θ23 > 0.5 Line Total