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}\theta_{ij}$, $\Delta m^{2}_{ij}$, and
$\delta_{\scriptscriptstyle\mathrm{CP}}$. T2K is not sensitive to the “solar”
oscillation parameters $\sin^{2}\theta_{12}$ and $\Delta{}m^{2}_{21}$,
therefore constraints from the world averages reported in PDG 2019
Tanabashi:2018oca  are imposed222$\sin^{2}\theta_{12}=0.307(\pm 0.013)$,
$\Delta{}m^{2}_{21}=7.53(\pm 0.18)\times 10^{-5}\leavevmode\nobreak\
\mathrm{eV^{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}\theta_{13}$ from the world average reported in PDG
2019 Tanabashi:2018oca , referred to as the “reactor
constraint”333$\sin^{2}\theta_{13}=2.18(\pm 0.07)\times 10^{-2}$. The reactor
constraint has a significant effect on the sensitivity to other oscillation
parameters of interest, notably $\delta_{\scriptscriptstyle\mathrm{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.

(a) $\nu$-mode 1R$\mu$

(b) $\nu$-mode 1R$e$

(c) $\nu$-mode 1R$e$1d$e$

(d) $\overline{\nu}$-mode 1R$\mu$

(e) $\overline{\nu}$-mode 1R$e$

Figure 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 App. B,
Tab. 17.

###  8.1 Bayesian results

Figure 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}\theta_{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 have
68%, 90%, and 95% credible intervals indicated in different shades of grey.
All credible regions are calculated from marginalising over both mass
orderings, although panels displaying $\Delta{}m^{2}_{32}$ show only the
portions of the distributions in the normal mass ordering
($\Delta{}m^{2}_{32}>0$).

The Bayesian results presented in this section are obtained by sampling the
posterior distributions through MCMC metropolis ; hastings  analysis, using
the ND and FD selections simultaneously. The MCMC analysis presented in Sec. 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). For the 1R$e$1d$e$ selection—which is dominated by
$1e^{-}1\pi^{+}$ final states—the nucleon mass is replaced by the
$\Delta(1232)$ mass. The $\mu$-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}\theta_{23}$,
$\Delta{}m^{2}_{32}$, $\delta_{\scriptscriptstyle\mathrm{CP}}$ (or
$\sin\delta_{\scriptscriptstyle\mathrm{CP}}$), and Gaussian priors are applied
on $\Delta{}m^{2}_{21}$ and $\sin^{2}\theta_{12}$. For $\sin^{2}\theta_{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%.

Fig. 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.

####  8.1.1 Atmospheric oscillation parameters

Figure 21: 68% and 90% credible intervals from the marginalised
$\sin^{2}\theta_{23}-\Delta{}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.

The effects of applying the reactor constraint on the
$\sin^{2}\theta_{23}-\Delta{}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}\theta_{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.

|  | $\sin^{2}\theta_{23}$ | Sum
---|---|---|---
|  | $<0.5$ | $>0.5$
$\Delta{}m^{2}_{32}$ |  $>0$ (NO) | 0.195 (0.260) | 0.613 (0.387) | 0.808 (0.647)
$<0$ (IO) | 0.035 (0.152) | 0.157 (0.201) | 0.192 (0.353)
Sum | 0.230 (0.412) | 0.770 (0.588) | 1.000

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

The results for the atmospheric parameters are summarised in Tab. 11, showing
the proportion of the posterior probability that lies in the different mass
orderings and $\theta_{23}$ octant, with and without the reactor constraint. A
flat prior distribution on both $\Delta{}m^{2}_{32}$ and $\sin^{2}\theta_{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 $\theta_{23}$ octant; 4.21 (1.83) for the normal over inverted mass
ordering; and a combined factor of 1.58 (0.63) for upper $\theta_{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 10.2307/2291091 , 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\sigma$”) 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.

Figure 22: The marginalised posterior probability density of
$\sin^{2}\theta_{23}$ with (red) and without (blue) the reactor constraint on
$\sin^{2}\theta_{13}$ applied. The shaded areas show the 68% and 95% regions
of highest posterior density, equivalent to the 1$\sigma$ and 2$\sigma$
credible intervals.

####  8.1.2 The CP-violating phase $\delta_{\scriptscriptstyle\mathrm{CP}}$,
and $\sin^{2}\theta_{13}$

A comparison of $\sin^{2}\theta_{13}-\delta_{\scriptscriptstyle\mathrm{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$\sigma$ regions
overlapping, comparable with the reactor constraint. A comparison of the
$\delta_{\scriptscriptstyle\mathrm{CP}}$ posterior distributions is shown in
Fig. 24, showing the impact of the reactor constraint on T2K’s
$\delta_{\scriptscriptstyle\mathrm{CP}}$ result. The external constraint
breaks the partially degenerate effects of $\sin^{2}\theta_{13}$ and
$\delta_{\scriptscriptstyle\mathrm{CP}}$ on the $\nu_{e}$ appearance, leading
to the $\nu$-mode 1R$e$ and 1R$e$1d$e$ selections having a larger sensitivity
to $\delta_{\scriptscriptstyle\mathrm{CP}}$.

Figure 23: 68% and 90% credible intervals from the marginalised
$\sin^{2}\theta_{13}-\delta_{\scriptscriptstyle\mathrm{CP}}$ posterior
distribution with (red) and without (blue) the reactor constraint (green band)
applied, marginalised over both mass orderings. Figure 24: The marginalised
posterior probability density of $\delta_{\scriptscriptstyle\mathrm{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$\sigma$ and 2$\sigma$ credible intervals.

####  8.1.3 The Jarlskog invariant

The sampled posterior probability density is in part a function of the PMNS
mixing angles and $\delta_{\scriptscriptstyle\mathrm{CP}}$, which means the
probability distribution for the Jarlskog invariant PhysRevLett.55.1039 ;
JARLSKOG2005323 ,

$J=\sin\theta_{13}\cos^{2}\theta_{13}\sin\theta_{12}\cos\theta_{12}\sin\theta_{23}\cos\theta_{23}\sin\delta_{\scriptscriptstyle\mathrm{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
$\theta_{12}$ and $\theta_{13}$ constraints. The preference for
$\sin^{2}\theta_{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 $\delta_{\scriptscriptstyle\mathrm{CP}}$ constraint, a
preference for negative values of $J$ emerges. The blue curve in Fig. 25 is
recreated in Fig. 26 showing the $1\sigma$, $2\sigma$, and $3\sigma$ credible
intervals. Two-dimensional credible regions for the Jarlskog invariant against
both $\sin^{2}\theta_{23}$ and $\delta_{\scriptscriptstyle\mathrm{CP}}$ are
included in App. C.

Figure 25: Posterior probability distributions for the Jarlskog invariant
using a prior distribution from the 2019 PDG reactor constraint on
$\theta_{13}$Tanabashi:2018oca  (red), prior from all parameters except
sampling $\theta_{23}$ from the T2K posterior (green), and the full T2K
posterior (blue). All three posterior probabilities used a prior probability
distribution flat in $\delta_{\scriptscriptstyle\mathrm{CP}}$.

Although this analysis does not rule out CP-conserving values of
$\delta_{\scriptscriptstyle\mathrm{CP}}$ at $2\sigma$, it does rule out $J=0$
at the $2\sigma$ level and excludes the $J>0.17$ region at $>3\sigma$ with a
flat prior probability in $\delta_{\scriptscriptstyle\mathrm{CP}}$. The
dependence of $J$ on the choice of a prior flat in
$\delta_{\scriptscriptstyle\mathrm{CP}}$ or flat in
$\sin\delta_{\scriptscriptstyle\mathrm{CP}}$ is shown in Fig. 26. The prior
flat in $\sin\delta_{\scriptscriptstyle\mathrm{CP}}$ flattens out the Jarlskog
distribution, which in turn slightly expands the 2$\sigma$ credible interval
to where $J=0$ is just included. These conclusions agree with previous studies
on the impact of the $\delta_{\scriptscriptstyle\mathrm{CP}}$ prior at T2K
Abe:2021gky ; T2K:2019bcf .

Figure 26: Posterior probability distributions for the Jarlskog invariant taken from posterior distributions with priors that are either flat in $\delta_{\scriptscriptstyle\mathrm{CP}}$ (blue) or flat in $\sin\delta_{\scriptscriptstyle\mathrm{CP}}$ (orange). 1$\sigma$, 2$\sigma$, and 3$\sigma$ credible intervals are shown as the region between the vertical black solid line and the specified vertical dashed lines. Selection |  $p$-value
---|---
Shape | Rate
1R$\mu$ |  $\nu$-mode | 0.48 | 0.18
$\overline{\nu}$-mode | 0.85 | 0.74
1R$e$ |  $\nu$-mode | 0.19 | 0.49
$\overline{\nu}$-mode | 0.61 | 0.39
1R$e$1d$e$ |  $\nu$-mode | 0.86 | 0.22
All | 0.73 | 0.30

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.

####  8.1.4 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 Gelman_Post  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 Tab. 12. Good $p$-values
are demonstrated for all cases.

(a) $\nu$-mode 1R$\mu$

(b) $\nu$-mode 1R$e$

(c) $\nu$-mode 1R$e$1d$e$

(d) $\overline{\nu}$-mode 1R$\mu$

(e) $\overline{\nu}$-mode 1R$e$

Figure 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}\theta_{13}$ is constrained
from T2K data alone with no reactor constraint applied.

###  8.2 Frequentist results

As in previous T2K analyses, the frequentist results are obtained using the
marginal likelihood $\mathcal{L}_{\mathrm{marg}}(\theta)$ $=$
$\int\mathrm{d}\eta\,p(\eta)$ $\mathcal{L}(\theta,\eta)$ as the test
statistic. Here, $\mathcal{L}(\theta,\eta)$ is the binned Poisson likelihood
for the parameter of interest, $\theta$, and the nuisance parameters, $\eta$.
The statistical treatment of nuisance parameters in the fit is thus identical
to the Bayesian analysis and assumes a prior probability distribution
$p(\eta)$. The numerical integration is performed by varying systematic
uncertainties with a Gaussian covariance matrix from the ND analysis in Sec. 6
as a constraint, and varying the other oscillation parameters with a flat
prior probability distribution on $\sin^{2}\theta_{23}$,
$\delta_{\scriptscriptstyle\mathrm{CP}}$, $\Delta{}m^{2}_{32}$, and
$\sin^{2}2\theta_{13}$, or a Gaussian prior on $\sin^{2}2\theta_{13}$.
Confidence intervals and regions are constructed with two different methods.
For critical parameters with known boundary effects, the Feldman–Cousins (FC)
method Feldman:1997qc  is utilised to calculate the coverage. This is
performed for the result using the reactor constraint, on the one-dimensional
confidence intervals in $\delta_{\scriptscriptstyle\mathrm{CP}}$ and
$\sin^{2}\theta_{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
$\Delta\chi^{2}(\theta)=\chi^{2}(\theta)-\min_{\theta^{\prime}}\chi^{2}(\theta^{\prime})$
values via Wilks’s theorem 10.1214/aoms/1177732360 , where
$\chi^{2}=-2\ln\mathcal{L}_{\mathrm{marg}}$. Whether $\Delta\chi^{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 analysis. The frequentist analysis bins the $e$-like FD samples in
reconstructed lepton angle and reconstructed lepton momentum, and the
$\mu$-like samples in reconstructed lepton angle and the reconstructed
neutrino energy, defined in the same way as in the Bayesian analysis presented
in Sec. 8.1. In previous analyses, the $\mu$-like samples were binned only in
reconstructed neutrino energy, and adding the lepton angle information
increases the $1\sigma$ expected sensitivity to $\Delta{}m^{2}_{32}$ by
$\mathcal{O}(1\%)$.

Parameter | With reactor constraint | Without reactor constraint
---|---|---
Normal ordering | Inverted ordering | Normal ordering | Inverted ordering
$\delta_{\scriptscriptstyle\mathrm{CP}}$ (rad.) | $-1.97_{-0.62}^{+0.97}$ | $-1.44_{-0.59}^{+0.56}$ | $-2.22_{-0.81}^{+1.25}$ | $-1.29_{-0.83}^{+0.72}$
$\sin^{2}\theta_{13}/10^{-3}$ | — | — | $28.0_{-6.5}^{+2.8}$ | $31.0_{-6.9}^{+3.0}$
$\sin^{2}\theta_{23}$ | $0.561_{-0.038}^{+0.019}$ | $0.563_{-0.032}^{+0.017}$ | $0.467_{-0.018}^{+0.106}$ | $0.466_{-0.019}^{+0.103}$
$\Delta{}m^{2}_{32}/10^{-3}\,(\mathrm{eV^{2}})$ | $2.494_{-0.058}^{+0.041}$ | — | $2.495_{-0.058}^{+0.041}$ | —
$|\Delta m^{2}_{31}|/10^{-3}\,(\mathrm{eV^{2}})$ | — | $2.463_{-0.056}^{+0.042}$ | — | $2.463_{-0.055}^{+0.043}$
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 $\Delta\chi^{2}$ method.

Global best-fit values are given in Tab. 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
$\delta_{\scriptscriptstyle\mathrm{CP}}$ and $\sin^{2}\theta_{23}$. All the
following results are from the fit to data using the reactor constraint.

####  8.2.1 The CP-violating phase $\delta_{\scriptscriptstyle\mathrm{CP}}$,
and mass ordering

Fig. 28 shows the $\Delta\chi^{2}$ distributions for
$\delta_{\scriptscriptstyle\mathrm{CP}}$ in both mass orderings with FC-
adjusted confidence intervals, which are also summarised in Tab. 14. A large
region of $\sin\delta_{\scriptscriptstyle\mathrm{CP}}>0$ is excluded at
$>3\sigma$ confidence level (CL), whereas the CP-conserving values
$\delta_{\scriptscriptstyle\mathrm{CP}}=0,\pi$ are excluded at 90% CL. In
particular, $\delta_{\scriptscriptstyle\mathrm{CP}}=\pi$ is just inside the
$2\sigma$ interval.

Confidence |  $\delta_{\scriptscriptstyle\mathrm{CP}}$ (rad.) | $\sin^{2}\theta_{23}$
---|---|---
Level | Normal ordering | Inverted ordering | Normal ordering | Inverted ordering
$1\sigma$ | $[-2.67,-1.00]$ | — | $[0.529,0.582]$ | —
$90\%$ | $[-3.01,-0.52]$ | $[-1.74,-1.07]$ | $[0.444,0.593]$ | $[0.536,0.584]$
$2\sigma$ | $[-\pi,-0.28]\cup{}[3.10,\pi]$ | $[-2.16,-0.74]$ | $[0.436,0.598]$ | $[0.512,0.592]$
$3\sigma$ | $[-\pi,0.33]\cup{}[2.59,\pi]$ | $[-2.83,-0.14]$ | N/A | N/A
Table 14: FC-corrected confidence intervals for
$\delta_{\scriptscriptstyle\mathrm{CP}}$ and $\sin^{2}\theta_{23}$ from the
fit to data in the frequentist analysis, using the reactor constraint on
$\sin^{2}2\theta_{13}$. The $3\sigma$ FC correction was not computed for
$\sin^{2}\theta_{23}$. Figure 28: The $\Delta\chi^{2}$ distribution in
$\delta_{\scriptscriptstyle\mathrm{CP}}$ from fitting to the data with the
reactor constraint applied. The confidence intervals in the shaded regions are
calculated using the FC method. Figure 29: The $\Delta\chi^{2}$ distribution
in $\delta_{\scriptscriptstyle\mathrm{CP}}$ for incremental modifications from
the previous analysis Abe:2021gky  to this result, for normal and inverted
mass orderings. “E” corresponds to this analysis, except that unlike the main
frequentist result, the $\mu$-like samples do not use the scattering angle
information for better compatibility with the previous T2K analysis.

As was also seen in the Bayesian analysis, the constraint on
$\delta_{\scriptscriptstyle\mathrm{CP}}$ is weaker compared to T2K’s previous
analysis Abe:2021gky ; T2K:2019bcf . Fig. 29 shows the impact on the
$\Delta\chi^{2}$ of $\delta_{\scriptscriptstyle\mathrm{CP}}$ after each update
introduced in this analysis, all of which weaken the
$\delta_{\scriptscriptstyle\mathrm{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
$\delta_{\scriptscriptstyle\mathrm{CP}}$ parameter 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\sigma$
for all $\delta_{\scriptscriptstyle\mathrm{CP}}$ values, mostly consistent
with the expected sensitivity at
$\sin\delta_{\scriptscriptstyle\mathrm{CP}}=-1$. Replacing each sample’s event
distribution in data by the expectation of the model shows that the stronger
$\delta_{\scriptscriptstyle\mathrm{CP}}$ constraint observed in data compared
to the expectation comes from the $\nu$-mode 1R$e$1d$e$ sample.

Figure 30: The $\Delta\chi^{2}$ distribution in
$\delta_{\scriptscriptstyle\mathrm{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 from an
ensemble of toy simulated experiments created with true normal ordering and
$\delta_{\scriptscriptstyle\mathrm{CP}}=-\pi/2$, showing the $\Delta\chi^{2}$
for 68% and 95% of the toys, and their median.

####  8.2.2 Atmospheric oscillation parameters

Figure 31: Confidence regions in $\sin^{2}\theta_{23}-\Delta{}m^{2}_{32}$
($|\Delta m^{2}_{31}|$ in inverted mass ordering) for the data fit with the
reactor constraint applied, obtained with the constant $\Delta\chi^{2}$
method, where in each mass ordering hypothesis a fixed mass ordering is
assumed.

The $\sin^{2}\theta_{23}-\Delta{}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. Fig. 32 shows the contour from a fit using
only the 1R$\mu$ samples, which shows that the constraint is dominated by the
1R$\mu$ samples, with the 1R$e$ samples providing the sensitivity to the
octant of $\theta_{23}$.

Figure 32: Comparison of confidence regions in
$\sin^{2}\theta_{23}-\Delta{}m^{2}_{32}$ for normal ordering, between a full
fit and a fit using only the $\mu$-like samples. The intervals are calculated
with the constant $\Delta\chi^{2}$ method, and applying the reactor
constraint.

The evolution of the $\sin^{2}\theta_{23}-\Delta{}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 $\Delta{}m^{2}_{32}$ constraint comes from
changing the cross-section model and updating the ND constraint. For
$\Delta{}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 Sec. 9. For $\sin^{2}\theta_{23}$ there is a
slight shift in shape from the latest data in the FD.

Figure 33: The $\Delta\chi^{2}$ distribution in $\Delta{}m^{2}_{32}$ and
$\sin^{2}\theta_{23}$ for normal ordering, showing incremental modifications
of the previous analysis Abe:2021gky ; T2K:2019bcf  to this result. “E”
corresponds to this analysis, except that unlike the main frequentist result,
the $\mu$-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 $\Delta\chi^{2}$ distribution for $\sin^{2}\theta_{23}$ is shown in Fig.
34 with the confidence intervals summarised in Tab. 14. The new data at the FD
has reduced compatibility with maximal mixing, which is now outside the FC-
corrected $1\sigma$ confidence interval. Whilst the upper octant is favoured
at $1\sigma$ CL, the data is still compatible with both octant hypotheses at
90% CL. These results are compatible with the sensitivity, shown in Fig. 35.

Figure 34: The $\Delta\chi^{2}$ distribution in $\sin^{2}\theta_{23}$ for
fitting to the data with the reactor constraint applied. The confidence
intervals in the shaded regions are calculated using the FC method. Figure 35:
The $\Delta\chi^{2}$ distribution in $\sin^{2}\theta_{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}\theta_{23}=0.56$.

###  8.3 Cross-fitter comparisons

To compare the consistency of the Bayesian analysis described in Sec. 8.1 and
the frequentist analysis described in Sec. 8.2, the Bayesian posterior
distributions were recast into frequentist $\Delta\chi^{2}$ distributions
comparable to the frequentist analysis. Fig. 36 shows comparisons of
$\sin^{2}\theta_{23}-\Delta{}m^{2}_{32}$ and
$\sin^{2}\theta_{13}-\delta_{\scriptscriptstyle\mathrm{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_{\mu}-\cos\theta_{\mu}$ 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 $\mu$-like samples. Additionally, the frequentist analysis also uses
the reconstructed lepton scattering angle for the $\mu$-like samples to
disentangle systematic uncertainties related to energy scale and mis-
reconstructed backgrounds at the FD. 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.

Figure 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 Sec. 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}\theta_{13}$.

###  8.4 Comparisons with other experiments

In the global context of neutrino oscillation experiments, these results
provide leading constraints on both the atmospheric oscillation parameters,
$\Delta{}m^{2}_{32}$ and $\sin^{2}\theta_{23}$, and the CP-violating phase,
$\delta_{\scriptscriptstyle\mathrm{CP}}$. Whereas the other experiments
profile over parameters to calculate the $\Delta\chi^{2}$, T2K instead
calculates the marginal likelihood for the $\Delta\chi^{2}$. Fig. 37 shows the
90% confidence regions in $\sin^{2}\theta_{23}-\Delta{}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. Fig. 38 compares the
90% confidence regions in
$\sin^{2}\theta_{23}-\delta_{\scriptscriptstyle\mathrm{CP}}$ for both
orderings to NOvA and SK. The confidence intervals on $\sin^{2}\theta_{23}$
significantly overlap, as do the intervals for
$\delta_{\scriptscriptstyle\mathrm{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
$\pi<\delta_{\scriptscriptstyle\mathrm{CP}}<2\pi$ region, with a weak
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.

Figure 37: Comparison of the 90% confidence regions in
$\sin^{2}\theta_{23}-\Delta{}m^{2}_{32}$ for normal ordering with NOvA
NOvA:2021nfi , Super-K Super-Kamiokande:2017yvm , IceCube
icecube_2018_10_21234_B4105H , and MINOS+ MINOS:2020llm . The NOvA and IceCube
constraints are 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 $\Delta\chi^{2}$ method, assuming
normal ordering. Figure 38: Comparison of 90% confidence regions in
$\sin^{2}\theta_{23}-\delta_{\scriptscriptstyle\mathrm{CP}}$ over both mass
orderings with NOvA NOvA:2021nfi  and Super-K Super-Kamiokande:2017yvm . The
T2K and NOvA confidence regions have been computed using the FC method,
whereas the Super-K results are obtained with the constant $\Delta\chi^{2}$
method.

##  9 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 Sec. 5.3, on the oscillation parameter constraints. The oscillation
analysis in Sec. 8 had these uncertainty inflation strategies applied, and
this section summarises the procedure, with details provided in App. B.

###  9.1 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
$\delta_{\scriptscriptstyle\mathrm{CP}}=0$ and non-maximal
$\sin^{2}\theta_{23}$, detailed in App. B, Tab. 17. This section only presents
results with the Asimov data near T2K’s best-fit point. The PDG 2019 reactor
constraint on $\sin^{2}\theta_{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 robust, 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\sigma$ confidence interval of an
oscillation parameter, relative to the $1\sigma$ from systematic uncertainties
in the expected sensitivity analysis. An example is discussed in Sec. B.1.

  * •

ND to FD extrapolation issues: Some alternate models may not produce a
significant bias on 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 predicted
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 Hyper-Kamiokande:2018ofw  and DUNE
DUNE:2020ypp . An example is discussed in Sec. B.2.

All the individual biases are summed in quadrature and are used to inflate the
confidence interval for $\Delta{}m^{2}_{32}$, due to its simple Gaussian
probability density. For $\delta_{\scriptscriptstyle\mathrm{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
$\delta_{\scriptscriptstyle\mathrm{CP}}$. This is done by calculating the
difference in the $\Delta\chi^{2}$ distribution for
$\delta_{\scriptscriptstyle\mathrm{CP}}$ for the Asimov data and the simulated
data, and adding the difference to the $\Delta\chi^{2}$ distribution for
$\delta_{\scriptscriptstyle\mathrm{CP}}$ from the data, where the 90%
confidence interval was calculated using the FC method mentioned in Sec. 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
mis-modelling 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 $\hskip 1.42262pt\hbox
to0.0pt{\hss\raisebox{4.30554pt}{\scalebox{0.5}[0.5]{(}}}\overline{\nu}\hbox
to0.0pt{\raisebox{4.30554pt}{\scalebox{0.5}[0.5]{)}}\hss}_{e}$ at the FD near
the oscillation maximum will likely lead to a stronger constraint on
$\delta_{\scriptscriptstyle\mathrm{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 mis-modelling 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 Sec. 5.3. Details of the simulated data study procedure and two
examples are provided in App. B.

###  9.2 Results

Tab. 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 App. B, Tab. 18. The impact of the simulated data studies
on $\sin^{2}\theta_{23}$ and $\delta_{\scriptscriptstyle\mathrm{CP}}$ was
found to be small compared to the impact of statistical and systematic
uncertainties. The largest bias on $\Delta{}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 quadrature444The non-CCQE, data-driven low pion
momentum, low $Q^{2}$ pion suppression from MINERvA, pion SI, the CCQE form
factor with the largest impact (z-expansion, upper variation), Martini 2p2h,
and the removal energy simulated data studies were selected. to avoid double
counting similar physics effects, leading to an overall smearing on
$\Delta{}m^{2}_{32}$ of $1.35\times 10^{-5}\mathrm{eV^{2}}$. For comparison,
the overall uncertainty on $\Delta{}m^{2}_{32}$ from the expected sensitivity
study, before the simulated data procedure, was $5.7\times
10^{-5}\mathrm{eV^{2}}$ and is dominated by the uncertainty from statistics,
which was $5.3\times 10^{-5}\mathrm{eV^{2}}$. Generally, the simulated data
studies had a smaller impact on $\delta_{\scriptscriptstyle\mathrm{CP}}$, due
to its uncertainty being dominated by the statistics in the electron-like
selections at the FD.

Simulated data set | Relative to | $\sin^{2}\theta_{23}$ | $\Delta{}m^{2}_{32}$ | $\delta_{\scriptscriptstyle\mathrm{CP}}$
---|---|---|---|---
CCQE z-exp high | Total | 0.3% | 2.1% | 0.4%
Syst. | 0.7% | 5.7% | 1.7%
CCQE removal energy | Total | 0.0% | 4.8% | 1.3%
Syst. | 0.0% | 13.4% | 5.2%
Non-CCQE | Total | 8.7% | 11.8% | 1.7%
Syst. | 21.3% | 32.7% | 6.9%
2p2h Martini | Total | 0.7% | 2.7% | 0.4%
Syst. | 1.6% | 7.3% | 1.6%
MINERvA pion tune | Total | 2.9% | 2.5% | 0.9%
Syst. | 7.2% | 6.8% | 3.5%
Data-driven pion | Total | 4.7% | 6.5% | 1.0%
Syst. | 11.6% | 17.9% | 3.9%
Pion SI | Total | 0.7% | 20.8 % | 1.0%
Syst. | 1.9% | 57.8 % | 4.6%
Table 15: Biases on the main oscillation parameters for each simulated data
set, calculated as the shift in the middle of the $1\sigma$ confidence
interval relative to the overall uncertainty from systematic sources (“Syst.”)
and the total (“Total”) to one decimal place.

In the previous T2K analyses Abe:2021gky ; T2K:2019bcf , the simulated data
study for the nucleon removal energy had a significant impact, especially on
$\Delta{}m^{2}_{32}$, and an additional uncertainty was introduced. In this
analysis, the updated nucleon removal energy uncertainty, described in Sec. 5,
has caused it to no longer be the dominant source of systematic uncertainty.

Tab. 16 shows the changes to the 90% confidence interval for
$\delta_{\scriptscriptstyle\mathrm{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 CP-conserving values presented in Sec. 8.

Simulated data set | Change to 90% C.L. of $\delta_{\scriptscriptstyle\mathrm{CP}}$
---|---
-3.01 | -0.52
CCQE z-exp high | 0.05 | 0.04
CCQE removal energy | 0.00 | 0.02
Non-CCQE | 0.06 | 0.09
2p2h Martini | 0.04 | 0.04
MINERvA pion tune | 0.05 | 0.04
Data-driven pion | 0.07 | 0.04
Pion SI | 0.00 | 0.01
Table 16: Shifts of the 90% confidence interval boundaries of
$\delta_{\scriptscriptstyle\mathrm{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)
$\delta_{\scriptscriptstyle\mathrm{CP}}$ interval edge from the data fit. Only
the absolute size of the shift is taken into account.

##  10 Conclusions

The T2K collaboration has measured the three-flavour PMNS neutrino oscillation
parameters $\Delta{}m^{2}_{32}$, $\sin^{2}\theta_{13}$, $\sin^{2}\theta_{23}$,
$\delta_{\scriptscriptstyle\mathrm{CP}}$, the Jarlskog invariant $J$, and the
mass ordering, using the statistics at the FD equivalent to $3.6\times
10^{21}$ POT. T2K continues to favour neutrino oscillations with near-maximal
CP violation, in the upper octant of $\sin^{2}\theta_{23}$, in the normal mass
ordering, with a $\sin^{2}\theta_{13}$ consistent with the measurements by
reactor experiments.

The analysis included $4.72\times 10^{20}$ POT more neutrino data at the FD,
and $5.73(4.48)\times 10^{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
$\Delta{}m^{2}_{32}$ and $\delta_{\scriptscriptstyle\mathrm{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 T2K:2019eao ; Oyama:2020kev  the neutrino beam for the Hyper-
Kamiokande experiment Hyper-Kamiokande:2018ofw . 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 T2K:2019eao ; Oyama:2020kev . 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 ultra-pure water Super-Kamiokande:2021the , drastically
increasing the efficiency in tagging 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 T2K:2019bbb ; Attie:2021yeh ;
Blondel:2020hml  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.

###### Acknowledgements.

We thank the J-PARC staff for superb accelerator performance. We thank the
CERN NA61/SHINE Collaboration for providing valuable particle production data.
We acknowledge the support of MEXT, JSPS KAKENHI and bilateral programs,
Japan; NSERC, the NRC, and CFI, Canada; the CEA and CNRS/IN2P3, France; the
DFG, Germany; the INFN, Italy; the Ministry of Science and Education and the
National Science Centre, Poland; the RSF, RFBR and the Ministry of Science and
Higher Education, Russia; MICINN and ERDF funds and CERCA program, Spain; the
SNSF and SERI, 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,
NII for SINET5, the BC DRI Group, Prairie DRI Group, ACENET, SciNet, and
CalculQuebec consortia in the Digital Research Alliance of Canada, and GridPP
in the United Kingdom. In addition, the participation of individual
researchers and institutions has been further supported by funds from the ERC
(FP7), “la Caixa” Foundation, the European Union’s Horizon 2020 Research and
Innovation Programme under the Marie Sklodowska-Curie grant; the JSPS, Japan;
the Royal Society, UK; French ANR; 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.

## References

  * (1) K. Abe et al., (The T2K collaboration), Phys. Rev. D 103 no. 11, (2021) 112008.
  * (2) K. Abe et al., (The T2K collaboration), Nature 580 no. 7803, (2020) 339–344. [Erratum: Nature 583, E16 (2020)].
  * (3) Y. Fukuda et al., (The Super-Kamiokande collaboration), Phys. Rev. Lett. 81 (1998) 1562–1567.
  * (4) Q. R. Ahmad et al., (The SNO collaboration), Phys. Rev. Lett. 89 (2002) 011301.
  * (5) K. Eguchi et al., (The KamLAND collaboration), Phys. Rev. Lett. 90 (2003) 021802.
  * (6) M. H. Ahn et al., (The K2K collaboration), Phys. Rev. D74 (2006) 072003.
  * (7) P. Adamson et al., (The MINOS collaboration), Phys.Rev.Lett. 112 (2014) 191801.
  * (8) B. Pontecorvo, Sov. Phys. JETP 26 (1968) 984–988.
  * (9) Z. Maki, M. Nakagawa, and S. Sakata, Prog. Theor. Phys. 28 (1962) 870–880.
  * (10) S. Fukuda et al., (The Super-Kamiokande collaboration), Phys. Lett. B 539 (2002) 179–187.
  * (11) M. Tanabashi et al., (The Particle Data Group collaboration), Phys. Rev. D 98 no. 3, (2018) 030001. [Updated 2019 tables used for $\sin^{2}\theta_{13}$ constraint].
  * (12) K. Abe et al., (The T2K collaboration), Phys. Rev. Lett. 107 (2011) 041801.
  * (13) F. P. An et al., (The Daya Bay collaboration), Phys. Rev. Lett. 108 (2012) 171803.
  * (14) J. K. Ahn et al., (The RENO collaboration), Phys. Rev. Lett. 108 (2012) 191802.
  * (15) Y. Abe et al., (The Double Chooz collaboration), Phys. Rev. Lett. 108 (2012) 131801.
  * (16) K. Abe et al., (The T2K collaboration), Phys.Rev.Lett. 112 (2014) 061802.
  * (17) P. Adamson et al., (The NOvA collaboration), Phys. Rev. Lett. 116 no. 15, (2016) 151806.
  * (18) M. A. Acero et al., (The NOvA collaboration), Phys. Rev. Lett. 123 no. 15, (2019) 151803.
  * (19) L. Wolfenstein, Phys. Rev. D 17 (1978) 2369–2374.
  * (20) S. P. Mikheev and A. Y. Smirnov, Nuovo Cim. C 9 (1986) 17–26.
  * (21) J. Arafune, M. Koike, and J. Sato, Phys. Rev. D 56 (1997) 3093–3099. [Erratum: Phys.Rev.D 60, 119905 (1999)].
  * (22) C. Jarlskog, Phys. Rev. Lett. 55 (Sep, 1985) 1039–1042.
  * (23) C. Jarlskog, Physics Letters B 609 no. 3, (2005) 323–329.
  * (24) K. Abe et al., (The T2K collaboration), Nucl. Instrum. Meth. A 659 (2011) 106–135.
  * (25) S. Bhadra et al., Nucl. Instrum. Meth. A 703 (2013) 45–58.
  * (26) K. Suzuki, Nucl. Phys. B, Proc. Suppl. 229-232 (2012) 453. Neutrino 2010.
  * (27) K. Abe et al., (The T2K, J-PARC Neutrino Facility Group collaboration), arXiv:1908.05141 [physics.ins-det].
  * (28) Y. Oyama, PoS NuFact2019 (2020) 054.
  * (29) K. Abe et al., Nucl. Instrum. Meth. A 694 (2012) 211 – 223.
  * (30) A. Astbury et al., CERN-SPSC-78-06 (1, 1978) .
  * (31) M. Barranco-Luque et al., (The UA1 collaboration), Nucl. Instrum. Meth. 176 (1980) 175.
  * (32) J. Altegoer et al., (The NOMAD collaboration), Nucl. Instrum. Meth. A 404 (1998) 96–128.
  * (33) S. Aoki et al., Nucl. Instrum. Meth. A 698 (2013) 135 – 146.
  * (34) D. Allan et al., Journal of Instrumentation 8 no. 10, (10, 2013) P10019–P10019.
  * (35) S. Assylbekov et al., Nucl. Instrum. Meth. A 686 (2012) 48 – 63.
  * (36) P.-A. Amaudruz et al., Nucl. Instrum. Meth. A 696 (2012) 1 – 31.
  * (37) N. Abgrall et al., Nucl. Instrum. Meth. A 637 no. 1, (2011) 25 – 46.
  * (38) Y. Fukuda et al., (The Super-Kamiokande collaboration), Nucl. Instrum. Meth. A 501 (2003) 418–462.
  * (39) K. Abe et al., (The Super-Kamiokande collaboration), Nucl. Instrum. Meth. A 1027 (2022) 166248.
  * (40) N. Abgrall et al., (The NA61/SHINE collaboration), Eur. Phys. J. C 76 no. 11, (2016) 617.
  * (41) N. Abgrall et al., (The NA61/SHINE collaboration), Eur. Phys. J. C 76 no. 2, (2016) 84.
  * (42) O. Benhar, A. Fabrocini, S. Fantoni, and I. Sick, Nucl. Phys. A579 (1994) 493–517.
  * (43) E. S. Pinzon Guerra et al., Phys. Rev. D 99 no. 5, (2019) 052007.
  * (44) K. Abe et al., (The T2K collaboration), Phys. Rev. D87 no. 1, (2013) 012001. [Addendum: Phys. Rev.D87,no.1,019902(2013)].
  * (45) T. T. Böhlen, F. Cerutti, M. P. W. Chin, A. Fassò, A. Ferrari, P. G. Ortega, A. Mairani, P. R. Sala, G. Smirnov, and V. Vlachoudis, Nucl. Data Sheets 120 (2014) 211–214.
  * (46) A. Ferrari, P. R. Sala, A. Fasso, and J. Ranft, “FLUKA: A multi-particle transport code (Program version 2005),” 10, 2005.  doi: 10.2172/877507.
  * (47) K. Abe et al., (The T2K collaboration), Phys. Rev. D 87 no. 1, (2013) 012001. [Addendum: Phys.Rev.D 87, 019902 (2013)].
  * (48) S. Agostinelli et al., (The GEANT4 collaboration), Nucl. Instrum. Meth. A 506 (2003) 250–303.
  * (49) M. Friend et al., “Neutrino beam flux prediction 2020,” Nov., 2021.  doi: 10.5281/zenodo.5734307.
  * (50) N. Abgrall et al., (The NA61/SHINE collaboration), Eur. Phys. J. C79 no. 2, (2019) 100.
  * (51) L. Alvarez-Ruso et al., (The NuSTEC collaboration), Prog. Part. Nucl. Phys. 100 (2018) 1–68.
  * (52) U. Mosel, Ann. Rev. Nucl. Part. Sci. 66 (2016) 171–195.
  * (53) T. Katori and M. Martini, J. Phys. G 45 no. 1, (2018) 013001.
  * (54) Y. Hayato and L. Pickering, The European Physical Journal Special Topics (6, 2021) .
  * (55) Y. Hayato, Acta Phys. Polon. B 40 (2009) 2477–2489.
  * (56) Y. Hayato, Nucl. Phys. B, Proc. Suppl. 112 no. 1, (2002) 171–176.
  * (57) A. Bodek and T. Cai, Eur. Phys. J. C79 no. 4, (2019) 293.
  * (58) R. Bradford et al., Nucl. Phys. B, Proc. Suppl. 159 (2006) 127 – 132.
  * (59) J. Nieves, I. Ruiz Simo, and M. J. Vicente Vacas, Phys. Lett. B707 (2012) 72–75.
  * (60) D. Rein and L. M. Sehgal, Annals Phys. 133 (1981) 79–153.
  * (61) K. M. Graczyk, J. Żmuda, and J. T. Sobczyk, Phys. Rev. D 90 no. 9, (2014) 093001.
  * (62) K. M. Graczyk and J. T. Sobczyk, Phys. Rev. D 77 (2008) 053001. [Erratum: Phys.Rev.D 79, 079903 (2009)].
  * (63) C. Berger and L. M. Sehgal, Phys. Rev. D 76 (2007) 113004.
  * (64) K. M. Graczyk and J. T. Sobczyk, Phys. Rev. D 77 (2008) 053003.
  * (65) K. S. Kuzmin, V. V. Lyubushkin, and V. A. Naumov, Phys. Part. Nucl. 35 (2004) S133–S138.
  * (66) C. Berger and L. M. Sehgal, Phys. Rev. D 79 (2009) 053003.
  * (67) D. Rein and L. M. Sehgal, Nucl. Phys. B 223 (1983) 29–44.
  * (68) D. Rein, Nucl. Phys. B 278 (1986) 61–77.
  * (69) M. Glück, E. Reya, and A. Vogt, Eur. Phys. J. C 5 no. 3, (1998) 461–470.
  * (70) A. Bodek and U. K. Yang, AIP Conf. Proc. 670 no. 1, (2003) 110–117.
  * (71) A. Bodek and U. K. Yang, AIP Conf. Proc. 792 no. 1, (2005) 257–260.
  * (72) T. Sjostrand, Computer Physics Communications 82 no. 1, (1994) 74 – 89.
  * (73) L. Aliaga et al., (The NuSTEC collaboration), arXiv:2011.07166 [hep-ph].
  * (74) L. L. Salcedo, E. Oset, M. J. Vicente-Vacas, and C. Garcia-Recio, Nucl. Phys. A 484 (1988) 557–592.
  * (75) E. Oset and L. L. Salcedo, Nucl. Phys. A 468 (1987) 631–652.
  * (76) P. Gueye et al., Phys. Rev. C60 (1999) 044308.
  * (77) D. Dutta et al., (The JLab E91013 collaboration), Phys. Rev. C68 (2003) 064603.
  * (78) M. Leuschner et al., Phys. Rev. C49 (1994) 955–967.
  * (79) P. K. A. de Witt Huberts, Journal of Physics G: Nuclear and Particle Physics 16 no. 4, (Apr, 1990) 507–544.
  * (80) D. Ruterbories et al., (The MINERvA collaboration), Phys. Rev. D 99 no. 1, (2019) 012004.
  * (81) P. A. Rodrigues et al., (The MINERvA collaboration), Phys. Rev. Lett. 116 (2016) 071802. [Addendum: Phys.Rev.Lett. 121, 209902 (2018)].
  * (82) K. Abe et al., (The T2K collaboration), Phys. Rev. D 101 no. 11, (2020) 112001.
  * (83) K. Abe et al., (The T2K collaboration), Phys. Rev. D 101 no. 11, (2020) 112004.
  * (84) P. Stowell et al., JINST 12 no. 01, (2017) P01016.
  * (85) K. L. Miller et al., Phys. Rev. D 26 (1982) 537–542.
  * (86) S. J. Barish et al., Phys. Rev. D 16 (1977) 3103.
  * (87) N. J. Baker et al., Phys. Rev. D 23 (1981) 2499–2505.
  * (88) D. Allasia et al., Nucl. Phys. B 343 (1990) 285–309.
  * (89) T. Kitagaki et al., Phys. Rev. D 28 (1983) 436–442.
  * (90) S. K. Singh, Nucl. Phys. B 36 (1972) 419–435.
  * (91) V. Bernard, L. Elouadrhiri, and U.-G. Meissner, J. Phys. G 28 (2002) R1–R35.
  * (92) B. Bhattacharya, R. J. Hill, and G. Paz, Phys. Rev. D 84 (2011) 073006.
  * (93) K. Abe et al., (The T2K collaboration), Phys. Rev. D 101 no. 11, (2020) 112001.
  * (94) K. Abe et al., (The T2K collaboration), Phys. Rev. D 101 no. 11, (2020) 112004.
  * (95) M. Martini, M. Ericson, G. Chanfray, and J. Marteau, Phys. Rev. C80 (2009) 065501.
  * (96) G. D. Megias, J. E. Amaro, M. B. Barbaro, J. A. Caballero, and T. W. Donnelly, Phys. Rev. D94 (2016) 013012.
  * (97) I. Ruiz Simo, J. E. Amaro, M. B. Barbaro, A. De Pace, J. A. Caballero, and T. W. Donnelly, J. Phys. G44 no. 6, (2017) 065105.
  * (98) K. Abe et al., (The T2K collaboration), Phys. Rev. D 96 no. 9, (2017) 092006. [Erratum: Phys.Rev.D 98, 019902 (2018)].
  * (99) G. M. Radecky et al., Phys. Rev. D 25 (1982) 1161–1173. [Erratum: Phys.Rev.D 26, 3297 (1982)].
  * (100) M. Derrick et al., Phys. Rev. D 23 (1981) 569.
  * (101) T. Kitagaki et al., Phys. Rev. D 34 (1986) 2554–2565.
  * (102) N. J. Baker et al., Phys. Rev. D 23 (1981) 2495–2498.
  * (103) C. Wilkinson, P. Rodrigues, S. Cartwright, L. Thompson, and K. McFarland, Phys. Rev. D 90 no. 11, (2014) 112017.
  * (104) A. A. Aguilar-Arevalo et al., (The MiniBooNE collaboration), Phys. Rev. D 83 (2011) 052007.
  * (105) A. A. Aguilar-Arevalo et al., (The MiniBooNE collaboration), Phys. Rev. D 81 (2010) 013005.
  * (106) A. A. Aguilar-Arevalo et al., (The MiniBooNE collaboration), Phys. Rev. D 83 (2011) 052009.
  * (107) B. Eberly et al., (The MINERvA collaboration), Phys. Rev. D 92 no. 9, (2015) 092008.
  * (108) T. Le et al., (The MINERvA collaboration), Phys. Lett. B 749 (2015) 130–136.
  * (109) O. Altinok et al., (The MINERvA collaboration), Phys. Rev. D 96 no. 7, (2017) 072003.
  * (110) C. L. McGivern et al., (The MINERvA collaboration), Phys. Rev. D 94 no. 5, (2016) 052005.
  * (111) A. Higuera et al., (The MINERvA collaboration), Phys. Rev. Lett. 113 no. 26, (2014) 261802.
  * (112) T. Yang, C. Andreopoulos, H. Gallagher, K. Hoffmann, and P. Kehayias, Eur. Phys. J. C 63 (2009) 1–10.
  * (113) C. Andreopoulos et al., Nucl. Instrum. Meth. A A614 (2010) 87–104.
  * (114) E. S. Pinzon Guerra et al., (The DUET collaboration), Phys. Rev. C95 no. 4, (2017) 045203.
  * (115) M. Day and K. S. McFarland, Phys. Rev. D 86 (2012) 053003.
  * (116) K. S. Kim and M. K. Cheoun, AIP Conf. Proc. 1189 no. 1, (2009) 163–168.
  * (117) J. Engel, Phys. Rev. C 57 (1998) 2004–2009.
  * (118) K. Abe et al., (The T2K collaboration), Phys. Rev. D 98 no. 3, (2018) 032003.
  * (119) X. G. Lu et al., (The MINERvA collaboration), Phys. Rev. Lett. 121 no. 2, (2018) 022504.
  * (120) S. Dolan, arXiv:1810.06043 [hep-ex].
  * (121) D. Ruterbories et al., (The MINERvA collaboration), Phys. Rev. Lett. 129 no. 2, (2022) 021803.
  * (122) C. Adamuscin, E. Tomasi-Gustafsson, E. Santopinto, and R. Bijker, Phys. Rev. C 78 (2008) 035201.
  * (123) R. A. Smith and E. J. Moniz, Nucl. Phys. B 43 (1972) 605. [Erratum: Nucl.Phys.B 101, 547 (1975)].
  * (124) P. Stowell et al., Phys. Rev. D 100 no. 7, (2019) 072005.
  * (125) N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, Journal of Chemical Physics 21 no. 6, (1970) .
  * (126) W. K. Hastings, Biometrika 57 no. 1, (04, 1970) 97–109.
  * (127) F. James and M. Roos, Comput. Phys. Commun. 10 (1975) 343–367.
  * (128) K. Abe et al., (The T2K collaboration), Phys. Rev. D 91 no. 7, (2015) 072010.
  * (129) R. J. Barlow and C. Beeston, Comput. Phys. Commun. 77 (1993) 219–228.
  * (130) J. S. Conway in PHYSTAT 2011, pp. 115–120.  2011\.
  * (131) C. A. Argüelles, A. Schneider, and T. Yuan, JHEP 06 (2019) 030.
  * (132) M. Jiang et al., (The Super-Kamiokande collaboration), PTEP 2019 no. 5, (2019) 053F01.
  * (133) R. E. Kass and A. E. Raftery, Journal of the American Statistical Association 90 no. 430, (1995) 773–795.
  * (134) A. Gelman, X.-L. Meng, and H. Stern, Statistica Sinica 6 (1996) 733–807.
  * (135) G. J. Feldman and R. D. Cousins, Phys. Rev. D 57 (1998) 3873–3889.
  * (136) S. S. Wilks, The Annals of Mathematical Statistics 9 no. 1, (1938) 60 – 62.
  * (137) M. A. Acero et al., (The NOvA collaboration), Phys. Rev. D 106 no. 3, (2022) 032004.
  * (138) K. Abe et al., (The Super-Kamiokande collaboration), Phys. Rev. D 97 no. 7, (2018) 072001.
  * (139) IceCube Collaboration, “Measurement of atmospheric neutrino oscillations with three years of data from the full sky,” 2018\.  doi: 10.21234/B4105H.
  * (140) P. Adamson et al., (The MINOS+ collaboration), Phys. Rev. Lett. 125 no. 13, (2020) 131802.
  * (141) K. Abe et al., (The Hyper-Kamiokande collaboration), arXiv:1805.04163 [physics.ins-det].
  * (142) B. Abi et al., (The DUNE collaboration), arXiv:2002.03005 [hep-ex].
  * (143) K. Abe et al., (The T2K collaboration), arXiv:1901.03750 [physics.ins-det].
  * (144) D. Attié et al., Nucl. Instrum. Meth. A 1025 (2022) 166109.
  * (145) A. Blondel et al., JINST 15 no. 12, (2020) P12003.
  * (146) K. Abe et al., (The T2K collaboration), Phys. Rev. D 101 no. 1, (2020) 012007.
  * (147) M. B. Avanzini et al., Phys. Rev. D 105 no. 9, (2022) 092004.

##  Appendix A Data release

A separate digital data release has been prepared for this publication,
containing all the results from the Bayesian and frequentist analyses
presented in Sec. 8. They are presented in ROOT format, containing TGraph and
TH1D objects. Detailed confidence and credible intervals are provided for the
parameters of interest, with and without the constraint on
$\sin^{2}\theta_{13}$ from reactor experiments applied, in both neutrino mass
orderings. A README is provided outlining the details and conventions of the
data release. The neutrino flux release is separately provided
megan_friend_2021_5734307 .

##  Appendix B Detailed simulated data studies

This section expands on the summary of the simulated data studies provided in
Sec. 9, providing details of two different studies and their effects. All
studies are done with parameters listed in Tab. 17, where one set is close to
the T2K best-fit point, and the other has
$\delta_{\scriptscriptstyle\mathrm{CP}}=0$ and non-maximal
$\sin^{2}\theta_{23}=0.45$ with the parameter set A being presented in this
paper.

First, the non-CCQE study is described in more detail as it produced the
largest bias on oscillation parameters. Then the MINERvA single-pion
suppression study is detailed—having a small impact at T2K statistics, but
will likely have larger impact in higher statistics experiments. The section
finishes with a summary of the impact of all the simulated data studies.

Parameter | Parameter set
---|---
A | B
$\delta_{\scriptscriptstyle\mathrm{CP}}$ | -1.601 | 0.0
$\sin^{2}\theta_{23}$ | 0.528 | 0.45
$\sin^{2}\theta_{13}$ | 0.0218
$\sin^{2}\theta_{12}$ | 0.307
$\Delta{}m^{2}_{32}(\text{eV}^{2})$ | $2.509\times 10^{-3}$
$\Delta{}m^{2}_{21}(\text{eV}^{2})$ | $7.530\times 10^{-5}$
Table 17: The oscillation parameters used in the expected sensitivity fits,
which are used to evaluate the relative size of the biases for the simulated
data studies.

###  B.1 Non-CCQE simulated data set

As described in Sec. 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 Abe:2020jbf ; Abe:2020uub ;
Ruterbories:2018gub ; Rodrigues:2015hik . 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 CC$0\pi$ interactions,
this is only done for the CC$0\pi$ samples.

  * •

The simulated non-CCQE events with true CC$0\pi$ topology are varied until the
prediction matches the data .

  * •

The change in the non-CCQE true CC$0\pi$ 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 $\nu_{\mu}$ and $\overline{\nu}_{\mu}$
samples to extract different weights for neutrino and anti-neutrino events,
where the effect was smaller for the $\overline{\nu}_{\mu}$ samples. The
weight extraction process and the resulting weights for neutrino events are
shown in Fig. 39.

Figure 39: The ND prediction after the fit to data for the FGD1 $\nu_{\mu}$
CC$0\pi$ 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.

The resulting distribution was used as the simulated data in the ND fit. As
shown in Fig. 40, the neutrino flux, 2p2h $\nu$ normalisation, and the $Q^{2}$
parameters shift to accommodate the model variation, notably present in the
higher $E_{\nu}$ region. This is largely expected, since non-CCQE events are
often from higher energy neutrinos due to the interaction cross section
increasing with $E_{\nu}$, whereas CCQE plateaus.

(a) FD $\nu$-mode $\nu_{\mu}$ flux parameters

(b) CC$0\pi$ cross-section parameters

(c) $Q^{2}$-dependent CCQE parameters

Figure 40:  Constraints on the $\nu$-mode $\nu_{\mu}$ flux (top) and CC$0\pi$
cross-section 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 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.
Fig. 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 $\Delta{}m^{2}_{32}$ and 21.3% on
$\sin^{2}\theta_{23}$ with respect to the size of the systematic uncertainty,
$\sigma_{syst}$. This is equivalent to an uncertainty inflation of $0.69\times
10^{-5}\leavevmode\nobreak\ \text{eV}^{2}$. This simulated data set has an
increased sensitivity to $\delta_{\scriptscriptstyle\mathrm{CP}}$, and no
significant bias in the best-fit value is observed.

Figure 41: Comparison of confidence regions in
$\sin^{2}\theta_{23}-\Delta{}m^{2}_{32}$ (top) for the expected sensitivity
fit (blue) and the non-CCQE simulated data study (red), and the
$\Delta\chi^{2}$ as a function of $\delta_{\scriptscriptstyle\mathrm{CP}}$
(bottom).

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.

###  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
Stowell:2019zsh . The effect was particularly pronounced at low 4-momentum
transfer, where a suppression of the cross section by $\sim 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 T2K:2019yqu  sees little need for a
suppression with recent NEUT versions Avanzini:2021qlx . Furthermore, MINERvA
sees the largest suppression in CC$1\pi^{0}$ selections, which are barely
selected in T2K’s FD. The CC$1\pi^{+}$ selection, 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
CC$0\pi$ interaction parameters. This is because they are largely constrained
by the dominant CC$0\pi$ selections at the ND, which are barely affected by
the simulated data. $C_{5}^{A}$ was pulled down, which acts to decrease the
overall SPP cross section, and $M_{A}^{RES}$ 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_{\mu}-\cos\theta_{\mu}$ to SPP events in the CC$0\pi$ selection. The
ratio of the prediction versus the simulated data at the ND is shown in Fig.
42, where the CC$1\pi$ sample is over-predicted by 20% in certain
$p_{\mu}-\cos\theta_{\mu}$ 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 CC$0\pi$ selections,
since SPP is a sub-leading contributor.

Figure 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 CC$1\pi$ sample.

The prediction for the SPP-dominated 1R$e$1d$e$ 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 1R$e$1d$e$ sample has
low statistics and occupies higher $E_{\nu}$ compared to the dominant
$\nu$-mode 1R$e$ selection, which provides most of the $\nu_{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.

Figure 43: The predictions for the 1R$e$1d$e$ selection at the FD using the
unchanged model (green) and including the MINERvA pion-suppression (blue). The
red bands show the propagated ND prediction from the fit to the simulated data
set, with the 1$\sigma$ uncertainty bands.

###  B.3 Summary of biases

As detailed in Sec. 9, a selection of the largest contributors to the bias
were used to inflate the confidence intervals for the two oscillation
parameters, $\Delta m^{2}$ and $\delta_{\scriptscriptstyle\mathrm{CP}}$. The
down-selection was used to avoid double-counting similar physics effects—for
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.

Tab. 18 presents the full bias table for the three oscillation parameters
$\sin^{2}\theta_{23}$, $\Delta{}m^{2}_{32}$, and
$\delta_{\scriptscriptstyle\mathrm{CP}}$. As with previous T2K analyses,
$\delta_{\scriptscriptstyle\mathrm{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 $\delta_{\scriptscriptstyle\mathrm{CP}}$ for
each simulated data study is shown in Tab. 19, where the shifts are dominated
by the three largest simulated data studies: the pion SI, the non-CCQE, and
data-driven pion tune.

Simulated data set | Relative to | $\sin^{2}\theta_{23}$ | $\Delta{}m^{2}_{32}$ | $\delta_{\scriptscriptstyle\mathrm{CP}}$
---|---|---|---|---
CCQE 3-comp nom. | Total | 1.0% | 0.4% | 0.8%
Syst. | 2.5% | 1.1% | 3.1%
CCQE 3-comp high | Total | 1.3% | 0.7% | 0.3%
Syst. | 3.2% | 1.8% | 1.1%
CCQE 3-comp low | Total | 0.7% | 0.2% | 0.2%
Syst. | 1.7% | 0.6% | 0.8%
CCQE z-exp nom. | Total | 2.5% | 0.2% | 0.6%
Syst. | 6.1% | 0.6% | 2.2%
CCQE z-exp high | Total | 0.3% | 2.1% | 0.4%
Syst. | 0.7% | 5.7% | 1.7%
CCQE z-exp low | Total | 3.1% | 0.2% | 0.1%
Syst. | 7.5% | 0.6% | 0.6%
CCQE removal energy | Total | 0.0% | 4.8% | 1.3%
Syst. | 0.0% | 13.4% | 5.2%
Non-CCQE | Total | 8.7% | 11.8% | 1.7%
Syst. | 21.3% | 32.7% | 6.9%
2p2h Martini | Total | 0.7% | 2.7% | 0.4%
Syst. | 1.6% | 7.3% | 1.6%
MINERvA pion tune | Total | 2.9% | 2.5% | 1.0%
Syst. | 7.2% | 6.8% | 3.5%
Data-driven pion | Total | 4.7% | 6.5% | 1.0%
Syst. | 11.6% | 17.9% | 3.9%
Pion SI | Total | 0.7% | 20.8% | 1.0%
Syst. | 1.9% | 57.8% | 4.6%
Table 18: Biases on the main oscillation parameters for each simulated data set, calculated as the shift in the middle of the $1\sigma$ confidence interval relative to the overall uncertainty from systematic sources (“Syst.”) and the total (“Total”) to one decimal place. Simulated data set | Change to 90% C.L. of $\delta_{\scriptscriptstyle\mathrm{CP}}$
---|---
-3.01 | -0.52
CCQE 3-comp nom. | 0.04 | 0.02
CCQE 3-comp high | 0.05 | 0.03
CCQE 3-comp low | 0.04 | 0.03
CCQE z-exp nom. | 0.01 | 0.01
CCQE z-exp high | 0.05 | 0.04
CCQE z-exp low | 0.00 | 0.00
CCQE removal energy | 0.00 | 0.02
Non-CCQE | 0.06 | 0.09
2p2h Martini | 0.04 | 0.04
MINERvA pion tune | 0.05 | 0.04
Data-driven pion | 0.07 | 0.04
Pion SI | 0.00 | 0.01
Table 19: Shifts of the 90% confidence interval boundaries of
$\delta_{\scriptscriptstyle\mathrm{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)
$\delta_{\scriptscriptstyle\mathrm{CP}}$ interval edge from the data fit. Only
the absolute size of the shift is taken into account.

##  Appendix C Two-dimensional Jarlskog credible regions

The measurements of $\delta_{\scriptscriptstyle\mathrm{CP}}$ and
$\sin^{2}\theta_{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. 25. T2K’s preference for values of
$\theta_{23}\sim 45\degree$ maximises the factor of
$\sin{\theta_{23}}\cos{\theta_{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
$\delta_{\scriptscriptstyle\mathrm{CP}}\approx-90\degree$ maximises the
$\sin\delta_{\scriptscriptstyle\mathrm{CP}}$ factor and selects negative
values of $J$. Fig. 45 illustrates the influence of the
$\delta_{\scriptscriptstyle\mathrm{CP}}$ posterior with values of
$\delta_{\scriptscriptstyle\mathrm{CP}}$ further away from maximal CP-
violation mapping to smaller preferred magnitudes of $J$.

Figure 44: 1$\sigma$, 2$\sigma$, and 3$\sigma$ credible regions in
$J-\sin^{2}\theta_{23}$ space extracted from the Bayesian analysis discussed
in Sec. 8.1. Figure 45: 1$\sigma$, 2$\sigma$, and 3$\sigma$ credible regions
in $J-\delta_{\scriptscriptstyle\mathrm{CP}}$ space extracted from the
Bayesian analysis discussed in Sec. 8.1.

##  Appendix D Additional critical values for $\sin^{2}\theta_{23}$ and
$\delta_{\scriptscriptstyle\mathrm{CP}}$

Critical values for the Feldman–Cousins confidence intervals and regions are
shown in Fig. 46 for one-dimensional fits in
$\delta_{\scriptscriptstyle\mathrm{CP}}$ and $\sin^{2}\theta_{23}$, and in
Fig. 47 for two-dimensional fits of
$\sin^{2}\theta_{23}-\delta_{\scriptscriptstyle\mathrm{CP}}$. In both cases,
the upper limit of the $1\sigma$ toy-statistical uncertainty interval on the
critical values is conservatively used for the confidence interval
computation. Compared to the asymptotic values from Wilks’s theorem, the
physical boundaries at $\delta_{\scriptscriptstyle\mathrm{CP}}=\pm\pi/2$ and
$\sin^{2}\theta_{23}\approx(2\cos^{2}\theta_{13})^{-1}\approx 0.513$ tend to
pull the critical values down, whereas degeneracies of
$\delta_{\scriptscriptstyle\mathrm{CP}}$, $\sin^{2}\theta_{23}$ around these
physical boundaries, and mass ordering, pull the critical values higher.
Further differences are also caused by the assumed distribution of true
oscillation parameter values for the nuisance oscillation parameters, i.e.
those not plotted here. These are thrown from a distribution close to their
posterior distributions from the fit to data, and the MC prediction with true
oscillation probabilities set to the T2K best-fit point is used as the “data”.
This is chosen instead of the actual data to avoid including the effects of
statistical fluctuations in data.

(a) Normal ordering, $\delta_{\scriptscriptstyle\mathrm{CP}}$.

(b) Inverted ordering, $\delta_{\scriptscriptstyle\mathrm{CP}}$.

(c) Normal ordering, $\sin^{2}\theta_{23}$.

(d) Inverted ordering, $\sin^{2}\theta_{23}$.

Figure 46: FC-corrected critical $\Delta\chi^{2}$ values for 1D-fits in
$\delta_{\scriptscriptstyle\mathrm{CP}}$ and $\sin^{2}\theta_{23}$, computed
at true values indicated with black dots and linearly interpolated in between.
The error bands show the $1\sigma$ toy-statistical uncertainty (binomial
confidence interval). The dotted lines show the asymptotic values from Wilks’s
theorem for reference.

(a) Normal ordering, 68%.

(b) Inverted ordering, 68%.

(c) Normal ordering, 90%.

(d) Inverted ordering, 90%.

Figure 47: 68% and 90% FC-corrected critical $\Delta\chi^{2}$ values for
2D-fits of $\sin^{2}\theta_{23}-\delta_{\scriptscriptstyle\mathrm{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.

##  Appendix E Systematic uncertainty parameters from the analysis of ND data

The parameters from the uncertainty model after the fit to ND data are
presented here. The carbon-only uncertainties are not propagated to the FD, as
no carbon is present there. The values before and after the fit to ND data,
and the uncertainty from the diagonal entry of the covariance matrix, are
shown. The values are tabulated in Tab. 20 and Tab. 21 for the flux
parameters, and Tab. 22 for the cross-section parameters. The fit results are
from the gradient-descent fitter, and the uncertainty is evaluated by scanning
around the best-fit point. These parameters are only valid in the context of
the exact model used in this analysis, and the context shouldn’t be
interpreted as an actual global constraint on the parameters.

Flux parameter | Pre-fit | Post-fit
---|---|---
FD $\nu_{\mu}$ [0.0, 0.4] | $1.00\pm 0.07$ | $1.11\pm 0.05$
FD $\nu_{\mu}$ [0.4, 0.5] | $1.00\pm 0.06$ | $1.10\pm 0.04$
FD $\nu_{\mu}$ [0.5, 0.6] | $1.00\pm 0.05$ | $1.08\pm 0.03$
FD $\nu_{\mu}$ [0.6, 0.7] | $1.00\pm 0.05$ | $1.07\pm 0.03$
FD $\nu_{\mu}$ [0.7, 1.0] | $1.00\pm 0.07$ | $1.06\pm 0.04$
FD $\nu_{\mu}$ [1.0, 1.5] | $1.00\pm 0.07$ | $1.03\pm 0.04$
FD $\nu_{\mu}$ [1.5, 2.5] | $1.00\pm 0.06$ | $1.03\pm 0.04$
FD $\nu_{\mu}$ [2.5, 3.5] | $1.00\pm 0.07$ | $1.01\pm 0.04$
FD $\nu_{\mu}$ [3.5, 5.0] | $1.00\pm 0.09$ | $0.98\pm 0.04$
FD $\nu_{\mu}$ [5.0, 7.0] | $1.00\pm 0.10$ | $0.92\pm 0.04$
FD $\nu_{\mu}$ [7.0, 30.0] | $1.00\pm 0.12$ | $0.91\pm 0.04$
FD $\overline{\nu}_{\mu}$ [0.0, 0.7] | $1.00\pm 0.09$ | $1.06\pm 0.08$
FD $\overline{\nu}_{\mu}$ [0.7, 1.0] | $1.00\pm 0.06$ | $1.04\pm 0.05$
FD $\overline{\nu}_{\mu}$ [1.0, 1.5] | $1.00\pm 0.07$ | $1.04\pm 0.06$
FD $\overline{\nu}_{\mu}$ [1.5, 2.5] | $1.00\pm 0.08$ | $1.05\pm 0.07$
FD $\overline{\nu}_{\mu}$ [2.5, 30.0] | $1.00\pm 0.08$ | $1.04\pm 0.06$
FD $\nu_{e}$ [0.0, 0.5] | $1.00\pm 0.06$ | $1.09\pm 0.04$
FD $\nu_{e}$ [0.5, 0.7] | $1.00\pm 0.05$ | $1.08\pm 0.04$
FD $\nu_{e}$ [0.7, 0.8] | $1.00\pm 0.05$ | $1.06\pm 0.04$
FD $\nu_{e}$ [0.8, 1.5] | $1.00\pm 0.06$ | $1.04\pm 0.04$
FD $\nu_{e}$ [1.5, 2.5] | $1.00\pm 0.08$ | $1.00\pm 0.04$
FD $\nu_{e}$ [2.5, 4.0] | $1.00\pm 0.09$ | $0.98\pm 0.04$
FD $\nu_{e}$ [4.0, 30.0] | $1.00\pm 0.09$ | $0.98\pm 0.05$
FD $\overline{\nu}_{e}$ [0.0, 2.5] | $1.00\pm 0.10$ | $1.02\pm 0.09$
FD $\overline{\nu}_{e}$ [2.5, 30.0] | $1.00\pm 0.13$ | $1.09\pm 0.11$
Table 20: FD $\nu$-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_{\nu}$ in units of GeV for each parameter. Flux parameter | Pre-fit | Post-fit
---|---|---
FD $\nu_{\mu}$ [0.0, 0.7] | $1.00\pm 0.09$ | $1.11\pm 0.06$
FD $\nu_{\mu}$ [0.7, 1.0] | $1.00\pm 0.06$ | $1.07\pm 0.05$
FD $\nu_{\mu}$ [1.0, 1.5] | $1.00\pm 0.06$ | $1.07\pm 0.04$
FD $\nu_{\mu}$ [1.5, 2.5] | $1.00\pm 0.07$ | $1.07\pm 0.04$
FD $\nu_{\mu}$ [2.5, 30.0] | $1.00\pm 0.07$ | $1.02\pm 0.04$
FD $\overline{\nu}_{\mu}$ [0.0, 0.4] | $1.00\pm 0.07$ | $1.09\pm 0.05$
FD $\overline{\nu}_{\mu}$ [0.4, 0.5] | $1.00\pm 0.06$ | $1.09\pm 0.04$
FD $\overline{\nu}_{\mu}$ [0.5, 0.6] | $1.00\pm 0.06$ | $1.07\pm 0.04$
FD $\overline{\nu}_{\mu}$ [0.6, 0.7] | $1.00\pm 0.05$ | $1.06\pm 0.03$
FD $\overline{\nu}_{\mu}$ [0.7, 1.0] | $1.00\pm 0.08$ | $1.09\pm 0.04$
FD $\overline{\nu}_{\mu}$ [1.0, 1.5] | $1.00\pm 0.08$ | $1.06\pm 0.04$
FD $\overline{\nu}_{\mu}$ [1.5, 2.5] | $1.00\pm 0.06$ | $1.01\pm 0.04$
FD $\overline{\nu}_{\mu}$ [2.5, 3.5] | $1.00\pm 0.07$ | $1.01\pm 0.05$
FD $\overline{\nu}_{\mu}$ [3.5, 5.0] | $1.00\pm 0.09$ | $0.95\pm 0.06$
FD $\overline{\nu}_{\mu}$ [5.0, 7.0] | $1.00\pm 0.09$ | $0.95\pm 0.06$
FD $\overline{\nu}_{\mu}$ [7.0, 30.0] | $1.00\pm 0.12$ | $0.93\pm 0.09$
FD $\nu_{e}$ [0.0, 2.5] | $1.00\pm 0.09$ | $1.03\pm 0.07$
FD $\nu_{e}$ [2.5, 30.0] | $1.00\pm 0.08$ | $1.03\pm 0.07$
FD $\overline{\nu}_{e}$ [0.0, 0.5] | $1.00\pm 0.06$ | $1.08\pm 0.04$
FD $\overline{\nu}_{e}$ [0.5, 0.7] | $1.00\pm 0.05$ | $1.07\pm 0.04$
FD $\overline{\nu}_{e}$ [0.7, 0.8] | $1.00\pm 0.06$ | $1.06\pm 0.04$
FD $\overline{\nu}_{e}$ [0.8, 1.5] | $1.00\pm 0.06$ | $1.04\pm 0.04$
FD $\overline{\nu}_{e}$ [1.5, 2.5] | $1.00\pm 0.08$ | $1.01\pm 0.06$
FD $\overline{\nu}_{e}$ [2.5, 4.0] | $1.00\pm 0.09$ | $1.01\pm 0.07$
FD $\overline{\nu}_{e}$ [4.0, 30.0] | $1.00\pm 0.15$ | $1.09\pm 0.13$
Table 21: FD $\overline{\nu}$-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_{\nu}$ in units of GeV for each parameter. Parameter | Pre-fit | Post-fit | Comment
---|---|---|---
$M_{A}^{QE}\leavevmode\nobreak\ (\text{GeV}/c^{2})$ | $1.03\pm 0.06$ | $1.17\pm 0.04$ |
$Q^{2}<0.05\leavevmode\nobreak\ \text{GeV}^{2}$ | $1.00\pm\infty$ | $0.78\pm 0.05$ |  Norm. on true CCQE events in true $Q^{2}$.
$0.05<Q^{2}<0.10\leavevmode\nobreak\ \text{GeV}^{2}$ | $1.00\pm\infty$ | $0.89\pm 0.04$
$0.10<Q^{2}<0.15\leavevmode\nobreak\ \text{GeV}^{2}$ | $1.00\pm\infty$ | $1.03\pm 0.05$
$0.15<Q^{2}<0.20\leavevmode\nobreak\ \text{GeV}^{2}$ | $1.00\pm\infty$ | $1.03\pm 0.08$
$0.20<Q^{2}<0.25\leavevmode\nobreak\ \text{GeV}^{2}$ | $1.00\pm\infty$ | $1.09\pm 0.10$
$0.25<Q^{2}<0.50\leavevmode\nobreak\ \text{GeV}^{2}$ | $1.00\pm 0.11$ | $1.26\pm 0.06$
$0.50<Q^{2}<1.00\leavevmode\nobreak\ \text{GeV}^{2}$ | $1.00\pm 0.18$ | $1.14\pm 0.08$
$Q^{2}>1.00\leavevmode\nobreak\ \text{GeV}^{2}$ | $1.00\pm 0.40$ | $1.26\pm 0.14$
$\Delta E_{rmv}^{C}\nu$ (MeV) | $2.00\pm 6.00$ | $-2.38\pm 1.75$ |
$\Delta E_{rmv}^{C}\overline{\nu}$ (MeV) | $0.00\pm 6.00$ | $1.64\pm 1.93$ |
$\Delta E_{rmv}^{O}\nu$ (MeV) | $4.00\pm 6.00$ | $2.55\pm 3.08$ |
$\Delta E_{rmv}^{O}\overline{\nu}$ (MeV) | $0.00\pm 6.00$ | $-1.26\pm 3.19$ |
2p2h norm. $\nu$ | $1.00\pm\infty$ | $1.06\pm 0.15$ |
2p2h norm. $\overline{\nu}$ | $1.00\pm\infty$ | $0.72\pm 0.16$ |
2p2h norm. C$\rightarrow$O | $1.00\pm 0.20$ | $1.05\pm 0.15$ |
2p2h shape C | $0.00\pm 3.00$ | $0.97\pm 0.46$ |  -1 is non-$\Delta$-like, 0 is Nieves et al. Nieves:2011yp , +1 is $\Delta$-like.
2p2h shape O | $0.00\pm 3.00$ | $0.00\pm 0.17$
2p2h low-$E_{\nu}$ $\nu$ | $1.00\pm 1.00$ | $1.00\pm 1.00$ |  +1 is Nieves-like Nieves:2011yp , 0 is Martini-like Martini:2009uj . Not fit at ND.
2p2h high-$E_{\nu}$ $\nu$ | $1.00\pm 1.00$ | $1.00\pm 1.00$
2p2h low-$E_{\nu}$ $\overline{\nu}$ | $1.00\pm 1.00$ | $1.00\pm 1.00$
2p2h high-$E_{\nu}$ $\overline{\nu}$ | $1.00\pm 1.00$ | $1.00\pm 1.00$
$C^{A}_{5}$ | $0.96\pm 0.15$ | $0.98\pm 0.06$ |
$M_{A}^{RES}\leavevmode\nobreak\ (\text{GeV}/c^{2})$ | $1.07\pm 0.15$ | $0.79\pm 0.05$ |
$I_{1/2}$ non-res norm. low-$p_{\pi}$ $\overline{\nu}_{\mu}$ | $0.96\pm 0.96$ | $0.96\pm 0.96$ | Not fit at ND.
$I_{1/2}$ non-res norm. | $0.96\pm 0.40$ | $0.87\pm 0.23$ |
CC coh. C norm. | $1.00\pm 0.30$ | $0.61\pm 0.22$ |
CC coh. O norm. | $1.00\pm 0.30$ | $0.61\pm 0.22$ |
Coulomb corr. $\nu$ | $1.00\pm 0.02$ | $1.00\pm 0.02$ |
Coulomb corr. $\overline{\nu}$ | $1.00\pm 0.01$ | $1.00\pm 0.01$ |
$\nu_{e}$/$\nu_{\mu}$ norm. | $1.00\pm 0.03$ | $1.00\pm 0.03$ |  No ND selection, poorly constrained.
$\overline{\nu}_{e}$/$\overline{\nu}_{\mu}$ norm. | $1.00\pm 0.03$ | $1.00\pm 0.03$
CC Bodek-Yang on/off DIS | $0.00\pm 1.00$ | $1.04\pm 0.19$ |  +1 is B-Y supp. off, 0 is B-Y supp. on. Bodek:2003wc ; Bodek:2005de
CC Bodek-Yang on/off multi-$\pi$ | $0.00\pm 1.00$ | $-0.03\pm 0.18$
CC multiplicity multi-$\pi$ | $0.00\pm 1.00$ | $0.14\pm 0.71$ |  +1 is AGKY-like Yang:2009zx , 0 is NEUT-like.
CC misc. norm. | $1.00\pm 1.00$ | $2.28\pm 0.43$ |
CC DIS+multi-$\pi$ norm. $\nu$ | $1.00\pm 0.04$ | $1.06\pm 0.03$ |
CC DIS+multi-$\pi$ norm. $\overline{\nu}$ | $1.00\pm 0.07$ | $0.94\pm 0.06$ |
NC coh. norm. | $1.00\pm 0.30$ | $1.02\pm 0.30$ |  No ND selection, poorly constrained.
NC $1\gamma$ norm. | $1.00\pm 1.00$ | $1.00\pm 1.00$
NC other ND norm. | $1.00\pm 0.30$ | $1.66\pm 0.13$ | Not propagated to FD.
NC other FD norm. | $1.00\pm 0.30$ | $1.00\pm 0.30$ | Not fit at ND.
Pion FSI Quasi-Elastic | $1.00\pm 0.29$ | $0.83\pm 0.09$ |  Scaling of pion scattering probabilities relative to the constraint from external data PinzonGuerra:2018rju  .
Pion FSI Quasi-Elastic $p_{\pi}>500\leavevmode\nobreak\ \text{MeV}/c$ | $1.00\pm 0.47$ | $0.75\pm 0.16$
Pion FSI Inelastic | $1.00\pm 1.10$ | $1.71\pm 0.31$
Pion FSI Absorption | $1.00\pm 0.31$ | $1.19\pm 0.12$
Pion FSI Charge Exchange | $1.00\pm 0.44$ | $0.78\pm 0.34$
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 Sec. 5. Parameters without external constraints are
labelled with an uncertainty $\pm\infty$.