Measurements of neutrino oscillation parameters from the T2K experiment using \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$3.6\times 10^{21}$$\end{document}3.6×1021 protons on target

The T2K experiment presents new measurements of neutrino oscillation parameters using \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$19.7(16.3)\times 10^{20}$$\end{document}19.7(16.3)×1020 protons on target (POT) in (anti-)neutrino mode at the far detector (FD). Compared to the previous analysis, an additional \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$4.7\times 10^{20}$$\end{document}4.7×1020 POT neutrino data was collected at the FD. Significant improvements were made to the analysis methodology, with the near-detector analysis introducing new selections and using more than double the data. Additionally, this is the first T2K oscillation analysis to use NA61/SHINE data on a replica of the T2K target to tune the neutrino flux model, and the neutrino interaction model was improved to include new nuclear effects and calculations. Frequentist and Bayesian analyses are presented, including results on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sin ^2\theta _{13}$$\end{document}sin2θ13 and the impact of priors on the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta _{\textrm{CP}}$$\end{document}δCP measurement. Both analyses prefer the normal mass ordering and upper octant of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sin ^2\theta _{23}$$\end{document}sin2θ23 with a nearly maximally CP-violating phase. Assuming the normal ordering and using the constraint on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sin ^2\theta _{13}$$\end{document}sin2θ13 from reactors, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sin ^2\theta _{23}=0.561^{+0.021}_{-0.032}$$\end{document}sin2θ23=0.561-0.032+0.021 using Feldman–Cousins corrected intervals, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varDelta {}m^2_{32}=2.494_{-0.058}^{+0.041}\times 10^{-3}~\text {eV}^2$$\end{document}Δm322=2.494-0.058+0.041×10-3eV2 using constant \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varDelta \chi ^{2}$$\end{document}Δχ2 intervals. The CP-violating phase is constrained to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta _{\textrm{CP}}=-1.97_{-0.70}^{+0.97}$$\end{document}δCP=-1.97-0.70+0.97 using Feldman–Cousins corrected intervals, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta _{\textrm{CP}}=0,\pi $$\end{document}δCP=0,π is excluded at more than 90% confidence level. A Jarlskog invariant of zero is excluded at more than \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2\sigma $$\end{document}2σ credible level using a flat prior in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta _{\textrm{CP}},$$\end{document}δCP, and just below \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2\sigma $$\end{document}2σ using a flat prior in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sin \delta _{\textrm{CP}}.$$\end{document}sinδCP. When the external constraint on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sin ^2\theta _{13}$$\end{document}sin2θ13 is removed, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sin ^2\theta _{13}=28.0^{+2.8}_{-6.5}\times 10^{-3},$$\end{document}sin2θ13=28.0-6.5+2.8×10-3, in agreement with measurements from reactor experiments. These results are consistent with previous T2K analyses.


Introduction
The Tokai to Kamioka (T2K) experiment produces a beam of predominantly muon neutrinos by impinging protons from an accelerator onto a target, using magnetic horns to direct the outgoing collision products which thereafter decay into the neutrinos that form the beam.A suite of near detectors, 280 m downstream of the production target, characterise the neutrinos before long-baseline oscillations take effect, and a far detector, 295 km away, measures the long-baseline oscillations.This paper first introduces the neutrino oscillation formalism in Sect. 1 and summarises the T2K experiment in Sect. 2. Section 3 outlines the updates to the previous analysis [1,2], with the systematic uncertainties presented in detail in Sect. 4 for the neutrino flux, and in Sect. 5 for the neutrino interaction model.The analysis of near-detector data, which constrains the majority of the systematic uncertainties in the oscillation analysis, is described in Sect.6.The far-detector selections are described in Sect.7, and the new constraints on the oscillation parameters are presented in Sect.8. Section 9 summarises the simulated data studies, which act to increase the uncertainty on the oscillation parameters by studying the impact of alternative interaction models.The results are summarised in Sect.10, and the data release, amongst other supplementary material, is provided in the appendices.
The observation of neutrino survival probabilities changing as a function of both flavour and distance travelled was established in the late 1990s by Super-Kamiokande (SK) [3].Their measurements of neutrinos produced by cosmic rays in the atmosphere found that muon neutrinos disappeared after travelling through the Earth, whereas electron neutrinos did not.A few years later, the Sudbury Neutrino Observatory (SNO) found evidence that neutrino flavour change was responsible for the measured deficit of electron neutrinos compared to what was predicted from the Sun [4].Neutrino flavour changing was also confirmed using artificial sources of neutrinos in the long-baseline reactor experiment KamLAND [5] which measured the disappearance of ν e , and accelerator experiments K2K [6] and MINOS [7] which measured the disappearance of ν μ and ν μ .These experiments additionally characterised the oscillation curve in the ratio of the distance travelled over the neutrino energy, L/E, which governs the oscillation probability.The results can be summarised in a framework with three active neutrinos, where at least two neutrinos have non-zero mass.The flavour and mass eigenstates of the neutrinos, |ν l and |ν i respectively, are separate and can be related by a 3 × 3 unitary mixing matrix U as |ν l = U |ν i .The mixing matrix is the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix, which can be parametrised by three mixing angles, θ 12 , θ 23 , θ 13 , and a CP-violating phase, δ CP [8,9].The probabilities for neutrino flavour oscillations can then be expressed as functions of these mixing angles and the mass-squared differences, Δm 2 i j = m 2 i − m 2 j where m i is the mass of the ith neutrino mass eigenstate.The m 2 > m 1 ordering was established by measurements of solar neutrinos across multiple experiments [10].The ordering of the remaining mass states is unknown, with m 3 > m 2 > m 1 referred to as the normal ordering (NO), and m 2 > m 1 > m 3 as the inverted ordering (IO).This analysis uses the Particle Data Group (PDG) [11] convention for the order of the mixing matrices, U = U 23 ⊗ U 13 ⊗ U 12 .
The results from SK, SNO, and KamLAND showed that both θ 23 and θ 12 were non-zero.The last mixing angle, θ 13 , was indicated to be non-zero through T2K's 2.5σ measurement of ν μ → ν e [12].It was later precisely measured by short-baseline experiments Daya Bay [13], RENO [14], and Double Chooz [15], observing the disappearance of ν e from nuclear reactors.The long-baseline accelerator experiments T2K and NOvA subsequently observed the appearance of ν e in a ν μ beam at high significance [16,17], and NOvA observed ν e appearance in a ν μ beam at 4.4σ [18].The nonzero θ 13 mixing angle implies that a measurement of δ CP is possible in long-baseline accelerator-based experiments by measuring the appearance of ν e and ν e in ν μ and ν μ beams, respectively.
On its way to the FD, the beam passes through matter and the presence of electrons modifies the oscillation probabilities as compared to those in vacuum.Namely, chargedcurrent elastic scattering on electrons is possible for ν e and ν e (hereafter referred to as ( ν ) e ), but not for the other flavours [19,20].The sign of the matter effect differs for neutrinos and anti-neutrinos, and the magnitude is a function of the density of electrons in the path of the neutrinos, n e , the weak interaction coupling strength, G F , and the neutrino energy, E. The matter effect in the sun was central to measuring m 2 > m 1 .The probability for ( ν )  e appearance as a function of neutrino energy, E, and baseline, L , including a first-order approximation of the matter effects, is [ The first term in Eq. 1 is proportional to sin 2 θ 23 , which renders the ( ν ) e appearance sensitive to whether θ 23 is above or below π/4, referred to as the octant of θ 23 .This in turn determines whether the ν 3 mass eigenstate has a larger admixture of ν μ or ν τ .The term containing sin δ CP in Eq. 1 has the opposite sign for neutrinos and anti-neutrinos, and allows for CP symmetry violation if δ CP is different from 0 or π.The term containing cos δ CP does not violate CP symmetry, but can change the shape of the ( ν )  e appearance energy spectrum, and is important for precisely measuring δ CP .In T2K, the term proportional to sin δ CP can change the appearance probability by as much as ±30% given the current knowledge of the other mixing angles.J = J 0 sin δ CP is referred to as the Jarlskog invariant [22,23] and is a basis-independent measure of the CP-violation.This analysis presents T2K's constraints on Δm 2  32 , sin 2 θ 23 , sin 2 θ 13 , δ CP , J, and the mass ordering.

The T2K experiment
To measure δ CP and the other oscillation parameters, T2K uses a beamline that produces predominantly muon-flavoured neutrinos or anti-neutrinos with a peak energy of E ν ≈ 0.6 GeV, and has been alternating between neutrino and antineutrino configurations since 2014.A suite of near detectors (NDs), approximately 280 m from the beam production target, characterise T2K's neutrino beam before long-baseline oscillations become likely.The far detector (FD) is 295 km away and measures the appearance of ( ν ) e and the disappearance of ( ν )  μ in the ( ν ) μ -dominated beam.The rate and directional stability of the neutrino beam are measured by the onaxis neutrino ND, INGRID.The second ND, ND280, and the FD, SK, are 2.5 • off-axis with respect to the upstream proton beam that impinges on the neutrino production target.By being placed off-axis, the detectors sample a narrower neutrino energy distribution, peaking near the maximum of the ( ν )  e appearance spectrum.

Beamline
The T2K neutrino beam is produced at the Japan Proton Accelerator Research Complex (J-PARC) in Tokai, Ibaraki, by a high-intensity proton beam, incident on a production target [24].At J-PARC, H − ions from an ion source are accelerated to an energy of 400 MeV in a linear accelerator.Charge-stripping foils convert the beam to H + at injection into the rapid-cycling synchrotron, which accelerates the proton beam to 3 GeV.These protons are then injected into the main ring (MR) synchrotron, where they are accelerated to 30 GeV.The proton beam from the MR consists of eight bunches with width ∼ 80 ns (3σ ), referred to as a "spill", produced every 2.48 s.
The protons are extracted from the MR to the neutrino beamline, which consists of a series of normal-and superconducting magnets that are used to bend the proton beam in the direction of the FD, and to focus the beam onto the neutrino production target.The proton beam power, as well as the position, angle, and size of the proton beam at the target, are precisely measured by a series of proton beam monitors [24,25] installed along the neutrino beamline.
The 30 GeV protons strike a 91.4 cm-long monolithic graphite target installed in the first of three electromagnetic focusing horns.Outgoing charged pions and kaons are focused by these horns, which have been operating at a current of ±250 kA for nearly the full T2K run to date.The polarity of the horns can be set to focus either positively or negatively charged outgoing particles, and a 96 m-long decay volume is located directly downstream of the focusing system.Positively charged pions decay into positively charged muons and muon neutrinos, whilst negatively charged pions decay into negatively charged muons and muon anti-neutrinos.The former is referred to as ν-mode and the latter as ν-mode.Kaon and muon decays are the primary contributors to the ν e contamination in the ν μ -dominated beam.
A beam dump is situated at the end of the decay volume and absorbs surviving hadrons.A muon monitor downstream of the beam dump, MUMON [26], measures the intensity and profile of muons that have more than 5 GeV of energy.This measurement is used as a proxy for stability of the associated neutrino beam.The predicted neutrino fluxes and uncertainties are described in detail in Sect. 4.
The MR proton beam power has reached a maximum of 515 kW, and the protons on target (POT) and power history are shown in Fig. 1.Scheduled upgrades will increase the beam power to 1.3 MW and operate the focusing horns at ±320 kA current.This will significantly increase the POT per run cycle and provide more neutrinos at the ND and FD per POT.It will also reduce the ν μ and ν μ backgrounds in ν-mode and ν-mode [27,28], respectively, referred to as the wrong-sign component.

Near detectors
Two NDs are used directly in the oscillation analysis: the on-axis INGRID, and the off-axis ND280.Both detectors are housed in the same pit underground, with the centres of ND280 and INGRID approximately 24 m and 33 m, respectively, below the surface.
INGRID [29] is designed to measure the profile and stability of the neutrino beam.It samples the beam spill-by-spill with a transverse cross section of 10 × 10 m 2 with 14 identical modules arranged as a cross, as shown in Fig. 2. Each of the modules alternates iron target plates of 6.5 cm thickness with tracking scintillator planes of 1 cm thickness, for a total of 9 iron plates and 11 scintillator planes, and is surrounded by scintillator planes acting as vetoes.A module exposes a 1.24 × 1.24 m 2 area facing the beam, and provides a 7.1 t target mass.INGRID measures the beam direction with an accuracy higher than 0.4 mrad, within the required precision of ±1 mrad for the oscillation analysis.
ND280, hereafter referred to as the ND, is used to constrain the uncertainties on the neutrino flux and interactions in the analysis.It is a magnetised detector consisting of dif-Fig. 1 The protons on target (POT) delivered to T2K by the MR over time, with the beam intensity overlaid.The ND280 analysis uses runs 2 to 9, and the INGRID and FD analyses use runs 1 to 10, with run-by-run POT listed in Table 1 Fig. 2 The INGRID on-axis ND, used to measure the neutrino beam profile and rate [29].The beam direction is shown as into the paper ferent sub-detectors as shown in Fig. 3.The ND measures 5.6 m × 6.1 m × 7.6 m (width × height × length) around its outer edges including the magnet with the coordinate convention being z pointing along the nominal neutrino beam axis, Fig. 3 The ND280 off-axis ND, used to measure the neutrino flux and interactions before long-baseline oscillations [24].The detector coordinates and beam direction are superimposed, with the sub-detectors are labelled accordingly with x and y being the horizontal and vertical directions, respectively.The refurbished magnet from the UA1 [30,31] and NOMAD [32] experiments at CERN provides a magnetic field of 0.2 T, and the magnet yoke is instrumented with layers of plastic scintillator called the Side Muon Range Detector (SMRD) [33].Inside the magnet enclosure there is an electromagnetic calorimeter (ECal) [34] surrounding the inner detector, which is used to distinguish track-like and showerlike objects, and is made of alternating layers of plastic scintillator and lead.
The inner detector region houses the π 0 detector (P∅D) [35] in the upstream portion, which is made of alternating layers of water bags, brass sheets, and triangular x−y scintillator planes.The water bags can be filled with either water or air.The P∅D has its own ECal modules upstream and downstream of the water target region, made from alternating scintillator planes and lead sheets.The P∅D, ECal and SMRD also act as vetoes for interactions originating outside the detector, e.g.cosmic-ray muons and neutrino interactions in the sand upstream of the detector hall.Downstream in the direction of the FD, there are two Fine-Grained Detectors (FGDs) [36], which are each sandwiched by Time Projection Chambers (TPCs) [37].These sub-detectors are together referred to as the "tracker".The most upstream FGD (FGD1) is made of 15 polystyrene scintillator modules.One module is 186.4 cm × 186.4 cm × 2.02 cm and consists of two scintillator layers oriented in x and y, with each layer containing 192 9.6 mm wide square bars approximately 2 m long, which are read out at one end.The second FGD (FGD2) contains six passive water modules, each sandwiched by polystyrene scintillator modules identical to those in FGD1.The TPCs use a Ar:CF 4 :iC 4 H 10 gas mixture in a 95:3:2 concentration, and have a space point resolution of approximately 1 mm.This analysis selects interactions occurring in either FGD, using the FGDs and TPCs for track reconstruction and particle identification.The selection is detailed in Sect.6.1.The FGDs are capable of tracking charged particles, performing particle identification, and calculating momentum-by-range for contained particles.The TPCs are three-dimensional trackers which measure momentum through the curvature of the tracks in the magnetic field, with a resolution of δp ⊥ p ⊥ ∼ 0.1 p ⊥ , where p ⊥ is the momentum perpendicular to the magnetic field.The TPCs also provide excellent particle identification.

Far detector
The Super-Kamiokande (SK) detector [24,38] is the far detector (FD) for T2K.SK is a large water Cherenkov detector located 295.3 km from the neutrino production target with a 2.7 km water-equivalent overburden.It is filled with 50 kt of ultrapure water that is optically separated into an inner detector, ID, which forms the primary target for neutrino interactions, and an outer detector, OD, which serves to veto external backgrounds.
The ID is instrumented with 11,129 inward-facing photomultiplier tubes (PMTs) with 20-inch diameter, providing a total photocathode coverage of 40%.The OD is instrumented with 1885 8-inch outward-facing PMTs, which are connected to wavelength shifting plates and are attached to the same stainless steel structure that houses the ID PMTs.
The structure is offset 2 m from the wall of the OD and there is a 55 cm dead region between the ID and OD surfaces.Charged particles are detected by their Cherenkov ring pattern, and events are classified by the number of primary rings, the ring pattern of each ring, and the number of timedelayed electron rings consistent with a muon decay, hereafter referred to as "Michel electrons".This analysis selects single-ring ("1R") events, where the ring is either electronlike (1Re) or muon-like (1Rμ), with a selection-dependent cut on the number of delayed Michel electrons ("de").The FD selections are detailed in Sect.7.
The data used in this analysis were taken over two different periods of the SK detector operations and span the years 2010-2020, during what is referred to as the SK-IV period.Of the 36.01 × 10 20 POT reported here, 31.29 × 10 20 (runs 1-9) were collected in 2010-2018.In June 2018, SK detector operations were stopped for refurbishment in preparation for the gadolinium (Gd) loading of the water target for the SK-Gd project [39,40].During this work the detector surfaces were cleaned to remove rust and other impurities, detector walls were repaired to fix minor leaks, and failed PMTs were replaced in the ID and OD.This SK detector period is referred to as SK-V.  2POT available for analysis in the ν μ and ν μ modes, respectively.For a detailed breakdown of the POT in each run period, consult Table 1.Gadolinium loading commenced in July 2020, and this analysis does not include such data.

Updates from previous analysis
This section provides an overview of the improvements to T2K's previously published oscillation analysis [1,2], which are detailed in the subsequent sections.
• Data at the FD: The data at INGRID and the FD increased by 4.73×10 20 POT (+33%) in ν-mode, increasing the overall amount by 15%, detailed in Sect.7. • Data at the ND: The data at the ND increased by 5.73 × 10 20 POT (+99%) in ν-mode, and by 4.48 × 10 20 POT (+116%) in ν-mode, increasing the overall amount by 106%, detailed in Sect.6. • Selections at the ND: The increased data allowed for refining the ν-mode selections and re-binning all existing selections, improving the constraints on the systematic uncertainties from the ND in the oscillation analysis, detailed in Sect.6.1.• FD reprocessing: An updated model for the dark rate and gain drift in the PMTs had a slight impact on the reconstruction and the number of observed data events.The processing introduced one more ν-mode electronlike event, and three fewer ν-mode muon-like events, and had no overall effect on the ν-mode samples, detailed in Sect.7.
• Neutrino flux model: The neutrino flux was constrained using charged pion production data on a replica of the T2K production target from NA61/SHINE [41].Data on a thin target [42] was also used when appropriate.This reduced the flux uncertainties before the ND analysis from ∼ 9% down to ∼ 5% in the neutrino flux peak, detailed in Sect. 4.
• Neutrino interaction model: Several changes to the neutrino interaction model were made.The largest changes were switching to a more sophisticated spectral-function based nuclear model [43] for charged-current quasielastic (CCQE) interactions, introducing an additional uncertainty due to nuclear effects in the four-momentum transferred to the nucleus (Q 2 ), and adding an uncertainty for the nucleon removal energy.The nuclearcascade model for pions was tuned to external data [44], and the FD parametrisation was constrained by the fit to ND data, whereas it was previously allowed to vary separately.The interaction model for pions re-scattering within the detector at the ND and FD were unified, and is identical to the pion final-state interaction model, detailed in Sect. 5.However, constraints of e-scattering within the ND were not propagated to re-scattering at the FD, as the uncertainties were kept uncorrelated.

Neutrino flux model
This is the first T2K oscillation analysis to use hadron production measurements made on a replica of the T2K target by the NA61/SHINE experiment at CERN [41].The method for predicting the neutrino flux and propagating the associated uncertainties remains the same as in previous results [1,2,45].FLUKA 2011.2x[46,47] is used to simulate interactions inside the target.The outgoing particles from the target, which later decay to neutrinos, are tracked through the horn field using the GEANT3-based JNUBEAM package [48].
The prediction for pions exiting the target's surface are tuned to π + and π − yields measured by the NA61/SHINE experiment, using data collected in 2009 with a replica of the T2K production target [41].Pions that leave the target and are within the phase space covered by the replica target data, which is about 90% of the neutrinos at the flux peak, are given a weight where dn is the POT-normalised differential yield for data ("NA61") and simulation Monte-Carlo ("MC"), with exiting momentum p, polar angle θ, and longitudinal position z along the target for an exiting particle of type i = {π + , π − }.
For the particles leaving the target, no additional tuning weight is applied for any of the interactions or trajectories inside the target.Simulations for particles that are not covered by the replica target data, and interactions occurring outside the target, are tuned to NA61/SHINE data on π ± , K ± , K 0 s , , and p yields from a thin target taken in 2009 [42], and other external measurements, applying the same method as previous T2K analyses [45].The percentage of hadronic interactions which are tuned by external data is shown in Table 2.
In the previous thin-target tuning, a large uncertainty on the cross section of proton production was assigned.In the replica-target based tuning, this uncertainty is no longer necessary for particles covered by the replica target data, because the exiting particle yields can be tuned directly without referring to the interaction history inside the target.The uncertainties from NA61/SHINE are then incorporated with the uncertainties associated with the proton beam profile and out-of-target interactions to give the total uncertainty.
For the unconstrained interactions not covered by thinor replica-target data, a systematic uncertainty is calculated by dividing the kinematic phase space parametrised by Feynman-x F and transverse momentum, p T , into six regions.A 50% fully correlated normalisation uncertainty and a 50% shape uncertainty uncorrelated between the regions is assigned.The size of the uncertainty is motivated by comparing the hadron interaction models in FLUKA 2011.2c[46,47] and the GEANT 4.10.03[49] FTFP_BERT and FTF_BIC physics lists.
The predicted flux distributions are provided in Ref. [50] and are shown for the FD in Fig. 4. The largest difference compared to the previous neutrino flux prediction is the reduction of the ν μ component in ν-mode, and the ν μ component in ν-mode ("right-sign"), by 5-10% around the flux peak.Due to the large uncertainty on the hadron interactions in the previous tuning, this difference was covered by the flux uncertainties.To more clearly see wrong-sign and background contributions, the predicted neutrino flux spectra are also shown in logarithmic scale and for a wider range of energies in Fig. 5.
Overall, tuning with the NA61/SHINE 2009 replica target data reduces the uncertainty from 9 to 5% near the flux peak, as shown in Fig. 6.In future T2K analyses, outgoing kaons will also be tuned using NA61/SHINE T2K replica target data from 2010, published in 2019 [51].This will reduce the flux uncertainty at higher energies to ∼ 5%.With a reduced uncertainty contribution from hadron production errors, uncertainties coming from other sources are now dominant in some energy regions.In particular, uncertainties on the proton beam profile and neutrino beam off-axis angle significantly contribute to the uncertainty on the high-energy edge of the flux peak, since the width of the energy spectrum is directly affected by shifts in the off-axis angle.

Neutrino interaction model
Measurements of neutrino oscillations at T2K rely on comparing the neutrino interaction rates at the ND and the FD as a function of the incoming neutrino energy and flavour.These are determined from the observed products of neutri-Fig. 4 The predicted unoscillated neutrino fluxes at the FD in ν-mode (top) and ν-mode (bottom).The ν e and ν e components are scaled by ×100.The solid lines show the predictions after tuning to NA61/SHINE data on the T2K replica target, and the dotted grey lines show the predictions in the previous T2K analysis [1,2], tuned to thin target hadron production data.The bottom inset shows the ratio of the flux from the replica target tuning to the flux from the thin target tuning nos interacting with the nuclei inside the detectors, which requires a model to translate what is observed in the detector to information about the neutrino that interacted.Neutrino interaction uncertainties impact the oscillation analysis by changing the expected rate of neutrino interactions, altering the accuracy of the neutrino energy reconstruction, and complicating the extrapolation of model constraints from the ND to the FD.More details can be found in Refs.[1,[52][53][54].The neutrino interaction model has been significantly improved since the last analysis [1].This section first provides an overview of the components of the model and then discusses the associated uncertainties and their parametrisations.As briefly mentioned in Sect. 2 and detailed further in Sects.6.1 and 7, this analysis selects charged-current (CC) neutrino interaction events and has no dedicated neutralcurrent (NC) selections.The oscillation analysis at the FD specifically selects single-ring events and the model focuses on the treatment of such interactions.In these interactions, CCQE and 2p2h are the main contributors and are discussed next.Neutrino interactions in which a single pion is produced and the pion is missed -either due to its kinematics or by it Fig. 6 Uncertainty on the right-sign flux in ν-mode (top) and right-(middle) and wrong-sign (bottom) fluxes in ν-mode, broken down by the sources of uncertainty.The solid black line shows the total flux uncertainty in this analysis, and the dashed black line shows the total uncertainty for the previous T2K analysis [1,2], which used NA61/SHINE thin target data.The grey shaded region shows the shape of the neutrino flux being absorbed in the nuclear medium -are also an important contributor.Fig. 7 Neutrino cross sections for muon neutrinos interacting on a water target in NEUT, broken down by interaction mode as a function of neutrino energy.The predictions have been modified from their default to reflect the input model used in the oscillation analysis.The surviving muon neutrino flux as seen by the FD is shown with a white line, and the unoscillated muon neutrino flux as seen by the ND is shown as the grey shaded region.The figure is adapted from Ref. [55] 5.1 Base interaction model Simulations of neutrino interactions are performed with version 5.4.0 of the NEUT neutrino-nucleus interaction event generator [55][56][57].NEUT takes inputs from a variety of theoretical models for separate neutrino interaction channels.The total cross sections for each channel as a function of neutrino energy, overlaid on the T2K oscillated and unoscillated muon neutrino fluxes, are shown in Fig. 7.An overview of the channels most relevant to this analysis is presented below.

1p1h
One-particle one-hole (1p1h) interactions describe chargedcurrent quasi-elastic (CCQE) and neutral-current elastic (NCE) neutrino interactions in which a single nucleon from inside a target nucleus is ejected.CCQE interactions, which usually produce single-ring electron-like or muon-like events, are the dominant contributor to the FD event samples, making up roughly 70% of the 1Rμ selection.In NEUT, 1p1h interactions are modelled according to the scheme presented in Refs.[43,55], sometimes referred to as the "Benhar Spectral Function" model.This approach relies on the plane wave impulse approximation to factorise the 1p1h cross-section calculation into an expression containing a single-nucleon factor alongside a spectral function (SF).The SF is a twodimensional distribution describing the probability of finding a nucleon with momentum, |p|, and removal energy, E rmv , which corresponds to the energy required to remove the nucleon from the nuclear potential.This formalism provides a realistic description of the nuclear ground state and is built largely from exclusive measurements of 1p1h interactions in electron scattering, with additional theory-based contributions to describe the role of initial-state correlations between neighbouring nucleons.As an example, the two-dimensional SF for oxygen is shown in Fig. 8, which exhibits the shell structure of the nucleus.
The single-nucleon component of the 1p1h cross section uses the BBBA05 [59] description for the vector part of the nucleon form factors, and a simple dipole form for the axial part.The nucleon axial mass parameter appearing in the form factor, M Q E A , is constrained using bubble chamber measurements of neutrino interactions on light nuclear targets, as detailed later in Sect.5.2.

2p2h
In two-particle two-hole (2p2h) interactions, a neutrino interacts with a correlated pair of nucleons, ejecting both from the nucleus.Although this is not a dominant process at T2K, it usually produces single-ring electron-like or muon-like events in the FD -making up about 12% of the 1Rμ selection at the FD -and is therefore important to the oscillation analysis.As T2K's neutrino energy estimator is based on the assumption that the interaction was CCQE, applying it to 2p2h events causes a natural bias.Thus it is crucial that the relative contribution of 2p2h events to the selections, and the bias they cause to the neutrino energy estimator, are well modelled.NEUT describes the charged-current 2p2h cross section and outgoing lepton kinematics with the Nieves et al. model [60].In this model, the 2p2h cross section peaks in two distinct regions of momentum and energy transfer, referred to as "Δ" and "non-Δ" excitation regions, which each cause distinctly different biases in neutrino energy reconstruction [1].Neutral-current 2p2h interactions are not simulated in NEUT.Their inclusion would have a negligible impact on the oscillation analysis as such interactions would make a small contribution to an already small NC background, which is prescribed large uncertainties.

Single-pion production
Single-pion production (SPP) processes are the dominant contributor for the T2K FD sample that requires a single electron-like ring with one delayed decay electron (referred to as 1Re1de in Sect.7).The events also contribute to the other event samples when the pion is not observed due to interactions in the detector or the nucleus, or due to reconstruction inefficiencies.SPP at T2K stems mostly from the neutrino-induced excitation of an initial-state nucleon to a Fig. 8 The two-dimensional probability density distribution for the spectral function for oxygen in NEUT [43] (left), and the projection onto the removal energy axis (right).On the left, the darker colour represents a higher probability of finding an initial-state nucleon with a particular removal energy and momentum.The two sharp p-shells at ∼ 12 MeV and E rmv ∼ 18 MeV, and the larger diffuse s-shell at E rmv ∼ 20−65 MeV and |p| < 100 MeV/c, are visible.The predictions for the shell positions from another model [58] are overlaid on the right with dashed lines, for protons (red) and neutrons (blue).The energy in MeV is labelled for each prediction baryon resonance that decays into a pion and a nucleon, and makes up about 13% of the 1Rμ selection.These processes are described in NEUT by the Rein-Sehgal (RS) model [61] in the outgoing hadronic mass region W < 2.0 GeV, with additional improvements to the nucleon axial form factors [62,63] and the inclusion of the final-state lepton mass in the calculation [64][65][66].Whilst Δ(1232) excitations are the dominant contributors to the SPP cross section, a total of 18 baryonic resonances are included in addition to a nonresonant process in the mixed isospin channels.Interference between the resonances is incorporated, but not between the resonant and non-resonant components.The initial-state model for SPP interactions in NEUT is a simple relativistic Fermi gas.
Coherent scattering off nuclei also contributes to the SPP cross section, especially at low four-momentum transfer.In this analysis, NEUT models coherent interactions with the Berger-Sehgal model [67], updated from the RS model [68], and includes Rein's model of diffractive pion production [69].

Deep inelastic scattering
Deep inelastic scattering (DIS) describes neutrino interactions with the quark constituents of nucleons.It is a subdominant process in T2K's oscillation analysis due to the neutrino energy and the single-ring event selections at the FD.The cross section in NEUT is calculated using the GRV98 [70] Parton Distribution Functions (PDFs), which describe the probability to find a quark of a given type with a given value of the Bjorken scaling variables, x and y, inside the target nucleon.Bodek-Yang (BY) modifications [71,72] are made to extend the validity of this approach to the rela-tively low four-momentum transfers, Q 2 1.5 GeV 2 , typical for interactions at T2K.
In NEUT, the modelling of DIS processes begins for interactions where the hadronic invariant mass W > 1.3 GeV.To avoid double counting the aforementioned non-resonant single-pion production, only DIS interactions that produce more than one pion in the final state are considered.The generation of the hadronic state is split depending on W : for interactions with W > 2 GeV PYTHIA 5.72 [73] is used, whilst for W < 2 GeV a custom model interpolating between the Δ(1232) and DIS interactions is employed, described in Sec.V C of Ref. [74].

Final-state interactions
The simulated neutrino interaction events produce an outgoing hadronic system at the interaction vertex inside the nucleus, in addition to the outgoing lepton.These hadrons can undergo final-state interactions (FSI) in the nuclear medium.In NEUT, pion FSI are described using the semi-classical intranuclear cascade model by Salcedo and Oset [75,76], tuned to modern π − A scattering data [44].Nucleon FSI are described in an analogous cascade model [56].Within the cascade, the outgoing hadrons are individually stepped through the remnant nucleus where they can elastically scatter, be re-absorbed, undergo charge-exchange processes, and/or emit additional hadrons which are also stepped through the cascade.Amongst other effects, such cascades allow for SPP events to have no observable pions in the final state after FSI, and for 1p1h interactions to appear as pion production interactions.
Table 3 The parameters included in the 1p1h uncertainty model with their values and uncertainties before the ND analysis.The uncertainties for the removal energy parameters are around their central value and contain the carbon-oxygen and ν-ν correlations described in the text.The first five Q 2 parameters are not externally constrained before the analysis, and are free to vary between 0 and ∞.The units of the Q 2 ranges are GeV 2

Parameter
Central value Uncertainty

Coulomb corrections
Following a charged-current neutrino interaction, the electrostatic interaction between the remnant nucleus and the outgoing charged lepton can cause a small shift in the lepton's momentum.The size of this Coulomb correction has been determined from the analysis of electron scattering data [77] and is implemented as a small nucleus and lepton-flavour dependent shift in the momentum of the outgoing lepton.The size of this shift is −3.6 MeV (+2.6 MeV) for outgoing μ − (μ + ) from a carbon target, and −4.3 MeV (+3.3 MeV) for outgoing μ − (μ + ) from an oxygen target.

Uncertainty parametrisation
Mismodelling of neutrino interactions can bias the measurements of oscillation parameters -for instance attributing an increase in single-ring events to an increase in 2p2h interactions instead of CCQE interactions.It is crucial to evaluate the impact that plausible variations of NEUT's interaction model can have on the neutrino oscillation analysis.This section describes the chosen parametrisation of such variations and the corresponding parameters' uncertainties.When possible, theory-driven uncertainties are used, but in many cases this offers insufficient freedom to describe available data, and additional empirically driven parameters are required.To cover the caveats of such an approach, and to consider plausible model variations not included in the model parametrisation, a variety of simulated data studies are per-formed.These are detailed in Sect.5.3, and applied to the oscillation analysis in Sect.9 and Appendix B.

1p1h uncertainties
The 1p1h uncertainty model is split into three categories: removal energy related to the initial state described by the SF, the neutrino-nucleon interaction, and ad hoc freedoms in Q 2 from nuclear effects, amongst others, inspired by external data.The central values and uncertainties are summarised in Table 3.

Removal energy:
A mismodelling of nucleon removal energy would directly bias the reconstructed neutrino energy, which would subsequently bias the extraction of the neutrino oscillation parameters, notably Δm 2 .This was identified as a leading source of uncertainty in a simulated-data study in the last T2K oscillation analysis [1,2].In this analysis, a more reliable modelling of removal energy with accompanying uncertainties was developed.
Unlike the simplistic Fermi-gas models used in the previous iterations of T2K's neutrino oscillation analyses, the SF model does not have a single fixed value for the nuclear binding energy that can be varied as a parameter.Instead, the SF removal energy distribution, extracted largely from exclusive electron scattering data, reflects the shell structure of the nucleus, shown earlier in Fig. 8.The positions of the removal energy peaks, used as an input to the SF model, are measured with a resolution of 2−6 MeV [78] and lower [79].Measurements of the peak positions for carbon differ by up to 2 MeV for the s-shell and 6 MeV for the p-shell [58].The relative strength of each peak also has an uncertainty of up to 10% for carbon [43,80].To extract a SF from (e, e p) data, the impact of nuclear effects such as FSI must be incorporated, and an uncertainty of 5 MeV in this correction is applied [58].In view of these uncertainties, a global removal energy shift uncertainty of 6 MeV is included in the analysis alongside a 3 MeV uncertainty on the difference between the carbon and oxygen removal energies.Further uncertainties are accounted for by the introduction of parameters that allow freedom as a function of Q 2 , described in more detail below.
The construction of the SF from (e, e p) data, and the associated uncertainties, can only be directly applied to modelling 1p1h neutrino interactions with initial-state protons, i.e. antineutrino CCQE interactions.The SF for initial-state neutrons cannot be directly constrained in the same way and the implementation in NEUT assumes that protons and neutrons have the same removal energy distributions.However, as can be seen in Fig. 8, nuclear shell models predict that this is not the case.Calculations suggest that proton and neutron ground states differ in their removal energy by 1−4 MeV, depending on the shell and target [58].For the sharper p-shells, where an energy shift is more important relative to the width of the shell, the offset between the SF and the model calculations for neutrons is around 4 MeV for oxygen and 2 MeV for carbon.To account for this, the central value removal energies of the SF for neutrino interactions are shifted by these amounts, and an uncertainty of 4 MeV is applied on the difference between neutrino and anti-neutrino removal energies.
The removal energy shifts are encoded in four parameters depending on whether they affect initial-state protons (ν CCQE interactions) or neutrons (ν CCQE interactions), and whether the target is carbon or oxygen: The removal energy parameters shift a CCQE event's outgoing lepton momentum and depends on the event's lepton kinematics, neutrino energy, and neutrino flavour.
"Low Q 2 "parameters: NEUT's cross section for chargedcurrent interactions leaving no mesons in the final state (CC0π ) interactions must be suppressed at low Q 2 to match recent measurements from MINERvA [81,82] and T2K [83,84].This is often applied as a suppression of the CCQE cross section via the inclusion of a nuclear screening effect using the Random Phase Approximation (RPA) [60].However, such effects are not included in the SF CCQE model used in this analysis.Since the SF model is built largely on the impulse approximation -which is expected to break down at low momentum transfers 400 MeV/c [54] -extra uncertainties are added in the region where discrepancies with measurements are observed.
The low Q 2 suppression is implemented as five parameters which alter the normalisation of the CCQE cross section in a particular Q 2 range.The parameters span Q 2 = {0, 0.25} GeV 2 and are split into sub-ranges of 0.05 GeV 2 .Since the origin of this low Q 2 suppression in SF predictions is poorly understood, these parameters do not have an external constraint.Whilst this free parametrisation is effective at facilitating a ND-driven modification to the CCQE cross section, the lack of a theoretical basis limits the model's overall predictive power.Several simulated data studies are therefore discussed in Sect.5.3 to evaluate the bias from this technique in the extraction of neutrino oscillation parameters.

M Q E
A and "high Q 2 " parameters: The nucleon axial mass, M Q E A , is tuned to neutrino-deuterium scattering data in NUI-SANCE [85].CCQE cross-section data from ANL [86,87], BNL [88], BEBC [89], and FNAL [90] is used, and deuterium nuclear effects [91] and flux uncertainties for ANL and BNL are included.The central value and its uncertainty are adjusted and inflated to cover the result and previous global fit results [92], giving M Q E A = 1.03 ± 0.06 GeV.Uncertainties on the higher Q 2 > 0.25 GeV 2 predictions of the SF model are driven by the axial component of the neutrino-nucleon interaction, where the dipole model may be inadequate [93].An additional three "high Q 2 " parameters are added to allow an ad hoc freedom, with the goal of Fig. 9 Cross-section predictions for ν μ (solid) and ν μ (dashed) 2p2h interactions on 12 C from Martini et al. [96], Nieves et al. [60], and SuSA v2 [97,98] lessening the extent to which M Q E A is used as an effective parameter to correct for deviations from the dipole model.The Q 2 ranges and uncertainties of the new high Q 2 parameters are based on comparisons of the Q 2 shape of the dipole and z-expansion models [93].

2p2h uncertainties
The uncertainties related to 2p2h interactions are similar to those in T2K's previous oscillation analysis [1,2].Parameters altering the 2p2h normalisation independently for neutrinos and anti-neutrinos, and for carbon and oxygen interactions, are used.The 2p2h normalisations are unconstrained, and the carbon-oxygen scaling parameter has a 20% prior uncertainty.A separate shape uncertainty is also applied, which allows shifts in the Δ and non-Δ contributions in the energy and momentum transfer to the nucleus, (q 0 , |q|), of the Nieves model, also separated for carbon and oxygen interactions.
This analysis also includes additional new uncertainties that reflect the shape of the energy dependence of 2p2h using three different plausible models of the process, also studied by T2K cross-section analyses [94,95].The uncertainties span the maximal difference in 2p2h predictions from Martini et al. [96], Nieves et al. [60], and SuSAv2 [97,98], shown in Fig. 9. Four parameters are added which control the shape of the energy dependence of 2p2h below and above E ν = 600 MeV, and are separately applied to neutrino and anti-neutrino events.

Single-pion production uncertainties
The uncertainty treatment for SPP remains almost identical to previous T2K analyses [1,2,99].There are three central parameters in the modified RS model: the resonant axial mass, M R E S A ; the value of the axial form factor at zero transferred four-momentum, C A 5 (Q 2 = 0); and the normalisation of the I 1/2 non-resonant component.As for M Q E A , the parameters have been tuned to deuterium bubble chamber data using NUISANCE [85], selecting SPP data from ANL [100,101] and BNL [102,103], including corrected data [104].The uncertainties are inflated so that the model adequately describes the SPP cross section in different hadronic mass regions from ANL and BNL, and SPP cross-section measurements on nuclear targets from Mini-BooNE [105][106][107] and MINERvA [108][109][110][111].
A new parameter was introduced for anti-neutrino interactions producing low momentum pions, which constitute a background for the single-ring ν-mode samples.This extra freedom is added through an I 1/2 non-resonant normalisation parameter that affects both ν μ and ν e single pion interactions with p π < 200 MeV/c in the Rein-Sehgal model.The parameter is not constrained by the ND and has an uncertainty of 100%.
Normalisation parameters on the CC and NC coherent cross sections are included separately, and each is assigned an uncorrelated 30% uncertainty based on comparisons to MINERvA data [112].The uncertainty on coherent scattering is fully correlated between carbon and oxygen.

Deep inelastic scattering uncertainties
DIS interactions make a small contribution to the samples in this oscillation analysis due to T2K's neutrino energy.Nevertheless, uncertainties that cover variations in muon kinematics from CC DIS interactions are needed for the ND fit, whose selections contain some multi-π events, and have been significantly updated from previous analyses [99].
As discussed in Sect.5.1, NEUT uses PDFs with BY corrections to calculate the DIS cross section.The uncertainty in the BY corrections is parametrised as a fraction of the difference between using the GRV98 PDFs with and without the BY corrections.At Q 2 > 1.5 GeV 2 the impact is marginal, but in the peak region at lower Q 2 the impact is large, altering the predicted cross section by ∼ 40%.This parameter is split for W < 2 GeV (multi-π ) and W > 2 GeV (DIS) interactions.
Another parameter is introduced to modify the generation of the hadronic state for W < 2 GeV DIS interactions, which uses a custom model [55] to choose the particle multiplicities in an event.This parameter accounts for the differences between the custom model and the AGKY model [113] used in the GENIE event generator [114].
Two normalisation uncertainties are also included, motivated by comparing the NEUT CC-inclusive cross section to the world average of measurements at higher neutrino energies [11].The uncertainties are 3.5% for neutrino interac-tions and 6.5% for anti-neutrino interactions, and the two are uncorrelated.

Final-state interactions uncertainties
The NEUT pion cascade model has been tuned to better match external π − A scattering data [115].The tuning procedure constrains the probability for different interaction processes to occur in the pion cascade (e.g.pion absorption or charge exchange), and is notably more robust than previous parametrisations.The constraints on the pion FSI cascade from the ND analysis are propagated to the FD in this analysis, which was not done before.Furthermore, the simulations at the ND and the FD now use a consistent model for pions from the interaction vertex propagating through the nucleus ("pion final-state interactions"), and for pions propagating through the detector ("pion secondary interactions"), mentioned later in Sect.6.2.The ND constraint on the FSI parameters is only used to constrain the FD modelling of FSI and not the FD modelling of secondary interactions.

Other uncertainties
Additional uncertainties are applied to processes with small contributions to the analysis.As in previous analyses, the NC1γ production cross section has a 100% normalisation uncertainty.The NC elastic, NC resonant kaon and eta production, and NC DIS interactions are grouped together and referred to as "NC other" interactions, which have a 30% normalisation uncertainty that is uncorrelated at ND and FD.There is one uncertainty controlling the normalisation of the electron neutrino cross section, and another controlling the electron anti-neutrino cross section.The uncertainties are composed of two parts: one 2% uncorrelated part and one 2% anti-correlated part, which connects the two parameters [116].The parameters only affect electron (anti-)neutrino interactions, and have no effect on the other neutrino flavours.The total cross sections of CC resonant single-photon production, CC resonant kaon production, CC resonant eta production, and CC diffractive pion production are controlled by a single new parameter referred to as "CC misc", which is a 100% normalisation uncertainty, and such interactions are not affected by other model parameters.Two new parameters are included to account for Coulomb corrections [117,118].They control the normalisation of the (anti-)neutrino cross section for E ν = 0.4−0.6GeV with a 2%(1%) uncertainty, and are 100% anti-correlated.

Simulated data studies
The systematic uncertainties in the analysis are constructed to account for known uncertainties in neutrino interaction physics, but can not possibly cover every model scenario.
For instance, cross-section measurements from T2K and other experiments have shown that no single 1p1h model describes the kinematic phase space in T2K and MIN-ERvA [81,94,95,[119][120][121][122].In addition, the ND analysis, presented later in Sect.6, may compensate for cross-section mis-modelling by varying the flux parameters instead of the cross-section parameters, leading to good agreement with the observed event spectrum in lepton kinematics.However, the fitted model may scale the effect incorrectly in other important physics variables, e.g.E ν .It is therefore crucial to test whether the uncertainty model is flexible enough to capture variations under alternative cross-section models which are not directly implemented in the default uncertainty model, and whether the subsequent extrapolation of model constraints to the FD has an effect on constraining the oscillation parameters.
Some of the simulated data sets are similar to those presented in T2K's previous analyses [1,2].The studies are updated due to the significant changes in the uncertainty model and ND analysis.The alternative models and tunes are selected to cover a number of interaction types and effects, listed next.
CC0π simulated data sets: The dominant CC0π samples at the ND and the single-ring samples at the FD are designed to select CCQE-like events.The larger statistics in these samples requires testing for a range of alternative models, and the robustness of the neutrino interaction model.
• Non-CC-Quasi-Elastic (non-CCQE) contributions -Before the fit to data, the prediction of the CC0π selection at the ND is underestimated by 0−20%, depending on the outgoing lepton kinematics.Projecting the data and prediction onto the reconstructed four-momentum transfer, Q 2 rec , defined as the Q 2 calculated for a CCQE interaction on a stationary nucleon, and with a binding energy E b , the discrepancy is less than 5% at Q 2 rec < 0.1 GeV 2 and approximately 20% for higher Q 2 rec .The CCQE cross section is modified after the fit to ND data to account for the difference.This simulated data tests the hypothesis that the underestimation of data is actually due to non-CCQE contributions, and does so by scaling up their predictions instead of the CCQE components.The study is given in detail in Appendix B.
• Alternative CCQE form factors -The baseline model used in this analysis assumes a dipole parametrisation of the nucleon form factor.There are other form factor models, of which the 3-component (an extension of Ref. [123]) and z-expansion [93] formalisms were tested.
The effect is largely expected to be covered by the Q 2related freedoms of the cross-section model.• Multi-nucleon (2p2h) model -The Nieves et al. model [60] was used to describe 2p2h interactions in this analysis.
An alternative 2p2h model from Martini et al. [96] was tested in the simulated data studies, because its 2p2h cross section is larger and evolves differently in E ν for neutrinos and anti-neutrinos, shown earlier in Fig. 9. Modelling the 2p2h spectrum is important in the ND to FD extrapolation, as it is one of the main sources of bias in the reconstructed neutrino energy spectrum of CCQE-like samples.The SuSAv2 model [97,98], also shown earlier in Fig. 9, is a less extreme variation compared to the Martini model, so was not included.• Removal energy -The nuclear removal energy in the relativistic Fermi gas (RFG) model [124] was the largest contributor to uncertainty in the previous T2K analysis [1,2].This analysis' spectral function (SF) model [43], mentioned earlier in Sect.5.1, introduced an improved parametrisation for the removal energy uncertainty, and simulated data sets were developed to study its impact.
CC1π simulated data sets: Single-pion events are a background for the single-ring selections at the FD and contribute to the bias in reconstructed neutrino energy.Additionally, the 1Re1de sample at the FD specifically targets single-pion events, which motivates the need to have a robust uncertainty model of these interactions.Three simulated data sets were produced: • ND data-driven pion momentum modification -The 1Re1de selection at the FD tags low momentum pions below Cherenkov threshold by the presence of a delayed Michel electron.The ND analysis in Sect.6 uses selections based on muon kinematics and pion tagging to constrain the uncertainties, and does not study the pion kinematics directly.As such, single-pion events may be modelled well in muon kinematics and poorly in pion kinematics.A data-driven simulated data set was created by studying the CC1π selections at the ND, using the model that was fit to ND data in lepton kinematics.The model was used to predict the reconstructed pion momentum spectrum, p reco π , in the single-pion ND selections, which was compared to the data in the p reco π < 200 MeV/c region.The number of events was underestimated by ∼ 20%, which was applied as an overall normalisation to the simulation of all single-pion events at the FD that had a pion with generated (true) momentum below 200 MeV/c.This is the only simulated data set that was not applied at the ND, and tested only at the FD.
• MINERvA pion suppression -A low-Q 2 suppression of the single-pion production cross section in GENIE [114] was needed to consistently describe neutrino interactions on plastic scintillator (CH) from MINERvA and bubble chamber data on nucleons [125].The function parametrising this discrepancy was used to create simulated data at both the ND and the FD and the study is presented in detail in Appendix B. • Pion secondary interactions -This analysis introduced a new model for pions rescattering in the ND.The GEANT4 model [49] was replaced with NEUT's Salcedo-Oset model [75,76] which was tuned to π − A scattering data [44].A hybrid simulated data set which blended features of the two models was used in the ND analysis to study the impact of choosing one model over the other.
A summary of the simulated data studies is presented in Sect.9 after the analysis sections, and the simulated data studies are detailed in Appendix B.

Near-detector analysis
The high statistics data at the ND are used to constrain many of the neutrino flux and neutrino-nucleus interaction models present in the neutrino oscillation analysis.Sampling the unoscillated neutrinos at a high rate and tuning the prediction to the ND data allows for significant reduction of the uncertainties of the FD prediction.The ND analysis targets CC0π events as these are the signal at the FD, and additionally constrains the background contributions such as CC1π and CC multi-π events.Separation of ν μ and ν μ events is possible due to the magnetised sign-selecting ND, and there are ν μ selections in the ν-mode which constrain the wrong-sign background.
As in previous T2K oscillation analyses [1,2], two complementary likelihood sampling methods are used and are cross-validated.One is based on Markov Chain Monte Carlo (MCMC) methods [126,127], and the other is based on minimising a test-statistic through gradient-descent methods in Minuit [128].The MCMC analysis is inherently Bayesian, and has the ability to run a simultaneous ND+FD analysis whose results are presented in Sect.8.1.The gradient-descent analysis instead fits the systematic uncertainties in the simulation to find the global minimum of the test statistic that best describes the data at the ND, discussed in Sect.6.3.The central value and covariance matrix of the systematic uncertainties around that best-fit point is then propagated to the FD.The MCMC framework directly implements the removal energy shift parameters described in Sect.5.2 which allows for discrete event migrations between bins, whereas the gradient-descent framework smooths the effect by an effective binned treatment to avoid discontinuous likelihoods.The MCMC analysis also implements a non-uniform rectangular binning, meaning the binning in the x variable is not uniform in the y variable, allowing the events to be binned finer and more effectively, generally leading to improved sensitivity to the systematic uncertainties.The gradient-descent framework instead uses a uniform rectangular binning.The effect of these differences is tested at the FD by propagating the results from the gradient-descent framework, which assumes correlated Gaussian parameter constraints, and comparing to propagating the constraints from the steps in the MCMC, which is detailed in Sect.8.3.This section shows the results from the gradient-descent based analysis.
This analysis of ND data uses 19.867 × 10 20 POT, with 11.531×10 20 collected in ν-mode and 8.336×10 20 collected in ν-mode, as listed earlier in Table 1.This is an overall POT increase of 106% compared to the previous analysis.

ND selections
The doubling of ν-mode data in the ND allowed for a refining of the anti-neutrino selections.Additionally, the ν-mode beam samples now match the ν-mode beam samples in the separation of events by their reconstructed pion multiplicity.Previous analyses only split the ν-mode selections into events with a single muon-like track (CC-1Track) and events with a single muon-like track with at least one charged or neutral pion candidate (CC-NTrack).
The events are categorised into 18 samples, split into nine equivalent FGD1 and FGD2 samples to separate neutrino interactions on plastic scintillator (FGD1), and plastic scintillator and water (FGD2).The samples first require a reconstructed muon to be present.They are then split by the sign of the muon candidate -which implies the identity of the incoming neutrino -classifying events as ν μ events in νmode, ν μ events in ν-mode, and ν μ events in ν-mode.Each of these charged-current inclusive selections are separated into three reconstructed topologies based on the number of reconstructed charged pions.An event with no reconstructed pions is classified as CC0π ; an event with a single charged pion with opposite charge to the muon is CC1π ; and an event with any other number of charged pions (e.g.1μ − 2π + or 1μ − 1π − in ν-mode), or at least one neutral pion, is classified as CC other.There is no requirement on the number of proton tracks and there is no dedicated ν e or ν e selection.
The pion tagging in the ν μ selections is the same as in previous T2K analyses [1,2].A pion is tagged by either a pion-like track in the TPC, a pion-like track contained in the FGD, or an isolated delayed Michel electron in an FGD.In the FGD and TPC tagging, the pion candidate is required to share its vertex with the muon candidate, and for the Michel tag it is required to be in the same FGD as the candidate Table 4 Efficiencies and purities for each of the selections at the ND in this analysis, including wrong-sign background components.The efficiency is defined as the number of events that have a reconstructed selection that matches the true selection, divided by the total number of events with that same true selection.The purity is defined as the number of events with the desired selection divided by the total number of events in the selection The efficiencies and purities are determined from reconstructed simulated events, and are provided in Table 4, which shows similar performance for the two FGDs.FGD2 has worse Michel tagging and FGD-contained track reconstruction than FGD1 due to the passive water layers, resulting in a lower efficiency for CC1π selections.The purity for CC0π selections for ν-mode and ν-mode is above 70%, and ∼ 55% for the ν μ in ν-mode due to the wrong-sign neutrino flux having a longer tail, which makes multi-particle final states more likely.The ν μ CC0π efficiency is higher than ν μ CC0π due to ν μ CCQE interactions usually producing a neutron in lieu of the proton from ν μ CCQE interactions.In ν μ CCQE interactions, the outgoing proton may produce a clear track in the detector, which has a probability of being mis-tagged for a π + (or μ + for ν μ selections), and so enters another selection; this is very unlikely when the outgoing particle is a neutron.Furthermore, ν μ interactions generally produce a larger proportion of forward-going events, where the ND has better acceptance.
The ν μ CC other selections' low purities compared to the ν μ in ν-mode and ν μ in ν-mode equivalents stem from the larger wrong-sign background that, for the reasons stated earlier, produces multiple pions which may be wrongly selected as the μ + candidate.In addition, the muon candidate in νmode can be incorrectly assigned as a high momentum proton around p ∼ 1 GeV/c, where the energy loss in the TPC for a proton is similar to that of a muon.This track confusion seldom happens in the ν μ selections, since it selects a negatively charged track.A π − is rarely selected as the μ − in ν μ selections since it requires a higher energy multi-π event or final-state interactions of a hadron from the primary interaction.
Generally, the mis-identification of the muon candidate is largest at low momentum, when it does not leave a long enough track to reliably assess the degree of bending in the ND's magnetic field.Almost all wrong-sign muons, pions and electrons selected as the muon candidate occupy this region.In the case of mis-identification, the muon candidate is otherwise a pion with same charge due to their similar energy loss in the FGDs and TPCs.Using the combined FGD+TPC detector system, there is a 94%, 86%, and 77% probability that the muon candidate is a muon in the CC0π , CC1π and CC other selections, respectively.

ND related uncertainties
Dedicated control samples have been developed to evaluate the response of the ND and to quantify systematic uncertainties [129].These uncertainties include the modelling of pion and proton secondary interactions in the detector, particle mis-identification probabilities in the TPCs and FGDs, magnetic field distortions, momentum resolutions and scales, efficiencies related to clustering, tracking and track matching, Michel-tagging efficiencies, pile-up, FGD mass, out of fiducial volume (OOFV) background events, and sand muon backgrounds.Sand muon backgrounds enter the selections when neutrinos from a beam spill interact in the sand surrounding the ND pit, creating a muon that enters the ND.These uncertainties can migrate events into or out of selections and change the reconstructed particles' kinematics.The uncertainties can either be efficiency-like (dependent on a particle's kinematics) or normalisation-like (independent of a particle's kinematics).This analysis is the first to use NEUT's semi-classical Salcedo-Oset cascade model [75,76], mentioned in Sect.5.1.5,for pion secondary interactions in the detector, where previous analyses used GEANT4 [49].The model was tuned to external π − A scattering data [44], and was found to agree better with data and be more consistent across the interaction channels and pion energy ranges compared to GEANT4.Additionally, T2K now uses the same model for pion final-state and secondary interactions in both the ND and the FD.The ND constraint on pion final-state interactions is propagated to the FD, whereas the constraint on the secondary interactions is not.
The uncertainties from the detector uncertainties are presented in Table 5, and are 1.2−2.1% for the CC0π selections, and 2.5−4.0%for the CC1π and CC other selections.The secondary interaction uncertainty for pions contribute 70−95% of the total detector-related uncertainties, depending on the selection.For reference, the statistical uncertainty on the number of events in the ND selections, presented later in Table 6, is 0.5−1.3% for the ν-mode selections, and 1.1−3.9%for the ν-mode selections.

Defining the likelihood
Each selection is binned in the reconstructed muon momentum, p μ , and the cosine of the muon angle with respect to the detector z-axis, cos θ μ , which nearly lines up with the average neutrino direction. 1The ND likelihood is constructed by calculating the −2 ln L total of the data and simulation (MC) across all bins in all samples at each set of the parameter values.The systematic uncertainties in the models for the ND response, neutrino interactions, and neutrino flux, detailed in previous sections, are encoded via a Gaussian penalty term, which includes the covariances between the systematic uncertainties, shown in Eq. 7. The treatment of statistical uncertainties in the simulation has been updated [130,131] and was validated against a complementary approach [132] and the previously used method.The total likelihood is defined as where L stat is the statistical likelihood, L MC stat is the MC statistical uncertainty likelihood, and L syst is the likelihood of the systematic uncertainties.The frequentist analysis maximises L total by finding the minimum of −2 ln L total , and the Bayesian analysis samples the −2 ln L total around the minimum in proportion to the posterior probability.The first two terms in Eq. 5 are linked, as the statistical uncertainty on the MC affects the number of MC events.The two statistical contributions read, where in each bin j of sample i, N Data (N MC ) is the number of events in data (MC), β j scales the unweighted MC events, and σ β j is a measure of the MC statistical uncertainty.The systematic uncertainties are parametrised as correlated Gaussian penalties, where x ( μ) are the values of the systematic uncertainty parameters during (before) the fit, and V is their covariance matrix.The ND constrains the flux uncertainty at the FD through such a covariance matrix.The low-momentum ν μ SPP, neutrino energy-dependent 2p2h, NC other, NC1γ, and ( ν ) e / ( ν ) μ parameters are barely constrained by the ND analysis, so their constraints are not propagated to the FD in the frequentist analysis.In the simultaneous ND+FD Bayesian analysis, both detectors are used to constrain these parameters.

Results of the ND analysis
The ND analysis sees large shape changes in the ν-mode ν μ flux parameters with roughly 10% enhancement at low E ν and 10% suppression at high E ν , as shown in Fig. 10.The neutrino flux parameters have strong correlations with each other and with some cross-section parameters, such as M Q E A and the Q 2 parameters, shown in Fig. 13.Moving the flux parameters by this amount incurs a penalty of −2 ln L flux /N dof ∼ 1 for this variation in flux parameters Fig. 11 Constraints on the CC0π parameters, excluding the CCQE Q 2 parameters, from the fit to ND data (black points, black lines), overlaid on the input uncertainty (red band).The parameters on the left-hand side of the figure are presented as a ratio to the generated value in NEUT, and the right side shows the removal energy parameters, E rmv , with shifts in units of MeV.CCQE interactions are generated in NEUT with 21 GeV, but a pre-fit value of 1.03 GeV was used after analysis of CCQE bubble chamber data.The absence of an uncertainty band reflects that the parameter was not constrained by external inputs before the analysis of ND data Fig. 12 Constraints on the CCQE Q 2 parameters as a function of Q 2 from the fit to ND data (black points, black lines), overlaid on the input uncertainty (red band).The absence of an uncertainty band reflects that the parameter was not constrained by external inputs before the analysis of ND data due to the large correlations, confirmed by p-value studies in Sect.6.6.
Figures 11 and 12 show the CC0π cross-section parameters after the fit.Despite the external constraint on M Q E A , the data prefers a larger value of M Q E A = 1.16 GeV.A complementary fit, changing the uncertainty on M Q E A to be unconstrained instead of informed by bubble chamber data, had little impact on the ND analysis and the predictions at the FD; hence the constraint on M Q E A is primarily driven by the ND data.The 2p2h normalisation is different for neutrinos and anti-neutrinos, which are both constrained to ∼ 15% uncertainty, with the 2p2h normalisation for neutrinos consistent with the prediction from Nieves et al.The 2p2h normalisation for carbon and oxygen is consistent with 1, although the shape parameter for oxygen agrees with the Nieves model, whereas the carbon parameter is pulled to be more Δ-like, differing by ∼ 1σ.The removal energy parameters are within their uncertainties before the fit, and are compatible for the carbon, oxygen, neutrino and anti-neutrino parameters.
The CCQE Q 2 parameters are shown in Fig. 12, where there is a suppression at low Q 2 until 0.2 GeV 2 , consistent with other cross-section data mentioned in Sect.5.2.At higher Q 2 the data prefers an enhancement of 20−30%.The Q 2 parameters have strong anti-correlations with the flux parameters, as shown in Fig. 13, and studies with fixed values of the Q 2 parameters showed that the flux parameters compensate for differences in Q 2 for CCQE events, a testament to the parameters' correlations.
The 2p2h normalisation has been given the freedom to independently vary for neutrino and anti-neutrinos, and differences in 2p2h neutrino and anti-neutrino parameters may reflect a more general mismodelling of CC0π interactions.This may allow deficiencies in the anti-neutrino CCQE model to be absorbed in the 2p2h normalisation parameters.Similarly, M Q E A and the CCQE Q 2 normalisation parameters may be absorbing effects from a different axial form factor parametrisation, which may evolve differently as a function of other variables, e.g.E ν , as mentioned in Sect.5.3.Both of these effects, amongst others, are studied through simulated data studies in Sect.9 and Appendix B. The full parameter set with their values before and after the analysis of ND data is provided in Appendix E.
The MCMC and gradient-descent analyses differed in the treatment of the removal energy uncertainty.The MCMC allows for discrete movement of events between bins, which may produce multi-modal posterior probability distributions (output constraint) of the removal energy parameters.The smoothed binned implementation in the gradient-descent framework prevents this from disrupting the ability to find the maximum likelihood, whilst still capturing the overall physics behaviour of the removal energy uncertainty.The impact of this and other differences between the analyses, such as the non-uniform rectangular binning scheme, were addressed by separately propagating the covariance matrix from the gradient-descent framework and the parameter variations sampled by the MCMC to the oscillation analysis in Sect.8.
In general, the constraints on the parameters and the impact of the ND analysis agrees with the expected sensitivity.Furthermore, compatible results are found between the MCMC and the gradient-descent analyses in the central value estimates, uncertainties, and correlations of the parameters, leading to consistent sample predictions at the ND and the FD.

ND predictions
The aforementioned selections in the data and simulation are compared before and after fitting to data, using the constraints on the systematic uncertainties from Sect.6.4.Table 6 shows the number of events in each selection, where the agreement between the post-fit prediction and the data is notably improved compared to that of the pre-fit prediction, especially for the CC0π events, which comprise the main signal at the FD.There is a consistent rise across all CC0π selections and a small suppression of ν-mode 1π + events, improving agreement with the data.This causes the smaller ν-mode 1π − prediction to also be suppressed, since they share parameters in the interaction model, with the neutrino flux and detector uncertainties being more loosely correlated, connected only through their input covariance matrices.
The observed and predicted ν-mode ν μ FGD1 CC0π events projected onto p μ are shown in Fig. 14 before and after the fit to data.Before the fit, there is a notable underprediction which is largest at low p μ .The fit increases the CCQE and 2p2h components and decreases the 1π components in the prediction to agree with the data.For comparison, the ν-mode ν μ FGD2 CC0π selection is also shown in Fig. 14, where there is agreement between the prediction and the data before the fit, which marginally improves after the fit.This showcases the ability of the systematic uncertainty treatment in the analysis to modify and constrain the modelling of neutrino and anti-neutrino interactions on carbon and oxygen separately, and the strength of having a signselecting ND.

Assessing model compatibility with data
A p-value is calculated to assess the probability of the model given the data, and represents the probability that a model with a test statistic equal to or larger than the observed data is found.Simulated data sets, referred to as "toys", are created by varying the systematic uncertainties in the model according to their input covariances before the ND analysis, and statistical fluctuations are applied.The model is fit to each toy and the (−2 ln L ) min is calculated.The pvalue is defined as the fraction of the simulated data sets with (−2 ln L ) Toy min ≥ (−2 ln L ) Data  min .An a priori criteria of p > 0.05 is required of the ND analysis for the results to be used in the oscillation analysis.Using a total of 895 simulated data sets, p = 0.74, demonstrating good agreement between the model and the data.
Breaking down the (−2 ln L ) min contributions by the likelihoods from the selected samples and systematic uncer-Fig.13 Correlations between selected ν-mode ν μ FD flux and CCQE cross-section parameters.The flux and Q 2 normalisation parameters' ranges are in units of GeV.The strong anti-correlations between the flux and cross-section parameters significantly reduce the uncertainties on the predictions at the FD tainties in Table 7, the selected samples are generally described well with p = 0.82, with individual p-values for the CC0π selections between p = 0.15−0.93.Splitting the neutrino flux contributions into ν-mode ν μ , ν-mode ν μ , ν-mode ν μ and ν-mode ν μ , p = 0.74, 0.74, 0.31, 0.37 respectively, showing good compatibility.The cross-section systematics are the worst contributor with p = 0.01, coming predominantly from parameters that are pulled away from their external constraints, e.g.
When instead varying the systematic uncertainty parameters with respect to their constraints after fitting to data, the crosssection model p-value improves to approximately p = 0.3.This indicates that the cross-section model before the fit to data is unfavourable, but after the fit to data is satisfactory.The near-detector analysis constrains the product of the neutrino flux, ND detector, and neutrino interaction uncertainties, leading to large correlations between the systematic uncertainties, as demonstrated in Fig. 13.Therefore, studying one group's p-value in isolation from the other is not exact.For this reason, the p-values from the uncertainty parameters do not have to follow the same strict criteria of p > 0.05.However, the low p-value does highlight the need for contin-ued effort in developing realistic neutrino interaction models and associated uncertainties.

Far-detector selection
The FD event selection in this analysis is the same as used in previous T2K results [1]; only the data have been updated, and the selection is briefly reviewed here.Similarly, the method of evaluating systematic uncertainties related to the FD is unchanged from previous analysis, where atmospheric events in SK are used to calculate the uncertainties using a MCMC-based approach.
The event reconstruction in SK uses both charge and timing information from hits in the PMTs, and particles are detected using their Cherenkov rings.The vertex position, momentum, and particle type of each ring is reconstructed [133].Muons and electrons are differentiated by their ring profiles, where muons generally produce "sharper" rings due to less scattering, and electrons produces "fuzzier" rings due to their electromagnetic showers.All samples in this analysis are based on observing one electron-like (1Re)  The right is a similar comparison showing the parent muon's particle ID parameter, which separates events into electron-like (positive values) and muon-like (negative values) categories.The uncertainty on the data points is statistical or muon-like (1Rμ) primary Cherenkov ring, and a specific number of delayed triggers relative to the primary interaction, consistent with a Michel electron from an unseen charged pion's decay chain (referred to as decay electron, or "de").Three samples are selected in the ν-mode data: a CCQE-like ν e sample (ν-mode 1Re with 0 de), a CCQE-like ν μ sample (ν-mode 1Rμ with 0 or 1 de), and a CC single pion-like ν e sample (ν-mode 1Re with 1 de).Similarly, there are two single-ring ν-mode data samples: a CCQE-like ν e sample (νmode 1Re with 0 de) and a CCQE-like ν μ sample (ν-mode 1Rμ with 0 or 1 de).Unlike the ND, the FD is not magnetised and can therefore not determine the charge of the outgoing particles.
Since the start of T2K operations in 2009, the gain of the SK inner detector's PMTs has increased at a rate of at most a few percent per year.In previous T2K analyses, this effect was corrected during the reconstruction stage using a run-byrun global correction factor for all PMTs.However, the gain drift differs based on the PMT production year, and the current analysis adopts a more detailed correction that accounts for these differences.All T2K FD data in this analysis have been reprocessed and reconstructed using the updated correction.The change to the gain correction results in a change in the observed charge available to the reconstruction algorithm relative to previous analyses, even when processing the same event.This may cause small shifts in an event's reconstructed parameters, including the number of rings, and each ring's particle type and momenta, which has caused some events to migrate into or out of the oscillation analysis samples with respect to the previous analyses.For the reprocessed run 1−9 data there are in total 1 more ν-mode 1Re, 1 fewer ν-mode 1Re1de, 1 more ν-mode 1Re, and 3 fewer ν-mode 1Rμ events compared to previous oscillation analysis.The Table 8 Summary of event migrations at the FD after reprocessing data from the previous T2K analysis [1,2]."Inward" refers to newly added events that were not present in the previous analysis, "outward" refers to events that were lost to the update, and "overlap" refers to the number of events that are common to the two analyses migration of the events is summarised in Table 8.As the gain correction is applied to data and not to the simulation, the event migration has been cross-checked in both atmospheric neutrino and cosmic-ray muon data samples, which are used to evaluate FD detector uncertainties in the T2K analysis.In both studies, the level of migration was found to be consistent with that observed in the T2K beam data.
Fig. 17 The events in the full data set for the five FD samples, shown in reconstructed lepton momentum and the angle between the neutrino beam and the lepton in the lab frame.The coloured background in the two-dimensional plot shows the expected number of events from the frequentist analysis, using the best-fit values for the oscillation and systematic uncertainty parameters, applying the reactor constraint on sin 2 θ 13 .The insets show the events projected onto each single dimension, and the red line is the expected number of events from the best-fit.The uncertainty represents the 1σ statistical uncertainty on the data Fig. 18 The number of ν-mode 1Re + 1Re1de versus ν-mode 1Re events (top, leading sin δ CP dependence) and ν-mode 1Re + 1Re1de + ν-mode 1Re events above and below E rec = 550 MeV (bottom, leading cos δ CP dependence), with the predicted number of events for various sets of oscillation parameters, as shown by the different coloured ellipses.The values for the neutrino mass splitting are from the frequentist analysis of data, where Δm 2 32 = 2.40 × 10 −3 eV 2 (Δm 2 31 = −2.46×10−3 eV 2 ) is the best-fit point in the normal (inverted) ordering.The uncertainties represent the 68% confidence interval for the mean of a Poisson distribution given the observed data point.The underlaid contours contain the predicted number of events 68% of simulated experiments, varying the systematic uncertainty parameters around the best-fit values from the fit to ND and oscillation parameters set to the best-fit values from a fit to data.The overlaid triangle point shows the predicted of events with both oscillation and systematic uncertainty parameters at their data best-fit values This analysis the first to include data following the refurbishment of the FD in after the detector had been prepared for the gadolinium phase [39] but still using the ultrapure water without gadolinium, referred to as the SK-V period.Following this work, T2K's run 10 was under slightly different detector conditions than that of the previous data sets.This period a background rate primarily at O(MeV) energies, irrelevant to T2K's analysis.During the run, the water's attenuation length, as measured by throughgoing cosmic-ray muons, was found to be stable above 90 m, consistent with data taken before the refurbishment, albeit slightly longer.This suggests event reconstruction and detector uncertainties should similarly be consistent between the data periods, and several cross-checks were performed to confirm this.
Figure 15 shows such a comparison between stopping cosmic-ray muon data and their Michel electrons taken during the run 9 and run 10 data periods at SK.The similarity of the distributions over both data sets highlights the stability of the detector and reconstruction algorithm following the refurbishment in 2018.Though only the reconstructed Michel momentum distribution and the parent muon's particle ID parameter are shown in the figure, distributions for other reconstructed parameters used in the T2K event selection showed similar high consistency.Kolmogorov-Smirnov tests of the expected events in run 10 confirmed this.This was true for other calibration data as well as for atmospheric neutrino data, and small differences in these distributions were within current uncertainties.
Good detector stability was also found for the timing and selection of events observed in the T2K beam.The distribution of event times relative to the start of the spill at J-PARC is shown in Fig. 16 for events with minimal outer detector activity, labelled fully-contained events.Events from run 10 showed a 34.2 ns RMS relative to their nearest expected bunch timing (dotted lines in the figure), consistent with that from previous runs.
Amongst the 354 selected fully-contained events in run 10, 75 were selected as 1Rμ, 18 as 1Re, and there were no new 1Re1de events for the analysis described in the next section.The number of events in each selections is presented in Sect.8, Table 9.

Oscillation analysis
This section presents the three-flavour oscillation analysis from the full data set presented in Fig. 17, including the constraints from the ND analysis in Sect.6.The analyses at the FD are first introduced, followed by the constraints on the oscillation parameters from the Bayesian and frequentist data analyses in Sects.8.1 and 8.2, respectively.The comparison of the Bayesian and frequentist analyses are presented in Sect.8.3, and the new result is put in the context of current world data in Sect.8.4.The results presented in this section include the uncertainty inflation procedure from simulated data studies mentioned in Sect.5.3, whose results are discussed in detail later in Sect.9 and Appendix B.
The impact of δ CP on the number of events in the selections is shown in Table 9, where there is a relatively small Table 10 Uncertainties on the number of events in each FD sample broken down by source after (before)the fit to ND data."FD + SI + PN" combines the uncertainties from the FD detector, secondary particle interactions (SI), and photo-nuclear (PN) effects."Flux⊗Interaction" denotes the combined effect from the ND constrained flux and inter-action parameters, and the unconstrained interaction parameters.The change in the "FD + SI + PN" uncertainties before and after the ND fit is an indirect effect due to the change of interaction mode fractions in the samples after the ND fit sensitivity in the ν-mode 1Re selection, and most sensitivity comes from the ν-mode 1Re selection, owing to the number of events in each sample.To summarise the results, the number of observed electron neutrino events are plotted against the observed anti-neutrino events in Fig. 18, where the data favours δ CP ∼ −π/2, Δm2 32 > 0, and sin 2 θ 23 > 0.50; i.e. near maximal CP violation, the normal mass ordering, and the upper octant in the PMNS paradigm.The 1Re +1Re1de events in ν-mode and the 1Re events in ν-mode are sensitive to sin δ CP , the neutrino mass ordering, and the octant of θ 23 , and their energy spectra has some sensitivity to cos δ CP , as illustrated in Fig. 18.Compared to T2K's previous analysis [1,2], the data are now closer to the best three-flavour fit prediction, resulting in a slightly weaker constraint on δ CP .The weaker constraint is, however, more compatible with the expected sensitivity of the experiment, discussed later Sect.8.2.
The systematic uncertainties on the predicted number of events before after the fit to ND data is given in Table 10.After the fit, the total uncertainty is reduced by a factor 2-5 depending on sample, with the impact from flux and interaction uncertainties reduced by more than 60%. the ND fit, the interaction uncertainties are of similar size to the FD detector, pion secondary interaction, photo-nuclear systematic uncertainties for all samples except the 1Re1de, which is dominated by FD detector uncertainties.The FD detector uncertainties characterise the performance of SK and its reconstruction, the pion secondary interaction uncertainties were discussed in Sect.5.2.5 and are informed by external π − A scattering data, and the photo-nuclear uncertainty comes from when Cherenkov photons are absorbed by the nuclei in the FD, causing particles to be mis-reconstructed or entirely missed due to the lack of any Cherenkov rings.Although the impact from uncertainties in the flux and interaction model are similar for the selections at about 3% when considered separately, they significantly correlate with each other after the fit to ND data, which causes the combined uncertainty from the ND-constrained interaction parameters and the neutrino flux to be smaller than the sum of their squares.
These constraints are used to build the predictions for the FD energy spectra including all uncertainties, as shown in Fig. 19.The five lower-Q 2 parameters have no external constraints, and the expected sensitivity from a FD-only fit (excluding the ND) is used as the uncertainty.This is solely for the purpose of providing a representative uncertainty on the events when an ND fit is not used, and this uncertainty is not used elsewhere in the analysis.
The degrees of freedom from the oscillation parameters are of the form sin 2 θ i j , Δm 2 i j , and δ CP .T2K is not sensitive to the "solar" oscillation parameters sin 2 θ 12 and Δm 2  21 , therefore constraints from the world averages reported in PDG 2019 [11] are imposed, 2 where the frequentist analysis fixes the parameters and the Bayesian analysis accounts for their uncertainties.An additional constraint may be imposed on sin 2 θ 13 from the world average reported in PDG 2019 [11], referred to as the "reactor constraint". 3The reactor constraint has a significant effect on the sensitivity to other oscillation parameters of interest, notably δ CP .Accordingly, results are presented with and without this constraint applied.The reactor constraint is applied as a Gaussian penalty to the test statistic for both the frequentist and Gaussian analyses.

Bayesian results
The Bayesian results presented in this section are obtained by sampling the posterior distributions through MCMC [126,127] analysis, using the ND and FD selections simultaneously.The MCMC analysis presented in Sect.6 is utilised for the ND.The e-like samples use both the reconstructed angle between the outgoing lepton and the mean neutrino direction, and the reconstructed neutrino energy assuming a CCQE interaction and a struck nucleon at rest Eq. ( 4).
19 Total uncertainty on the reconstructed neutrino energy spectrum in the FD selections before and after the ND analysis of data.The oscillation parameters are set to values near the T2K best-fit point, specified in Appendix B, Table 17 For the 1Re1de selection -which is dominated by 1e − 1π + final states -the nucleon mass is replaced by the Δ(1232) mass.The μ-like samples only use the reconstructed neutrino energy assuming a CCQE interaction.The posterior probability at the FD first includes the product of Poisson probabilities for observing the number of events in the data given the model prediction per bin across all samples.A Gaussian multivariate distribution is used to include the effect of external constraints on the systematic uncertainty parameters.The general form of the likelihood is the same as the ND analysis, presented in Eq. 6, but excludes the statistical uncertainty on the simulation for the FD.
Credible regions are extracted from lower dimensional marginalised posterior distributions for parameters of interest by adding up the highest probability density region until a certain fraction of the distribution is captured.Flat priors are used over the entire ranges of sin 2 θ 23 , Δm 2  32 , δ CP (or sin δ CP ), and Gaussian priors are applied on Δm 2  21 and sin 2 θ 12 .For sin 2 θ 13 either a flat or a Gaussian prior is applied via the aforementioned reactor constraint.The priors for normal and inverted orderings are the same, namely 50%. Figure 20 shows several marginalised posterior distributions for oscillation parameters of interest.Two-dimensional distributions for every combination of the four oscillation parameters of interest are shown with the 68% and 90% credible intervals in dashed and solid lines, respectively.Each two-dimensional posterior distribution also shows the point of highest probability density.Marginalised one-dimensional posterior probability distributions are also given for each of the four oscillation parameters with 68%, 90%, and 95% credible intervals in different shades of grey.

Atmospheric oscillation parameters
The effects of applying the reactor constraint on the sin 2 θ 23 − Δm 2  32 contours is shown in Fig. 21.Applying the constraint increases the probability density in the upper octant and the normal neutrino mass ordering.The marginalised posterior probability distribution of sin 2 θ 23 with and without the reactor constraint is shown in Fig. 22.The posterior probabilities are largely overlapping, with a preference for the upper octant when using the reactor constraint, and there is barely any octant preference without the reactor constraint.
The results for the atmospheric parameters are summarised in Table 11, showing the proportion of the poste-Fig.21 68% and 90% credible intervals from the marginalised sin 2 θ 23 − Δm 2  32 posterior distribution with (red) and without (blue) the reactor constraint applied.The top (bottom) shows the proportion of probability density in the normal (inverted) mass ordering rior probability that lies in the different mass orderings and θ 23 octant, with and without the reactor constraint.A flat prior distribution on both Δm 2  32 and sin 2 θ 23 is equivalent to comparing the likelihood that T2K's data is described by the different choices of hypotheses.The analysis with (without) the reactor constraint sees a Bayes factor (BF) of 3.35 (1.43) for the upper over the lower θ 23 octant; 4.21 (1.83) for the normal over inverted mass ordering; and a combined factor of 1.58 (0.63) for upper θ 23 octant and normal ordering.When calculating the BFs, the alternate hypothesis is any other combination of octant and mass ordering.Interpreting the largest BFs with the Jeffreys' scale, there is substantial evidence for the normal ordering when marginalising over the octant, and substantial evidence for the upper octant when marginalising over the mass ordering.In the more recent interpretation of BFs by Kass and Raftery [134], these both correspond to positive evidence.Importantly, the Jeffreys and Kass-Raftery definitions of "evidence" do not equate to the criteria often used in particle physics.For instance, a probability of 95.4% ("2σ ") is equivalent to a BF of 20.7, which is deemed as "decisive" on the Jeffreys' scale, and as "strong" on the Kass-Raftery scale.

8.1.2
The CP-violating phase δ CP , and sin 2 θ 13 A comparison of sin 2 θ 13 −δ CP contours with and without the reactor constraint is shown in Fig. 23.The regions are in good agreement, with a majority of the 1σ regions overlapping, comparable with the reactor constraint.A comparison of the δ CP posterior distributions is shown in Fig. 24, showing the Fig. 22 The marginalised posterior probability density of sin 2 θ 23 with (red) and without (blue) the reactor constraint on sin 2 θ 13 applied.The shaded areas show the 68% and 95% regions of highest posterior density, equivalent to the 1σ and 2σ credible intervals impact of the reactor constraint on T2K's δ CP result.The external constraint breaks the partially degenerate effects of sin 2 θ 13 and δ CP on the ν e appearance, leading to the ν-mode 1Re and 1Re1de selections having a larger sensitivity to δ CP .

The Jarlskog invariant
The sampled posterior probability density is in part a function of the PMNS mixing angles and δ CP , which means the probability distribution for the Jarlskog invariant [22,23], J = sin θ 13 cos 2 θ 13 sin θ 12 cos θ 12 sin θ 23 cos θ 23 sin δ CP (8) can be extracted directly from the steps in the MCMC.The posterior distribution for J is presented in Fig. 25, which favours a near-maximal negative J. The prior probability distribution is largely flat in the range J = [−0.035,0.035], with the fall-off beyond that coming from external θ 12 and θ 13 constraints.The preference for sin 2 θ 23 values near maximal mixing has the effect of picking out the more extreme values of J.When sampling the full posterior probability, which incorporates the δ CP constraint, a preference for negative values of J emerges.The blue curve in Fig. 25 is recreated in Fig. 26 showing the 1σ, 2σ, and 3σ credible intervals.Two-dimensional credible regions for the Jarlskog invariant against both sin 2 θ 23 and δ CP are included in Appendix C.
Although this analysis does not rule out CP-conserving values of δ CP at 2σ, it does rule out J = 0 at the 2σ level and excludes the J > 0.17 region at > 3σ with a flat prior probability in δ CP .The dependence of J on the choice of a prior flat in δ CP or flat in sin δ CP is shown in Fig. 26.The prior flat in sin δ CP flattens out the Jarlskog distribution, which Fig. 23 68% and 90% credible intervals from the marginalised sin 2 θ 13 − δ CP posterior distribution with (red) and without (blue) the reactor constraint (green band) applied, marginalised over both mass orderings in turn slightly expands the 2σ credible interval to where J = 0 is just included.These conclusions agree with previous studies on the impact of the δ CP prior at T2K [1,2].

Goodness-of-fit analysis
Predictions for the five samples at the FD are formed in Fig. 27, using the posterior probability distributions for the systematic uncertainties and oscillation parameters from the fit to data.By eye, the predictions agree well with the data, which are plotted as orange data points with statistical uncertainties applied.To quantify the model agreement with the data, the posterior predictive p-values [135] are calculated.These p-values can be calculated using either the total number of events per sample (rate-based) or the events per bin of each sample (shape-based).It can also be split by sample, or calculated as a total p-value.When including all samples, the shape-based and rate-based approach give p = 0.73 and p = 0.30 respectively.The p-values from both shape-and rate-based calculations broken down by sample and in total are tabulated in Table 12.Good p-values are demonstrated for all cases.
Fig. 24 The marginalised posterior probability density of δ CP with (red) and without (blue) the reactor constraint applied.The shaded areas show the 68% and 95% regions of highest posterior density, equivalent to the 1σ and 2σ credible intervals Fig. 25 Posterior probability distributions for the Jarlskog invariant using a prior distribution from the 2019 PDG reactor constraint on θ 13 [11] (red), prior from all parameters except sampling θ 23 from the T2K posterior (green), and the full T2K posterior (blue).All three posterior probabilities used a prior probability distribution flat in δ CP

Frequentist results
As in previous T2K analyses, the frequentist results are obtained using the marginal likelihood L marg (θ ) = dη p(η) L (θ, η) as the test statistic.Here, L (θ, η) is the binned Poisson likelihood for the parameter of interest, θ, and the nuisance parameters, η.The statistical treatment of nuisance parameters in the fit is thus identical to the Bayesian analy-Fig.26 Posterior probability distributions for the Jarlskog invariant taken from posterior distributions with priors that are either flat in δ CP (blue) or flat in sin δ CP (orange).1σ, 2σ, and 3σ credible intervals are shown as the region between the vertical black solid line and the specified vertical dashed lines sis and assumes a prior probability distribution p(η).The numerical integration is performed by varying systematic uncertainties with a Gaussian covariance matrix from the ND analysis in Sect.6 as a constraint, and varying the other oscillation parameters with a flat prior probability distribution on sin 2 θ 23 , δ CP , Δm 2 32 , and sin 2 2θ 13 , or a Gaussian prior on sin 2 2θ 13 .Confidence intervals and regions are constructed with two different methods.For critical parameters with known boundary effects, the Feldman-Cousins (FC) method [136] is utilised to calculate the coverage.This is performed for the result using the reactor constraint, on the one-dimensional confidence intervals in δ CP and sin 2 θ 23 , and their joint confidence region.For generating the ensemble of experiments for FC evaluation, the nuisance oscillation parameters are varied from the posterior distribution obtained by fitting a representative simulated data set, sometimes referred to as "Asimov data".This simulated data set is generated at the global best-fit point using the reactor constraint.Since the FC method is computationally intensive, the remaining confidence regions are constructed using constant Δχ 2 (θ ) = χ 2 (θ ) − min θ χ 2 (θ ) values via Wilks's theorem [137], where χ 2 = −2 ln L marg .Whether Δχ 2 is computed with respect to the minimum over both mass orderings, or the minimum in each mass ordering separately, is indicated in each of the results from the frequentist Fig. 27 The reconstructed neutrino energy distributions of each FD sample.Data with Poisson uncertainties are shown in orange and the distributions of the predictions are shown in the coloured background, with the mean of those distributions overlaid in red.The z-axis represents the number of MCMC samples that had a prediction in a specific bin, and its intensity is directly proportional to the probability.The predictions are built by sampling both the nuisance and oscillation parameters from the posterior probability distribution in the Bayesian analysis.sin 2 θ 13 is constrained from T2K data alone with no reactor constraint applied analysis.The frequentist analysis bins the e-like FD samples in reconstructed lepton angle and reconstructed lepton momentum, and the μ-like samples in reconstructed lepton angle and the reconstructed neutrino energy, defined in the same way as in the Bayesian analysis presented in Sect.8.1.
In previous analyses, the μ-like samples were binned only in reconstructed neutrino energy, and adding the lepton angle information increases the 1σ expected sensitivity to Δm 2

by O(1%).
Global best-fit values are given in Table 13.As noted in the Bayesian section, the results with and without the reactor constraint are compatible, with the former resulting in stronger constraints on δ CP and sin 2 θ 23 .All the following results are from the fit to data using the reactor constraint.

The CP-violating phase δ CP , and mass ordering
Figure 28 shows the Δχ 2 distributions for δ CP in both mass orderings with FC-adjusted confidence intervals, which are also summarised in Table 14.A large region of sin δ CP > 0 is excluded at > 3σ confidence level (CL), whereas the CPconserving values δ CP = 0, π are excluded at 90% CL.In particular, δ CP = π is just inside the 2σ interval.
As was also seen in the Bayesian analysis, the constraint on δ CP is weaker compared to T2K's previous analysis [1,2]. Figure 29 shows the impact on the Δχ 2 of δ CP after each Fig. 28 The Δχ 2 distribution in δ CP from fitting to the data with the reactor constraint applied.The confidence intervals in the shaded regions are calculated using the FC method Fig. 29 The Δχ 2 distribution in δ CP for incremental modifications from the previous analysis [1] to this result, for normal and inverted mass orderings."E" corresponds to this analysis, except that unlike the main frequentist result, the μ-like samples do not use the scattering angle information for better compatibility with the previous T2K analysis update introduced in this analysis, all of which weaken the δ CP constraint, with the largest contribution being the addition of the latest data in run 10 at the FD.The data is now more consistent with the expectation, shown in Fig. 30 for both normal and inverted ordering.In most of the δ CP parameter Table 13 Results for the oscillation parameters from the fit to data with and without the reactor constraint in the frequentist analysis, with the confidence intervals estimated using the constant Δχ   | in inverted mass ordering) for the data fit with the reactor constraint applied, obtained with the constant Δχ 2 method, where in each mass ordering hypothesis a fixed mass ordering is assumed space, T2K is below the upper limit of the 68% expectation band of ensemble experiments at maximal CP violation.The inverted mass ordering is disfavoured at more than 1σ for all δ CP values, mostly consistent with the expected sensitivity at sin δ CP = −1.Replacing each sample's event distribution in data by the expectation of the model shows that the stronger δ CP constraint observed in data compared to the expectation comes from the ν-mode 1Re1de sample.

Atmospheric oscillation parameters
The sin 2 θ 23 − Δm 2  32 confidence intervals are presented in Fig. 31.The contours are compatible for both mass orderings, with a slight shape change compared to the previous T2K analysis, to now marginally prefer the upper octant.Figure 32 shows the contour from a fit using only the 1Rμ samples, which shows that the constraint is dominated by the 1Rμ samples, with the 1Re samples providing the sensitivity to the octant of θ 23 .
The evolution of the sin 2 θ 23 −Δm 2 32 contour from the fit to data after introducing each update in the analysis is shown in Fig. 33.The most significant impact on the Δm 2  32 constraint comes from changing the cross-section model and updating the ND constraint.For Δm 2  32 , improvements in the removal energy uncertainty have significantly reduced the uncertainty before the smearing based on simulated data studies has been applied.Thanks to increased robustness of the uncertainty model, the size of the smearing has also been reduced by a factor of 2.8, discussed in detail in Sect.9.For sin 2 θ 23 , there is a slight shift in shape from the latest data in the FD.The new data favour a slightly larger sin 2 θ 23 , and so pushes less against the boundary of maximal mixing.This results in a slightly weaker constraint in the lower octant and similar constraint in the upper octant compared to the previous analysis.32 for normal ordering, between a full fit and a fit using only the μ-like samples.The intervals are calculated with the constant Δχ 2 method, and applying the reactor constraint Fig. 33 The Δχ 2 distribution in Δm 2  32 and sin 2 θ 23 for normal ordering, showing incremental modifications of the previous analysis [1,2] to this result."E" corresponds to this analysis, except that unlike the main frequentist result, the μ-like samples do not use the scattering angle information for better compatibility with the previous analysis.The best-fit point for "C" in orange is the same as "D" in green The Δχ 2 distribution for sin 2 θ 23 is shown in Fig. 34 with the confidence intervals summarised in Table 14.The new data at the FD has reduced compatibility with maximal mixing, which is now outside the FC-corrected 1σ confidence interval.Whilst the upper octant is favoured at 1σ CL, the data is still compatible with both octant hypotheses at 90% CL.These results are compatible with the sensitivity, shown in Fig. 35.

Cross-fitter comparisons
To compare the consistency of the Bayesian analysis described in Sect.8.1 and the frequentist analysis described in Sect.8.2, the Bayesian posterior distributions were recast into frequentist Δχ 2 distributions comparable to the frequentist analy-Fig.34 The Δχ 2 distribution in sin 2 θ 23 for fitting to the data with the reactor constraint applied.The confidence intervals in the shaded regions are calculated using the FC method Fig. 35 The Δχ 2 distribution in sin 2 θ 23 for fitting to the data with the reactor constraint applied, overlaid with the expected distributions from an ensemble of simulated experiments created with true NO and sin 2 θ 23 = 0.56 sis. Figure 36 shows comparisons of sin 2 θ 23 − Δm 2  32 and sin 2 θ 13 − δ CP contours from fits to data from both frameworks.The minor differences between the resulting contours can be attributed to two distinct analysis choices: the way in which the constraints from the near-detector analysis on the systematic uncertainties are applied, and the choice of kinematic variables in which the far detector samples are binned.The Bayesian analysis uses a ND analysis with irregular rectangular binning for the ND samples to better adapt to differences in p μ − cos θ μ phase space density, whereas the frequentist analysis' ND constraint uses regular binning.Both analyses use the lepton scattering angle for the e-like samples, but differ in the use of reconstructed energy (Bayesian) or reconstructed lepton momentum (frequentist) for the μ-like samples.Additionally, the frequentist analysis also uses the reconstructed lepton scattering angle for the μlike samples to disentangle systematic uncertainties related to energy scale and mis-reconstructed backgrounds at the FD.Fig. 36 Comparison of the 68% and 90% confidence intervals from fits to data from the Bayesian analysis ("Analysis A") and the frequentist analysis ("Analysis B"), discussed in Sect.8.3."Analysis B, A-like" configures the frequentist analysis in the same way as the Bayesian analysis, using the same binning at the FD and the same MCMC-based ND analysis.The contours are extracted from fits that fix the neutrino mass ordering to the normal ordering and apply the reactor constraint on sin 2 θ 13 Other differences, like the non-Gaussian nature of parameters included in the ND constraint or event-by-event vs. binned oscillation probability calculation, had little effect.When the frequentist analysis is configured to impose the constraints from the ND MCMC analysis and bin the FD samples similarly to the Bayesian analysis, these minor differences abate, shown in Fig. 36.The uncertainty models of both analyses were validated against each other for consistency and were found to agree.

Comparisons with other experiments
In the global context of neutrino oscillation experiments, these results provide leading constraints on both the atmospheric oscillation parameters, Δm 2  32 and sin 2 θ 23 , and the CP-violating phase, δ CP .Whereas the other experiments profile over parameters to calculate the Δχ 2 , T2K instead calculates the marginal likelihood for the Δχ 2 .Figure 37 shows the 90% confidence regions in sin 2 θ 23 −Δm 2  32 for the normal ordering from the frequentist analysis, compared to NOvA, SK and IceCube.There is general agreement between the experiments, with T2K providing the strongest constraints on both parameters.Figure 38 compares the 90% confidence regions in sin 2 θ 23 − δ CP for both orderings to NOvA and SK.The confidence intervals on sin 2 θ 23 significantly overlap, as do the intervals for δ CP .In the normal ordering, T2K excludes large regions of the NOvA constraint at 90% confidence interval, and NOvA excludes parts of T2K's 90% confidence interval.In the inverted ordering, the experiments consistently favour the π < δ CP < 2π region, with a weak for normal ordering with NOvA [138], Super-K [139], IceCube [140], and MINOS+ [141].The NOvA and IceCube constraints obtained with the FC method, but with different treatment of the mass ordering: NOvA takes the minimum over both mass orderings, whereas the IceCube contours assume normal ordering.The T2K, Super-K, and MINOS+ contours are computed with the constant Δχ 2 method, assuming normal ordering preference for the upper octant.Importantly, there is no significant tension between the experiments, and more data is needed to elucidate the matter.Furthermore, the joint oscillation analyses with the NOvA and SK collaborations will help address this.

Simulated data studies
Simulated data studies with the frequentist analysis were used to investigate the impact of alternative model predictions and data-driven tunes, discussed in Sect.5.3, on the oscillation parameter constraints.The oscillation analysis in Sect.8 had these uncertainty inflation strategies applied, and this section summarises the procedure, with details provided in Appendix B.

Methodology
In the simulated data studies, the prediction from an alternative model is treated as the data at the ND, and is fit with the usual systematic uncertainty model.The parameters are fit to simulated data at the ND and are propagated to the FD, and the reconstructed energy spectrum and oscillation parameter constraints are compared to an "Asimov" data set.In an Asimov analysis, the parameters for the systematic uncertainties are set to specific values and the predicted spectra at each detector is treated as the data, giving an expectation of the sensitivity if no statistical fluctuations were present.For the oscillation parameters, two separate Asimov points were tested: one close to T2K's best-fit point, and one with δ CP = 0 and non-maximal sin 2 θ 23 , detailed in Appendix B, Table 17.This section only presents results with the Asimov data near T2K's best-fit point.The PDG 2019 reactor constraint on sin 2 θ 13 is applied in the following studies, but had little impact on the overall conclusions.Although simulated data sets can result in both weaker and stronger constraints on the oscillation parameters than the expected sensitivity, they are only used to inflate the uncertainties in this analysis.
The simulated data set procedure mainly identifies two types of effects: • Systematic uncertainty model shortcomings: If the systematic uncertainty model is or if the effect of the alternative model is small, the oscillation parameter contours obtained with the simulated data sets will not see a bias with respect to the expected sensitivity.The bias is quantified as the percentage change of the middle of the 1σ confidence interval of an oscillation parameter, relative to the 1σ from systematic uncertainties in the expected sensitivity analysis.An example is discussed in Appendix B.1.• ND to FD extrapolation issues: Some alternate models may not produce a significant bias the oscillation parameters, often due to the low sensitivity of the samples they affect.Issues in the extrapolation process can be exposed by comparing three distributions: (i) the predicted spectrum at the FD from fitting the Asimov data set at the ND, (ii) the predicted spectrum at the FD from fitting to the simulated data at the ND, (iii) and the pre-dicted spectrum at the FD when applying the alternative model directly.Even though the bias on the oscillation parameters at T2K statistics may be small, simulated data studies may guide which of the systematic uncertainties are important to address in future T2K analyses and upcoming high-statistics experiments, such as Hyper-Kamiokande [142] and DUNE [143].An example is discussed in Appendix B.2.
All the individual biases are summed in quadrature and are used to inflate the confidence interval for Δm 2  32 , due to its simple Gaussian probability density.For δ CP , the effect of systematic uncertainties is much smaller and the probability density is non-Gaussian, so a different method is applied.Each simulated data set is studied to see if it impacts any major claims in the analysis; in this analysis the 90% confidence interval of δ CP .This is done by calculating the difference in the Δχ 2 distribution for δ CP for the Asimov data and the simulated data, and adding the difference to the Δχ 2 distribution for δ CP from the data, where the 90% confidence interval was calculated using the FC method mentioned in Sect.8.2.
Simulated data sets can drastically increase or decrease the number of events at both detectors.In such cases, comparing the constraints on the oscillation parameters to the expected sensitivity conflates the effects of propagating mismodelling from the ND analysis with the impact of increased or decreased statistics from the simulated data set.For instance, an alternative model that increases the number of ( ν ) e at the FD near the oscillation maximum will likely lead to a stronger constraint on δ CP due to the measurement being dominated by statistical uncertainties in those samples.To gauge this effect, the three predictions from the ND to FD extrapolation studies, outlined earlier, are used.If the model from the ND simulated data analysis predicts the spectrum of the alternative model well at the ND, and correctly predicts the spectrum at the FD compared to when directly applying the alternative model, a "scaled Asimov" approach is utilised.In these cases, two changes to the procedure are made: (a) the propagated ND constraint is the expected sensitivity to the systematic parameters if the real data is as predicted by the pre-fit model, and (b) the variation to the model that was used to build the simulated data set is also applied to the simulation that is being fit at the FD.This removes most of the statistical effect and better captures the features of propagating a mismodelling in the ND analysis.For this analysis, the scaled approach was only used for the 2p2h Martini simulated data set.
The simulated data studies all concerned the interaction model and were detailed in Sect.38 Comparison of 90% confidence regions in sin 2 θ 23 − δ CP over both mass orderings with NOvA [138] and Super-K [139].The T2K and NOvA confidence regions have been computed using the FC method, whereas the Super-K results are obtained with the constant Δχ 2 method

Results
Table 15 summarises the observed biases on the oscillation parameters, showing the simulated data set with the highest impact from each category.The full results are shown in Appendix B, Table 18.The impact of the simulated data studies on sin 2 θ 23 and δ CP was found to be small compared to the impact of statistical and systematic uncertainties.The largest bias on Δm 2 32 relative to the systematic uncertainty was found to be 57.8% from the pion SI simulated data set, and 20.8% relative to the overall uncertainty.Selected simulated data studies were added in quadrature 4 to avoid double counting similar physics effects, leading to an overall smearing on Δm 2 32 of 1.35 × 10 −5 eV 2 .For comparison, the overall uncertainty on Δm 2 32 from the expected sensitivity study, before the simulated data procedure, was 5.7 × 10 −5 eV 2 and is dominated by the uncertainty from statistics, which was 5.3 × 10 −5 eV 2 .Generally, the simulated data studies had a smaller impact on δ CP , due to its uncertainty being dominated by the statistics in the electron-like selections at the FD.
In the previous T2K analyses [1,2], the simulated data study for the nucleon removal energy had a significant impact, especially on Δm 2  32 , and an additional uncertainty was introduced.In this analysis, the updated nucleon removal energy uncertainty, described in Sect.5, has caused it to no longer be the dominant source of systematic uncertainty.
Table 16 shows the changes to the 90% confidence interval for δ CP for each of the simulated data studies.The non-CCQE and the data-driven pion momentum modification simulated data sets had the largest impact, shifting the 90% CL by 0.09 and 0.07 respectively.The change to the 90% confidence limits does not alter the conclusions on the exclusion of CPconserving values presented in Sect.8.

Conclusions
The T2K collaboration has measured the three-flavour PMNS neutrino oscillation parameters Δm 2 32 , sin 2 θ 13 , sin 2 θ 23 , δ CP , the Jarlskog invariant J, and the mass ordering, using the  20 POT more (anti-)neutrino data at the ND.For the first time, a neutrino flux constraint using charged pion data from a T2K replica target at NA61/SHINE was used, which approximately halves the flux uncertainty before the ND analysis.An updated neutrino interaction model with a refined initial-state and removal-energy model with associated uncertainties was also employed, amongst others.High statistics ND data was used to constrain the neutrino flux and interaction model uncertainties at the FD, which also constrains the wrong-sign background of the neutrino beam with the magnetised ND.Biases from unmodelled systematic uncertainties were studied through simulated data studies, which acted to inflate the Δm 2  32 and δ CP confidence intervals.These results present the strongest constraints on several neutrino oscillation parameters, and are more consistent with the expected sensitivities to the oscillation parameters compared to T2K's previous analysis.
The results are limited by statistics, and T2K will continue to take data as J-PARC upgrades [27,28] the neutrino beam for the Hyper-Kamiokande experiment [142].In preparation, the T2K beamline has recently undergone a long shutdown, and will be operating the magnetic horns at 320 kA current, with beam power in excess of 700 kW in the near future [27,28].Upcoming analyses at the FD will expand selections to include multiple Cherenkov rings, increasing statistics by approximately 30%.The FD has also begun collecting data with gadolinium doped in the ultrapure water [39], drastically increasing the efficiency in tag-ging interactions producing neutrons.At the ND, selections are being developed to improve the understanding of nuclear effects in neutrino interactions, such as 2p2h, in-medium corrections, and the initial state, thus addressing the larger systematic uncertainties in this analysis.Furthermore, the ND280 upgrade [144][145][146] will be ready to take data in 2023, providing significantly improved reconstruction capabilities for low momentum protons and pions with full angular coverage, which will allow for detailed study of nuclear effects, in addition to measurements of neutron kinematics.Moreover, the T2K collaboration is actively working with the NOvA and SK collaborations on combined neutrino oscillation analyses, taking advantage of synergies in experiment design to lift current degeneracies and increase statistics.

Data Availability Statement
This manuscript has data included as electronic supplementary material.The online version of this article contains supplementary material, which is available to authorized users.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.The images or other third party material in this article are included in the article's Creative Commons licence, unless indi-Fig.39 The ND prediction after the fit to data for the FGD1 ν μ CC0π sample is shown with the data, where the non-CCQE contribution is shown in shaded green.The dotted red line shows the non-CCQE contribution needed for the overall prediction to match the data.The weight is the ratio of the dotted red line to the shaded green B.1: Non-CCQE simulated data set As described in Sect.5, the Q 2 normalisation parameters included in this analysis are introduced to provide ad hoc freedom in the NEUT CCQE model, based on known discrepancies with data [81][82][83][84].The robustness of this parametrisation is tested and its impact on oscillation parameter measurements by assuming that the underestimation of data is entirely due to non-CCQE interactions.The simulated data set to replicate this effect is devised by • The prediction after the ND fit to data is weighted to restore the CCQE Q 2 parameters to their nominal values (i.e.unity).• The reconstructed Q 2 rec distributions are built for data and the simulation from the previous step using Eq. 3, with bins matching the ranges of the Q 2 parameters.Since the study concerns CC0π interactions, this is only done for the CC0π samples.
• The simulated non-CCQE events with true CC0π topology are varied until the prediction matches the data .• The change in the non-CCQE true CC0π events is then applied as a set of weights to the nominal prediction as a function of true Q 2 rec , resulting in the simulated data set.
This process is applied separately to ν μ and ν μ samples to extract different weights for neutrino and anti-neutrino events, where the effect was smaller for the ν μ samples.The weight extraction process and the resulting weights for neutrino events are shown in Fig. 39.
The resulting distribution was used as the simulated data in the ND fit.As shown in Fig. 40, the neutrino flux, 2p2h ν normalisation, and the Q 2 parameters shift to accommodate Fig. 40 Constraints on the ν-mode ν μ flux (top) and CC0π crosssection parameters (middle, bottom) from the fit to the non-CCQE simulated data set (black points, black lines), overlaid on the input uncertainty (red band).The parameters are presented as a ratio to the generated value in NEUT, except for the removal energy parameters which show the shifts in units of MeV the model variation, notably present in the higher E ν region.This is largely expected, since non-CCQE events are often from higher energy neutrinos due to the interaction cross section increasing with E ν , whereas CCQE plateaus.
The same set of weights were used to create a corresponding simulated data set for the FD samples.The ND constraints from this simulated data study were propagated to the FD simulation, and used to extract oscillation parameters.Figure 41 shows the oscillation parameter contours for an expected sensitivity fit and the "non-CCQE" simulated data set.Quantitatively, the simulated data set changes the interval width of 32.7% on Δm 2 32 and 21.3% on sin 2 θ 23 with respect to the size of the systematic uncertainty, σ syst .This is equivalent to an uncertainty inflation of 0.69 × 10 −5 eV 2 .This simulated data set has an increased sensitivity to δ CP , and no significant bias in the best-fit value is observed.
This study attempts to address the differences between neutrinos and anti-neutrinos by applying different weights for the two.Another study may also take into account the differences between interactions on carbon and oxygen targets, which some models predict should have different behaviour at low Q 2 .However, the ability for the ND analysis to distinguish behaviour at low Q 2 on carbon and oxygen is limited by statistics, and such a detailed study is currently not warranted.
Fig. 42 Ratio between the prediction after the fit to the simulated data from the MINERvA pion suppression to the simulated data itself, for the ND FGD1 CC1π sample Fig. 43 The predictions for the 1Re1de selection at the FD using the unchanged model (green) and including the MINERvA pionsuppression (blue).The red bands show the propagated ND prediction from the fit to the simulated data set, with the 1σ uncertainty bands B.2: MINERvA single-pion suppression Studies of MINERvA's neutrino-induced single-pion production (SPP) data found that the GENIE event generator over-predicted the cross section [125].The effect was particularly pronounced at low 4-momentum transfer, where a suppression of the cross section by ∼ 60% was needed to match the data.The data-driven tune from the paper was applied to T2K SPP events to gauge the impact of this low-Q 2 modification on the oscillation analysis.Importantly, the SPP model used in the GENIE event generator by MINERvA differs from that in NEUT for this analysis, especially in the low-Q 2 region, and T2K SPP data [148] sees little need for a suppression with recent NEUT versions [149].Furthermore, MINERvA sees the largest suppression in CC1π 0 selections, which are barely selected in T2K's FD.The CC1π + selec-Table 18 Biases on the main oscillation parameters for each simulated data set, calculated as the shift in the middle of the 1σ confidence interval relative to the overall uncertainty from systematic sources ("Syst.")and the total ("Total") to one decimal place

Simulated data set
Relative to sin 2 θ 23 (%) Δm tion, which has a larger overlap with selections in T2K's FD, has a significantly weaker suppression in the publication.The low-Q 2 suppression should therefore be considered as a conservative variation, designed to build an extreme simulated data study to analyse "worst-case scenario" mismodelling.
After the fit to the simulated data, the majority of affected parameters are related to SPP, with little change to the neutrino flux parameters and the CC0π interaction parameters.This is because they are largely constrained by the dominant CC0π selections at the ND, which are barely affected by the simulated data.C A 5 was pulled down, which acts to decrease the overall SPP cross section, and M R E S A is pulled up, which increases the high Q 2 SPP cross section.Since CC coherent events primarily occupy the low-Q 2 region, they too are suppressed by the CC coherent normalisation parameter, which is pulled lower.The 2p2h normalisation and highest Q 2 CCQE normalisations are also pulled low, as such events occupy a similar space in p μ − cos θ μ to SPP events in the CC0π selection.The ratio of the prediction versus the simulated data at the ND is shown in Fig. 42, where the CC1π sample is over-predicted by 20% in certain p μ − cos θ μ regions.Hence, the current uncertainty model at T2K can not account for a sizeable low-Q 2 suppression of SPP events.Importantly, this has a low impact on the dominant ND CC0π selections, since SPP is a sub-leading contributor.
The prediction for the SPP-dominated 1Re1de selection is shown in Fig. 43 and echoes the ND prediction; the prediction from fitting the ND model to the simulated data can not replicate the alternative model.However, the impact on the oscillation parameters is relatively small as the 1Re1de sample has low statistics and occupies higher E ν compared to the dominant ν-mode 1Re selection, which provides most of the ν e appearance at T2K.This issue may present significant biases for experiments with higher statistics of SPP, such as NOvA, DUNE and Hyper-Kamiokande.It highlights the need for high-quality single-pion neutrino data, and the development of sufficiently robust and flexible models to describe them.

B.3: Summary of biases
As detailed in Sect.9, a selection of the largest contributors to the bias were used to inflate the confidence intervals for the two oscillation parameters, Δm 2 and δ CP .The down-selection was used to avoid double-counting similar physics effectsfor instance the CCQE axial form factor described by either the 3-component or z-expansion models -by including the simulated data set with the largest bias.This section gives the details of the biases for all the simulated data sets on the oscillation parameters.
Table 18 presents the full bias table for the three oscillation parameters sin 2 θ 23 , Δm 2  32 , and δ CP .As with previous T2K analyses, δ CP is significantly less sensitive to systematic uncertainties, owing to the relative dominance of statistical uncertainty in the single ring electron-like samples at the FD.The shifts in the 90% confidence interval of δ CP for each simulated data study is shown in Table 19, where the shifts are dominated by the three largest simulated data studies: the pion SI, the non-CCQE, and data-driven pion tune.

Appendix C: Two-dimensional Jarlskog credible regions
The measurements of δ CP and sin 2 θ 23 are the main constraints on the Jarlskog invariant from T2K's analysis.Both of these parameter measurements have the effect of preferring the more extreme values of the Jarlskog invariant that are otherwise allowed, as shown earlier in

Fig. 5
Fig.5 The predicted unoscillated neutrino fluxes at the FD in ν-mode (top) and ν-mode (bottom) in logarithmic scale with an extended E ν range, after the tuning to NA61/SHINE data on the T2K replica target

Fig. 10
Fig.10 Constraints on the ν-mode ν μ flux uncertainty parameters at the FD from the fit to ND data (black points, black lines), overlaid on the input uncertainty (red band)

Fig. 15
Fig.15 Reconstruction performance at the FD of stopping cosmicray muons and the Michel electrons from their decays.The left panel shows the reconstructed momentum distribution of those electrons for data taken during the SK-IV (blue) and SK-V (red) detector periods.

Fig. 20
Fig.20 Marginalised posterior probability densities from the Bayesian analysis for oscillation parameters of interest from a fit to data with the reactor constraint on sin 2 θ 13 applied.The two-dimensional posteriors have 68% (dashed) and 90% (solid) credible levels indicated and the point with highest posterior probability.The one-dimensional posteriors

Fig. 30 TheFig. 31
Fig.30 The Δχ 2 distribution in δ CP from fitting to the data assuming normal (left) and inverted (right) neutrino mass ordering, with the reactor constraint applied.The distribution is overlaid with the expectations

Fig. 32
Fig.32 Comparison of confidence regions in sin 2 θ 23 − Δm2  32 for normal ordering, between a full fit and a fit using only the μ-like samples.The intervals are calculated with the constant Δχ 2 method, and applying the reactor constraint

5 . 3 .
Details of the simulated data study procedure and two examples are provided in Appendix B.

Fig.
Fig.38 Comparison of 90% confidence regions in sin 2 θ 23 − δ CP over both mass orderings with NOvA[138] and Super-K[139].The T2K and NOvA confidence regions have been computed using the FC method, whereas the Super-K results are obtained with the constant Δχ 2 method CEA and CNRS/IN2P3, France; the DFG (RO 3625/2), Germany; the INFN, Italy; the Ministry of Education and Science (DIR/WK/2017/05) and the National Science Centre (UMO-2018/30/E/ST2/00441), Poland; the RSF19-12-00325, RSF22-12-00358 and the Ministry of Science and Higher Education (075-15-2020-778), Russia; MICINN (SEV-2016-0588, PID2019-107564GB-I00, PGC2018-099388-BI00, PID2020-114687GB-I00) Government of Andalucia (FQM160, SOMM17/6105/ UGR) and the University of Tokyo ICRR's Inter-University Research Program FY2023 Ref. J1, and ERDF funds and CERCA program, Spain; the SNSF and SERI (200021_185012, 200020_188533, 20FL21_ 186178I), Switzerland; the STFC and UKRI, UK; and the DOE, USA.We also thank CERN for the UA1/NOMAD magnet, DESY for the HERA-B magnet mover system, the BC DRI Group, Prairie DRI Group, ACENET, SciNet, and CalculQuebec consortia in the Digital Research Alliance of Canada, GridPP and the Emerald High Performance Computing facility in the United Kingdom, and the CNRS/IN2P3 Computing Center in France.In addition, the participation of individual researchers and institutions has been further supported by funds from the ERC (FP7), "la Caixa" Foundation (ID 100010434, fellowship code LCF/BQ/IN17/11620050), the European Union's Horizon 2020 Research and Innovation Programme under the Marie Sklodowska-Curie grant agreement numbers 713673 and 754496, and H2020 grant numbers RISE-GA822070-JENNIFER2 2020 and RISE-GA872549-SK2HK; the JSPS, Japan; the Royal Society, UK; French ANR grant number ANR-19-CE31-0001; the SNF Eccellenza grant number PCEFP2_203261; and the DOE Early Career programme, USA.For the purposes of open access, the authors have applied a Creative Commons Attribution licence to any Author Accepted Manuscript version arising.

Fig. 41
Fig. 41 Comparison of confidence regions in sin 2 θ 23 − Δm 2 32 (top) for the expected sensitivity fit (blue) and the non-CCQE simulated data study (red), and the Δχ 2 as a function of δ CP (bottom)

Fig. 46 FC
Fig. 46 FC-corrected critical Δχ 2 values for 1D-fits in δ CP and sin 2 θ 23 , computed at true values indicated with black dots and linearly interpolated in between.The error bands show the 1σ toy-statistical uncertainty (binomial confidence interval).The dotted lines show the asymptotic values from Wilks's theorem for reference

Fig.
Fig. 68% and 90% FC-corrected critical Δχ 2 values for 2D-fits of sin 2 θ 23 − δ CP over both mass orderings, computed at the points indicated by circles and bi-linearly interpolated for the regions in between the points.The colour scheme is chosen to show the asymptotic values from Wilks's theorem in white

Fig. 25 .
T2K's preference for values of θ 23 ∼ 45 • maximises the factor of sin θ 23 cos θ 23 .Its best-fit value being slightly above that mark allows for smaller values of J as illustrated in Fig. 44.Similarly, T2K's preference for values of δ CP ≈ −90 • maximises the sin δ CP factor and selects negative values of J. Figure 45 illustrates the influence of the δ CP posterior with values of δ CP further away from maximal CP-violation mapping to smaller preferred magnitudes of J.

Table 1
Collected protons-on-target (POT) for each T2K run included in the analysis of T2K data at the ND and FD.The recorded POT at INGRID closely follows that of the FD

Table 2
Percentage of hadronic interactions in the target and downstream beam line for which external measurements are used in the tuning or uncertainty evaluation.The interactions are weighted by their contribution to the neutrino flux at the FD, separated into different horn focusing modes and neutrino flavours SK-V resumed data taking in January 2019 with ultrapure water and collected 4.73 × 10 20 POT during October 2019-February 2020 (run 10).These data were collected entirely in ν-mode, resulting in a total of 19.66 × 10 20 and 16.34 × 10

Table 5
Uncertainties on the total number of events in the ND analysis from detector uncertainties only, broken down by selection

Table 6
Number of events in each of the ND selections for data and the ratio to the prediction before and after the fit to data

Table 7 p
-values comparing the variations of the model before the ND analysis and the model fit to the data, broken down by likelihood contributors, and showing the p-value for all samples, and the total p-value including all samples and all systematic uncertainties

Table 9
Predictions for the number of events at the FD using oscillation parameters and systematic uncertainty parameters at their best-fit values whilst varying δ CP

Table 11
Fractions of posterior probability in different combinations of the mass ordering and θ 23 octant from fit to T2K data with (without) the reactor constraint on sin 2 θ 13 .NO (IO) refers to the normal (inverted)

Table 12
Breakdown of posterior predictive p-values by sample, quoted separately using a shape or rate based calculation, demonstrating good compatibility between the model and the data

Table 14
FC-corrected confidence intervals for δ CP and sin 2 θ 23 from the fit to data in the frequentist analysis, using the reactor constraint on sin 2 2θ 13 .The 3σ FC correction was not computed for sin 2 θ 23

Table 15
Biases on the main oscillation parameters for each simulated data set, calculated as the shift in the middle of the 1σ confidence interval relative to the overall uncertainty from systematic sources ("Syst.")and the total ("Total") to one decimal place

Table 16
Shifts of the 90% confidence interval boundaries of δ CP , in radians, as a result of the simulated data studies.The values in the top row correspond to the results of the data fit, assuming normal ordering.The values for each simulated data set are added to (subtracted from) the right (left) δ CP interval edge from the data fit.Only the absolute size of the shift is taken into accountThe simulated data sets with the largest impact are typed in bold statistics at the FD equivalent to 3.6 × 10 21 POT.T2K continues to favour neutrino oscillations with near-maximal CP violation, in the upper octant of sin 2 θ 23 , in the normal mass ordering, with a sin 2 θ 13 consistent with the measurements by reactor experiments.The analysis included 4.72 × 10 20 POT more neutrino data at the FD, and 5.73(4.48)× 10

Table 19
Shifts of the 90% confidence interval boundaries of δ CP , in radians, as a result of the simulated data studies.The values in the top row correspond to the results of the data fit, assuming normal ordering.The values for each simulated data set are added to (subtracted from) the right (left) δ CP interval edge from the data fit.Only the absolute size of the shift is taken into account The simulated data sets with the largest impact are typed in bold

Table 20
FD ν-mode flux parameters before and after the fit to ND data, including the uncertainty from the diagonal of the covariance matrix.The values in brackets show the range of E ν in units of GeV for each parameter

Table 21
FD ν-mode flux parameters before and after the fit to ND data, including the uncertainty from the diagonal of the covariance matrix.The values in brackets show the range of E ν in units of GeV for each parameter

Table 22
Cross-section parameters before and after the fit to ND data, including the uncertainty from the diagonal of the covariance matrix.The parameters are detailed in Sect. 5. Parameters without external constraints are labelled with an uncertainty ±∞