# proEQUIB: IDL Library for Plasma Diagnostics and Abundance Analysis

Ashkbiz Danehkar  Research Centre in Astronomy, Astrophysics and
Astrophotonics, Macquarie University, Sydney, NSW 2109, Australia  Harvard-
Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138,
USA

††margin:  DOI: 10.21105/joss.00899 Software • Review • Repository • Archive
Submitted: 30 June 2018
Published: 04 December 2018 License
Authors of papers retain copyright and release the work under a Creative
Commons Attribution 4.0 International License (CC-BY).

## Summary

The emission lines emitted from gaseous nebulae carry valuable information
about the physical conditions and chemical abundances of ionized gases in
these objects, as well as the interstellar reddening. We determine the
electron temperature, the electron density, and the ionic abundances from the
dereddened fluxes of _collisionally excited lines_ (CEL) and _recombination
lines_ (RL) identified in nebular spectra (see e.g. Danehkar, Parker, and
Ercolano (2013); Danehkar et al. (2014); Danehkar, Parker, and Steffen (2016);
Danehkar (2018b)).

proEQUIB is a library including several application programming interface
(API) functions developed in the Interactive Data Language (IDL), which can be
used to determine temperatures, densities, and chemical abundances from
emission lines of ionized nebulae. This IDL library can also be used by the
GNU Data Language (GDL) (Arabas et al. (2010); Coulais et al. (2010)), which
is a free and open-source alternative IDL compiler. This IDL/GDL package
employs the IDL library AtomNeb _Atomic Data for Ionized Nebulae_ (Danehkar
(2018a)), which contains collision strengths and transition probabilities for
collisional excitation calculations, and recombination coefficients for
recombination calculations. This package includes several API functions to
determine physical conditions and chemical abundances from CEL and RL, derive
interstellar extinctions from Balmer lines, and deredden the observed fluxes:

  * •

The API functions for the _CEL analysis_ were developed in the IDL programming
language based on the algorithm of the FORTRAN program EQUIB (Howarth and
Adams (1981); Howarth et al. (2016)), which calculates atomic level
populations and line emissivities in statistical equilibrium in multi-level
atoms for the given physical conditions. These API functions can be used to
determine the electron temperature, the electron density, and the ionic
abundances from the dereddened fluxes of _collisionally excited lines_ emitted
from ionized gaseous nebulae.

  * •

The API functions for the _RL analysis_ were developed in IDL according to the
algorithm of the recombination scripts by X. W. Liu and Y. Zhang included in
the FORTRAN program MOCASSIN (Ercolano et al. (2003); Ercolano, Barlow, and
Storey (2005)). These API functions can be used to determine the ionic
abundances from the dereddened fluxes of _recombination lines_ emitted from
ionized nebulae.

  * •

The API functions for the _reddening analysis_ were developed based on the
methods of the reddening functions from the Space Telescope Science Data
Analysis System (STSDAS) IRAF Package (Bushouse and Simon (1994); Shaw and
Dufour (1994)). These API functions can be employed to obtain _interstellar
extinctions_ for different reddening laws from the observed fluxes of Balmer
lines detected in nebular spectra, and deredden the measured fluxes of
emission lines.

proEQUIB has recently been used for plasma diagnostics and abundance analysis
of some planetary nebulae (Danehkar, Parker, and Steffen (2016); Danehkar
(2018b)). This IDL/GDL package heavily relies on the IDL Astronomy User’s
library (Landsman (1993); Landsman (1995)) and the IDL library AtomNeb
(Danehkar (2018a)). The API functions of this IDL library can easily be
utilized to generate spatially-resolved maps of extinction, temperature,
density, and chemical abundances from integral field spectroscopic
observations (see e.g. Danehkar, Parker, and Ercolano (2013); Danehkar et al.
(2014); Danehkar (2014)).

## Acknowledgements

A.D. acknowledges the receipt of a Macquarie University Research Excellence
Scholarship.

## References

Arabas, S., M. Schellens, A. Coulais, J. Gales, and P. Messmer. 2010. “GNU
Data Language (GDL) - a free and open-source implementation of IDL.” In _EGU
General Assembly Conference_ , 12:924. Geophysical Research Abstracts.

Bushouse, H., and B. Simon. 1994. “The IRAF/STSDAS Synthetic Photometry
Package.” In _Astronomical Data Analysis Software and Systems Iii_ , edited by
D. R. Crabtree, R. J. Hanisch, and J. Barnes, 61:339. Astronomical Society of
the Pacific Conference Series.

Coulais, A., M. Schellens, J. Gales, S. Arabas, M. Boquien, P. Chanial, P.
Messmer, et al. 2010. “Status of GDL - GNU Data Language.” In _Astronomical
Data Analysis Software and Systems Xix_ , edited by Y. Mizumoto, K.-I. Morita,
and M. Ohishi, 434:187. Astronomical Society of the Pacific Conference Series.

Danehkar, A. 2014. “Evolution of Planetary Nebulae with WR-type Central
Stars.” PhD thesis, Macquarie University, Australia. doi:10.5281/zenodo.47794.

Danehkar, A. 2018a. “AtomNeb: Atomic Data for Ionized Nebulae.” _The Journal
of Open Source Software_ submitted.

Danehkar, A. 2018b. “Bi-Abundance Ionisation Structure of the Wolf-Rayet
Planetary Nebula PB 8.” _Publications of the Astronomical Society of
Australia_ 35: e005. doi:10.1017/pasa.2018.1.

Danehkar, A., Q. A. Parker, and B. Ercolano. 2013. “Observations and three-
dimensional ionization structure of the planetary nebula SuWt 2.” _Monthly
Notices of the Royal Astronomical Society_ 434: 1513–30.
doi:10.1093/mnras/stt1116.

Danehkar, A., Q. A. Parker, and W. Steffen. 2016. “Fast, Low-ionization
Emission Regions of the Planetary Nebula M2-42.” _The Astronomical Journal_
151: 38. doi:10.3847/0004-6256/151/2/38.

Danehkar, A., H. Todt, B. Ercolano, and A. Y. Kniazev. 2014. “Observations and
three-dimensional photoionization modelling of the Wolf-Rayet planetary nebula
Abell 48.” _Monthly Notices of the Royal Astronomical Society_ 439: 3605–15.
doi:10.1093/mnras/stu203.

Ercolano, B., M. J. Barlow, and P. J. Storey. 2005. “The dusty MOCASSIN: fully
self-consistent 3D photoionization and dust radiative transfer models.”
_Monthly Notices of the Royal Astronomical Society_ 362: 1038–46.
doi:10.1111/j.1365-2966.2005.09381.x.

Ercolano, B., M. J. Barlow, P. J. Storey, and X.-W. Liu. 2003. “MOCASSIN: a
fully three-dimensional Monte Carlo photoionization code.” _Monthly Notices of
the Royal Astronomical Society_ 340: 1136–52.
doi:10.1046/j.1365-8711.2003.06371.x.

Howarth, I. D., and S. Adams. 1981. “Program EQUIB.” University College
London.

Howarth, I. D., S. Adams, R. E. S. Clegg, D. P. Ruffle, X.-W. Liu, C. J.
Pritchet, and B. Ercolano. 2016. “EQUIB: Atomic level populations and line
emissivities calculator.” Astrophysics Source Code Library, ascl:1603.005.

Landsman, W. B. 1993. “The IDL Astronomy User’s Library.” In _Astronomical
Data Analysis Software and Systems Ii_ , edited by R. J. Hanisch, R. J. V.
Brissenden, and J. Barnes, 52:246. Astronomical Society of the Pacific
Conference Series.

Landsman, W. B. 1995. “The IDL Astronomy User’s Library.” In _Astronomical
Data Analysis Software and Systems Iv_ , edited by R. A. Shaw, H. E. Payne,
and J. J. E. Hayes, 77:437. Astronomical Society of the Pacific Conference
Series.

Shaw, R. A., and R. J. Dufour. 1994. “The FIVEL Nebular Modelling Package in
STSDAS.” In _Astronomical Data Analysis Software and Systems Iii_ , edited by
D. R. Crabtree, R. J. Hanisch, and J. Barnes, 61:327. Astronomical Society of
the Pacific Conference Series.


reduce the uncertainties on the neutrino fluxes before oscillation, and on the
neutrino interaction model parameters. The best-fit values of the neutrino and
antineutrino flux and cross-section parameters together with their
correlations are extrapolated to SK, where the oscillation probabilities are
measured. In this section, the general model that describes the oscillations
between the three standard neutrino flavors, $\nu_{e}$, $\nu_{\mu}$ and
$\nu_{\tau}$, are described, with a focus on the relevant oscillation channels
for the T2K experiment. Finally, external measurements of neutrino oscillation
parameters, which are used as prior constraints in some fits, are summarized.

###  VII.1 Oscillation probabilities in the PMNS framework

As anticipated by Pontecorvo, the neutrino flavor eigenstates do not
correspond to mass eigenstates, but are linear superpositions of them:

$\displaystyle|\nu_{\alpha}\rangle=\sum_{k}U_{\alpha k}|\nu_{k}\rangle\hskip
28.45274pt\left(\alpha=e,\mu,\tau\right)$ (6)

where only the three active neutrinos, $k=1,2,3$ are considered. $U_{\alpha
k}$ is an element of the $3\times 3$ unitary matrix, called the Pontecorvo-
Maki-Nakagawa-Sakata (PMNS) matrix Maki:1962mu ; Pontecorvo:1967fh , which can
be parametrized as:

$\displaystyle U=$ $\displaystyle\left(\begin{array}[]{ccc}1&0&0\\\
0&c_{23}&s_{23}\\\
0&-s_{23}&c_{23}\end{array}\right)\left(\begin{array}[]{ccc}c_{13}&0&s_{13}e^{-i\delta_{CP}}\\\
0&1&0\\\
-s_{13}e^{i\delta_{CP}}&0&c_{13}\end{array}\right)\left(\begin{array}[]{ccc}c_{12}&s_{12}&0\\\
-s_{12}&c_{12}&0\\\ 0&0&1\end{array}\right)$ (16)

where $s_{ij}\equiv\sin\theta_{ij}$ and $c_{ij}\equiv\cos\theta_{ij}$,
$\theta_{ij}$ are the three mixing angles and $\delta_{CP}$ is the CP-
violating phase. The Majorana phases are neglected here as the three-flavor
oscillation probability is invariant under their rotation.

The neutrino oscillation probability is determined by six parameters: the
mixing angles $\theta_{12}$, $\theta_{13}$, $\theta_{23}$, which define the
amplitude of the oscillation probability; the differences in the squared
masses of the eigenstates, $\Delta m^{2}_{21}=m^{2}_{2}-m^{2}_{1}$ and $\Delta
m^{2}_{32}=m^{2}_{3}-m^{2}_{2}$, that define the oscillation frequency and the
position of the oscillation maxima as a function of $L/E$; and the CP-
violating phase, $\delta_{CP}$. While it is known that $\Delta m^{2}_{21}>0$,
the neutrino mass ordering (MO) and hence the sign of $\Delta m^{2}_{32}$ has
not yet been determined.

Effects due to charged-current coherent scattering of $\nu_{e}$ with the
electrons in the Earth matter can affect the neutrino event rate. The matter
effect potential is proportional to $+N_{e}$ for $\nu_{e}$ and $-N_{e}$ for
$\bar{\nu}_{e}$, where $N_{e}$ is the matter electron density.

Neglecting matter effects, the $\nu^{\hbox{\kern-8.0pt\raise
1.0pt\hbox{{\tiny(}{\lower 1.4pt\hbox{\bf--}}{\tiny)}}}}_{\mu}$ survival
probability can be written as

$\displaystyle P(\;\nu^{\hbox{\kern-8.0pt\raise 1.0pt\hbox{{\tiny(}{\lower
1.4pt\hbox{\bf--}}{\tiny)}}}}_{\mu}\rightarrow\nu^{\hbox{\kern-8.0pt\raise
1.0pt\hbox{{\tiny(}{\lower
1.4pt\hbox{\bf--}}{\tiny)}}}}_{\mu})=1-4\left(s^{2}_{12}c^{2}_{23}+s^{2}_{13}s^{2}_{23}c^{2}_{12}+2s_{12}s_{13}s_{23}c_{12}c_{23}\cos\delta_{CP}\right)s^{2}_{23}c^{2}_{13}\sin^{2}\phi_{31}$
$\displaystyle-4\left(c^{2}_{12}c^{2}_{23}+s^{2}_{13}s^{2}_{23}s^{2}_{12}-2s_{12}s_{13}s_{23}c_{12}c_{23}\cos\delta_{CP}\right)s^{2}_{23}c^{2}_{13}\sin^{2}\phi_{32}-4\left(s^{2}_{12}c^{2}_{23}+s^{2}_{13}s^{2}_{23}c^{2}_{12}+2s_{12}s_{13}s_{23}c_{12}c_{23}\cos\delta_{CP}\right)$
$\displaystyle\times\left(c^{2}_{12}c^{2}_{23}+s^{2}_{13}s^{2}_{23}s^{2}_{12}-2s_{12}s_{13}s_{23}c_{12}c_{23}\cos\delta_{CP}\right)\sin^{2}\phi_{21},$
(18)

where

$\displaystyle\phi_{kj}=\frac{\Delta m^{2}_{kj}L}{4E}.$ (19)

$|\Delta m^{2}_{32}|$ and $\sin^{2}2\theta_{23}$ can be determined from Eq.
18, but the $\theta_{23}$ octant, i.e., whether $\theta_{23}>\pi/4$ or
$\theta_{23}<\pi/4$, cannot be distinguished. At the T2K baseline and peak
neutrino energy, the impact of the matter effect is small, and a full
calculation including them would only modify the
$P(\nu_{\mu}\rightarrow\nu_{\mu})$ and
$P(\bar{\nu}_{\mu}\rightarrow\bar{\nu}_{\mu})$ probabilities by $\sim$0.1%.

The $\nu^{\hbox{\kern-8.0pt\raise 1.0pt\hbox{{\tiny(}{\lower
1.4pt\hbox{\bf--}}{\tiny)}}}}_{e}$ appearance probability, approximated to
first order in the matter effect (see Arafune:1997hd ), is

$\displaystyle P(\;\nu^{\hbox{\kern-8.0pt\raise 1.0pt\hbox{{\tiny(}{\lower
1.4pt\hbox{\bf--}}{\tiny)}}}}_{\mu}\rightarrow\nu^{\hbox{\kern-8.0pt\raise
1.0pt\hbox{{\tiny(}{\lower
1.4pt\hbox{\bf--}}{\tiny)}}}}_{e})=4c^{2}_{13}s^{2}_{13}s^{2}_{23}\sin^{2}\phi_{31}\left(1+\frac{2\alpha}{\Delta
m^{2}_{31}}\left(1-2s^{2}_{13}\right)\right)+8c^{2}_{13}s_{12}s_{13}s_{23}\left(c_{12}c_{23}\cos\delta_{CP}-s_{12}s_{13}s_{23}\right)\cos\phi_{23}\sin\phi_{31}\sin\phi_{21}$
$\displaystyle\accentset{(+)}{-}8c^{2}_{13}c_{12}c_{23}s_{12}s_{13}s_{23}\sin\delta_{CP}\sin\phi_{32}\sin\phi_{31}\sin\phi_{21}\
+4s^{2}_{12}c^{2}_{13}\left(c^{2}_{12}c^{2}_{23}+s^{2}_{12}s^{2}_{23}s^{2}_{13}-2c_{12}c_{23}s_{12}s_{23}s_{13}\cos\delta_{CP}\right)\sin^{2}\phi_{21}$
$\displaystyle-8c^{2}_{13}s^{2}_{13}s^{2}_{23}\left(1-2s^{2}_{13}\right)\frac{{\accentset{(-)}{+}}\alpha
L}{4E}\cos\phi_{32}\sin\phi_{31},$ (20)

where $``(-)$” corresponds to the change of sign required for antineutrino
oscillations. $\rho$ is the average density of the Earth matter through which
the neutrinos travel and
$\alpha\left[\text{eV}^{2}/\text{c}^{4}\right]=7.56\times
10^{-5}\rho\left[\text{g}/\text{cm}^{3}\right]E_{\nu}\left[\text{GeV}\right]$.
In all oscillation analyses described here, the exact three-flavor oscillation
probabilities are used, rather than the approximations shown in this section.

As well as being proportional to $\sin^{2}\theta_{13}$, the first, leading,
term of Eq. 20 is proportional to $\sin^{2}\theta_{23}$, which makes the
appearance probability (in contrast to the disappearance probability)
sensitive to the octant of $\theta_{23}$. The precise determination of
$\theta_{23}$ and the solution of the octant degeneracy are important for the
determination of $\delta_{\mathrm{CP}}$. Because of the precise measurement of
$\theta_{13}$ from reactor experiments (see Tab. 18), the largest uncertainty
in the oscillation model which affects the determination of
$\delta_{\mathrm{CP}}$ is due to the uncertainties on the value of
$\theta_{23}$.

By looking at the other terms in Eq. 20, it is easy to understand why the
$\nu^{\hbox{\kern-8.0pt\raise 1.0pt\hbox{{\tiny(}{\lower
1.4pt\hbox{\bf--}}{\tiny)}}}}_{e}$ appearance transition is the golden channel
in the search for CP violation. The term containing $\sin\delta_{CP}$ in Eq.
20 has an opposite sign for neutrinos and antineutrinos, and is the only term
that can violate CP symmetry. Since all the mixing angles have been measured
as nonzero, the requirement for violating CP symmetry in neutrino oscillations
is $\sin\delta_{CP}\neq 0$, i.e. $\delta_{CP}\neq 0,\pi$.

In T2K, the term proportional to $\sin\delta_{CP}$ can change the appearance
probability by as much as $\pm 30\%$. This means that an extreme value of
$\delta_{\mathrm{CP}}$ such as $\delta_{\mathrm{CP}}=-\pi/2$ would increase
(decrease) the $\nu_{e}$ ($\overline{\nu}_{e}$) by 30% with respect to
$\delta_{\mathrm{CP}}=0$. In Eq. 20 there are also other terms containing
$\cos\delta_{CP}$. While these terms do not violate CP symmetry, they change
the shape of the appearance spectrum and are important for a precise
determination of the value of $\delta_{\mathrm{CP}}$.

In summary, at first order, the term in $\sin\delta_{CP}$ defines the
normalization of the $\nu^{\hbox{\kern-8.0pt\raise 1.0pt\hbox{{\tiny(}{\lower
1.4pt\hbox{\bf--}}{\tiny)}}}}_{\mu}\rightarrow\nu^{\hbox{\kern-8.0pt\raise
1.0pt\hbox{{\tiny(}{\lower 1.4pt\hbox{\bf--}}{\tiny)}}}}_{e}$ spectrum, while
the term in $\cos\delta_{CP}$ provides information on the shape and allows the
discrimination of the CP-conserving phases $\delta_{CP}=0$ and $\pi$.

As mentioned above the appearance probability is also affected by the matter
effects. Matter effects can enhance either ${\nu}_{\mu}\rightarrow{\nu}_{e}$,
or $\bar{\nu}_{\mu}\rightarrow\bar{\nu}_{e}$ depending on the sign of $\Delta
m^{2}_{32}$. For instance if the mass ordering is normal (inverted)
${\nu}_{\mu}\rightarrow{\nu}_{e}$ oscillations are enhanced (suppressed) while
$\bar{\nu}_{\mu}\rightarrow\bar{\nu}_{e}$ oscillations are suppressed
(enhanced).

With its low neutrino energy and baseline of 295 km, the T2K experiment is not
sensitive to the matter effect. The two possible solutions for the mass
ordering (MO) change the appearance probability by approximately $\pm 10\%$
and a definitive measurement of the MO cannot be made by exploiting the matter
effects in T2K beam. The interplay between effects due to
$\delta_{\mathrm{CP}}$ and effects due to the MO on the
$\nu^{\hbox{\kern-8.0pt\raise 1.0pt\hbox{{\tiny(}{\lower
1.4pt\hbox{\bf--}}{\tiny)}}}}_{e}$-appearance probabilities are shown in Fig.
27. Each ellipse shows the effect due to $\delta_{\mathrm{CP}}$ on the two
appearance probabilities and the two ellipses represent the two possible
solution for the mass ordering. The two regions overlap for more than half of
the $\delta_{\mathrm{CP}}$ range, indicating the large degeneracy among the MO
and $\delta_{\mathrm{CP}}$ at the T2K baseline.

Figure 27: The oscillation probability for ${\nu}_{\mu}\rightarrow{\nu}_{e}$
and $\bar{\nu}_{\mu}\rightarrow\bar{\nu}_{e}$ transitions as a function of
$\delta_{CP}$ and MO for $E_{\nu}=0.6$ GeV and a baseline of 295 km. The other
oscillation parameters are fixed to the values shown in Tab. 13. The dashed
line indicates
$P({\nu}_{\mu}\rightarrow{\nu}_{e})=P(\bar{\nu}_{\mu}\rightarrow\bar{\nu}_{e})$.

The goal of T2K is the precise measurement of $\theta_{13}$, $\theta_{23}$,
$|\Delta m^{2}_{32}|$ and to search for evidence of CP violation. Since T2K
has excellent sensitivity to the absolute value of $\Delta m^{2}_{32}$ but not
to its sign, the absolute value of $\Delta m^{2}_{32}$ and the MO are
considered as two separate parameters. As mentioned above, the violation of CP
symmetry would manifest itself as a difference between neutrino and
antineutrino oscillation probabilities. In order to find a hint of this
violation it is necessary to directly compare
${\nu}_{\mu}\rightarrow{\nu}_{e}$ and
$\bar{\nu}_{\mu}\rightarrow\bar{\nu}_{e}$ oscillation probabilities, as also
shown in Fig. 27, where the anti-correlation as a function of both
$\delta_{CP}$ and MO is visible. The measurement of the
$\nu^{\hbox{\kern-8.0pt\raise 1.0pt\hbox{{\tiny(}{\lower
1.4pt\hbox{\bf--}}{\tiny)}}}}_{\mu}$ survival probabilities is important for
determining $|\Delta m^{2}_{32}|$, i.e. the position of the first oscillation
maximum, and $\sin^{2}2\theta_{23}$, drastically reducing the uncertainty on
$\delta_{CP}$. The degeneracy between $\delta_{CP}$ and the $\theta_{23}$
octant can be solved by comparing $\nu_{e}$ and $\bar{\nu}_{e}$ appearance
probabilities. Indeed, in contrast to $\delta_{CP}$, $\sin^{2}\theta_{23}$ has
the same effect on both the ${\nu}_{\mu}\rightarrow{\nu}_{e}$ and the
$\bar{\nu}_{\mu}\rightarrow\bar{\nu}_{e}$ oscillation probabilities. For this
reason, T2K measures the $\nu_{\mu}\rightarrow\nu_{\mu}$
($\bar{\nu}_{\mu}\rightarrow\bar{\nu}_{\mu}$) disappearance and the
$\nu_{\mu}\rightarrow\nu_{e}$ ($\bar{\nu}_{\mu}\rightarrow\bar{\nu}_{e}$)
appearance probabilities in $\nu$\- and $\overline{\nu}$-mode simultaneously
in the analyses described in this paper. The joint analysis also allows all of
the correlations between the oscillation parameters to be properly taken into
account.

###  VII.2 External constraints on the oscillation parameters

For some of the oscillation parameters, Gaussian priors, based on external
experiments, are used in the oscillation analyses. These parameters are the
solar oscillation parameters, to which T2K is insensitive, and, in some cases,
the value of $\textrm{sin}^{2}\theta_{13}$ as measured by reactor experiments.
In the latter case, the weighted average of Daya Bay, RENO and Double Chooz as
reported in Ref. PDG2015  is used. The Gaussian priors used in the analysis
are shown in Tab. 18. For all the results shown in this paper, it will be
specified whether reactor constraints on $\textrm{sin}^{2}\theta_{13}$ are
used.

All the parameters that are measured by T2K are included in the oscillation
analyses with flat priors which extend beyond the allowed regions for those
parameters. The boundaries on the flat priors used in the analysis are
$0.3<\textrm{sin}^{2}\theta_{23}<0.7$, 2$<$$\Delta m^{2}_{32}$ ($\times
10^{-3}~{}eV^{2}/c^{4}$)$<$3, $-\pi<\delta_{\mathrm{CP}}<\pi$, and
$0<\textrm{sin}^{2}(2\theta_{13})<0.4$ (when no reactor measurements are
used).

Table 18: Gaussian priors used for the neutrino oscillation parameters, taken from Ref. PDG2015 . Oscillation parameter |  Mean and standard deviation
---|---
$\sin^{2}2\theta_{12}$ |  $0.846\pm 0.021$
$\Delta m^{2}_{21}$ |  $(7.53\pm 0.18)\times 10^{-5}~{}\text{eV}^{2}/\text{c}^{4}$
$\sin^{2}2\theta_{13}$ |  $0.085\pm 0.005$

##  VIII Oscillation analysis method

In this section, the oscillation analyses which have been developed to provide
oscillation parameter estimations are described: three different analyses have
been developed within frequentist and Bayesian frameworks.

For all analyses, the appearance and disappearance channels in both $\nu$\-
and $\overline{\nu}$-modes are analysed simultaneously, incorporating
constraints from several T2K samples and external experiments. The oscillation
analyses compare the rate and distribution of events in analysis bins defined
by combinations of the reconstructed neutrino energy, $E_{\nu}^{rec}$ (see
Equation 4), the reconstructed lepton momentum, and the reconstructed angle of
the outgoing lepton with model predictions. Systematic uncertainties are
characterized by systematic model parameters, which tune the model prediction,
and by a covariance matrix, which determines the sizes and correlations of
their uncertainties. The fit to the near detector data is applied as a
multivariate Gaussian penalty term to constrain flux and cross-section
uncertainties common to the near and far detectors.

For a given set of systematic and oscillation parameters, the expected
distribution is computed in the following way. The expected SK MC samples are
produced assuming nominal systematic parameters without oscillation. The
effect of oscillations are then included by applying the neutrino oscillation
probability to events defined by neutrino energy and interacting flavor. In
order to do this all the analyses use the full three-flavor oscillation
probabilities. The number of events is weighted by the oscillation probability
and for the systematic parameters, as will be described in Sec. VIII.4. Then
the expected number of events in each bin of the far detector samples are
extracted from a Poisson distribution, leading to a likelihood of the form

$-2\;\textrm{ln}\;\mathcal{L}(\vec{o},\vec{f})=2\sum_{i=0}^{N-1}\Big{[}n^{obs}_{i}\cdot\textrm{ln}(n^{obs}_{i}/n^{exp}_{i})+(n^{exp}_{i}-n^{obs}_{i})\Big{]},$
(21)

where $\vec{o}$ and $\vec{f}$ are particular choices of the vector of the
neutrino oscillation parameters (free in the fit), and of the nuisance
parameters. $n^{obs}_{i}$ is the observed number of events in the $i$th
analysis template bin, and $n^{exp}_{i}=n^{exp}_{i}(\vec{o},\vec{f})$ is the
expected number of events in the $i$th bin.

The method described in Ref. Cousins:1991qz  is used to project out the
nuisance parameters. The high dimensionality of the likelihood function is
solved by computing the marginal likelihood. This marginal likelihood for the
parameters of interest is found by integrating the product of the likelihood
function and the priors over all the parameters that are not of interest, and
defined as:

$\displaystyle\mathcal{L}_{marg}(\vec{o};x)$
$\displaystyle=\int_{F}\mathcal{L}(\vec{o},\vec{f};x)\pi{(\vec{f})}d\vec{f}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\mathcal{L}(\vec{o},\vec{f}_{i};x).$
(22)

This defines the likelihood which is a function only of the parameters of
interest $\vec{o}$, given a dataset, $x$, and a model prediction defined by a
set of systematic and oscillation parameters. The marginal likelihood of Eq.
22 can be used to perform both a hybrid frequentist-Bayesian analysis as well
as a full Bayesian analysis. In both cases the prior p.d.f. for flux, cross
section and SK detector parameters is a multivariate Gaussian, with mean $\mu$
and covariance $V$. All the analyses have been validated with simulated
datasets using different neutrino interaction models and different sets of
oscillation parameters for the generation of the simulated samples. Examples
of the validations will be shown in Sect. IX.

###  VIII.1 Frequentist analyses

Two analyses have been developed in the frequentist framework. The main
difference between the two analysis is that one uses 2D templates in
{$E_{\nu}^{rec}$, $\theta_{lep}$} for the e-like samples and the other uses
{$p_{lep}$, $\theta_{lep}$} templates. The reasons for these choices will be
described in Sec. VIII.3. Both analyses construct confidence intervals or
regions by fitting the SK dataset, using the likelihood ratio
$-2\Delta\ln\mathcal{L}(\vec{o})=-2\ln[\mathcal{L}(\vec{o};x)/\mathcal{L}_{max}(x)]$
as the test statistic, where $\mathcal{L}(\vec{o};x)$ is as defined in Eq. 22.
If the Gaussian approximation holds, Eq. 22 converges to a $\chi^{2}$
distribution, and the negative log-likelihood ratio is sometimes denoted as
$\Delta\chi^{2}(\vec{o})$.

Due to the difficulties in providing confidence intervals with high
dimensionality, a subset of the oscillation parameters is projected out by
using the marginalization method described above, with the priors shown in
Tab. 18.

Two approaches can be used in order to calculate the confidence intervals. In
the Gaussian approximation, where the parameter dependence is linear, the so-
called “constant $\Delta\chi^{2}$ method” provides the correct coverage: the
confidence region consists of the set of oscillation parameter values that
satisfies the condition $-2\Delta\ln\mathcal{L}(\vec{o})\leq
X_{crit}(\vec{o})$. The precomputed values of $X_{crit}(\vec{o})$ PDG2015 ,
that define the $\alpha\%$ confidence level for different dimensionality, do
not depend on the value of $\vec{o}$.

Since in the neutrino oscillation analysis the gaussian approximation
condition does not always hold, a toy MC method, recommended by Feldman and
Cousins feldmancousins , should be used to define the values of
$X_{crit}(\vec{o})$ as a function of $\vec{o}$, in order to provide the
correct coverage. An ensemble of 10,000 toy datasets was produced for each
point of the oscillation parameter grid, fine enough to describe the variation
of $X_{crit}(\vec{o})$ as a function of the free parameters. The flux, cross-
section and detector systematic parameters were sampled from the same
multivariate Gaussian with mean $\mu$ and covariance $V$ used as prior p.d.f.
in the marginalization method, as described above. The oscillation parameters
constrained by external measurements, such as $\sin^{2}\theta_{12}$, $\Delta
m^{2}_{21}$ and $\sin^{2}\theta_{13}$, were sampled by using the p.d.f.’s with
parameters given in Tab. 18. The remaining oscillation parameters were sampled
using the distribution of $-2\Delta\ln\mathcal{L}(\vec{o})$ obtained from the
fit of the T2K data. This method will be referred as “posterior-predictive”.

Since the Feldman-Cousins method is computationally intensive, it is used only
for the most important result of the oscillation analysis, i.e. the confidence
intervals as a function of $\delta_{CP}$ and MO. For all the other results the
constant $\Delta\chi^{2}$ method, which is more practical and provides good
coverage to a first order approximation, is used.

###  VIII.2 Bayesian analysis

The Bayesian analysis performs the integration of Eq. 22 using a Markov Chain
Monte Carlo (MCMC) method. The Metropolis-Hastings algorithm Hastings:1970a
is used to populate the space of the oscillation and nuisance parameters,
distributed according to the posterior probability density function.

This algorithm uses a weighted random walk to explore the parameter space with
a chain of many points, $X_{i}$. A new proposed step to $X_{i+1}$ is accepted
with a probability equal to the minimum between 1 and $P(X_{i+1})/P(X_{i})$,
i.e., it is accepted with 100% probability when $P(X_{i+1})>P(X_{i})$ and with
a probability equal to $P(X_{i+1})/P(X_{i})$ when $P(X_{i+1})<P(X_{i})$. When
a step to $X_{i+1}$ is rejected, the point $X_{i}$ is repeated in the chain,
and a a new step is proposed again from the previous step $X_{i}$. Typically,
a chain consists of $10^{6}-10^{7}$ steps. The population of accepted points
in parameter space thus constructed is distributed according to the posterior
probability density.

Credible intervals are built by selecting the highest density region in the
space of the parameters of interest, typically one or two, that contains
$\alpha\%$ of the points. A kernel density estimator (KDE) Mills:KDE ;
Cranmer:KDE , is used to extract the most probable values of the oscillation
parameters of interest, A KDE PDF is built by smearing the step points with a
Gaussian function which has a width inversely proportional to the number of
points. The most probable values of the parameters of interest were found by
maximizing the KDE PDF using MINUIT minuit .

###  VIII.3 Analysis templates

The three analyses use different far detector templates. All of them use the
reconstructed neutrino energy, $E_{\nu}^{rec}$, for the $\nu_{\mu}$ and
$\bar{\nu}_{\mu}$ candidate samples, while different templates are used for
the $\nu_{e}$ and $\bar{\nu}_{e}$ candidate samples. One of the two
frequentist analyses and the MCMC Bayesian analysis use $E_{\nu}^{rec}$ and
the angle between the lepton and the neutrino beam direction, $\theta_{lep}$.
The choice of these two parameters is due to the fact that $E_{\nu}^{rec}$ is
important to infer the oscillation probability function, which depends on the
true neutrino energy. $\theta_{lep}$ is used because it provides additional
separation between $\nu$ and $\bar{\nu}$ candidate events. The $\nu$ and
$\bar{\nu}$ cross sections differ as a function of $Q^{2}$, and leptons
produced in $\nu$ interactions tend to be emitted at larger angles than those
produced in $\overline{\nu}$ interactions.

As was shown in Tab. 14, while the number of
$\bar{\nu}_{\mu}\rightarrow\bar{\nu}_{e}$ in $\nu$-mode is negligible, in
$\overline{\nu}$-mode, the number of oscillated $\nu_{e}$ is roughly 30% of
the number of oscillated $\overline{\nu}_{e}$ and thus it is crucial to use
templates which are able to distinguish between the two cases. In Fig. 28, the
expected {$E_{\nu}^{rec}$, $\theta_{lep}$} distribution in the single ring
e-like sample in $\overline{\nu}$-mode for the signal
($\bar{\nu}_{\mu}\rightarrow\bar{\nu}_{e}$) and for the oscillated $\nu_{e}$
background is shown.

In Fig. 29 the 2D templates for the three $\nu_{e}$ candidate samples are
shown. The $\nu_{e}$ CC1$\pi^{+}$ sample has no events with reconstructed
energy below 400 MeV due to the $\Delta^{++}$ production requirement in Eq. 5
used to reconstruct the neutrino energy.

Figure 28:  Expected {$E_{\nu}^{rec}$, $\theta_{lep}$} distributions for the
signal $\bar{\nu}_{\mu}\rightarrow\bar{\nu}_{e}$ (top) and wrong-sign
background $\nu_{\mu}\rightarrow\nu_{e}$ (bottom) in the e-like CCQE-like
selection in $\overline{\nu}$-mode. The superimposed grey dots correspond to
the data. The expectation is based on the parameters of Tab. 13.

Figure 29:  {$E_{\nu}^{rec}$,$\theta_{lep}$} templates for the $\nu_{e}$ CCQE-
like (top), $\nu_{e}$ CC1$\pi^{+}$ (middle) in $\nu$-mode and $\bar{\nu}_{e}$
CCQE-like (bottom) in $\overline{\nu}$-mode candidate samples. Both signal and
background events are included in the expected distributions based on the
oscillation parameters of Tab. 13. The superimposed grey dots correspond to
the data.

The second frequentist analysis also uses 2D templates but instead of using
$E_{\nu}^{rec}$, it uses lepton momentum, $p_{lep}$, and $\theta_{lep}$. This
choice provides similar $\nu$–$\overline{\nu}$ separation as {$p_{lep}$,
$\theta_{lep}$} templates are directly related to $E_{\nu}^{rec}$. Indeed, a
particular value of $E_{\nu}^{rec}$ corresponds to a slice in the {$p_{lep}$,
$\theta_{lep}$} phase space. For this reason, {$E_{\nu}^{rec}$,
$\theta_{lep}$} and {$p_{lep}$, $\theta_{lep}$} analysis templates provide
very similar sensitivities to the oscillation parameters.

Finally, another difference among the three analyses is that the Bayesian
analysis performs a simultaneous fit of the ND280 and SK datasets in order to
validate the extrapolation of the flux, cross-section and detector systematic
parameters from the near to far detector. The other two analyses perform a fit
of the far detector data using the priors obtained by the near detector data
fit for the flux and cross section systematic parameters. The two methods show
very good agreement in the estimation of flux and cross-section parameters as
was shown in Sec. V.

###  VIII.4 Systematic uncertainties on the oscillation analyses

The systematic uncertainties of the oscillation analysis are separated into
four different categories: flux, cross section, and SK detector and momentum
scale parameters. The flux and many of the cross-section parameters (mentioned
in Tab. 9) are correlated by the ND280 measurements as shown in Sec. V,
effectively reducing the systematic uncertainties of the measurements at SK.

The flux parameters assign the systematic uncertainties on the neutrino and
antineutrino flux. These parameters are applied as normalization factors for
different neutrino flavor subcategories and true energy bins. The flux
uncertainties at SK are constrained by the ND280 measurements.

The cross-section parameters are applied based on the true neutrino
interaction category of each event. Most of these parameters are normalization
factors, but some also contain shape information, modifying the kinematic
distributions in $p_{lep}-\theta_{lep}$ and $E_{\nu}^{rec}$ probability
density functions for corresponding oscillation analysis samples.

There are several categories of SK uncertainty: detector efficiency, final
state interaction (FSI), secondary interaction (SI), and photo-nuclear (PN)
effect uncertainties. The update on the SK detector efficiency uncertainty
related to the addition of the CC1$\pi$ e-like sample is described in Sec.
VI.3. In addition, the FSI, SI, and PN effect uncertainties on the
CC1$\pi^{+}$ e-like sample are evaluated using the same method introduced in
Ref. Abe:2015awa . The increased uncertainties for this sample is mainly due
to the larger backgrounds affecting it, and to larger uncertainties on the
pion FSI and SI as shown in Tab. 19, where the SK systematic uncertainties are
separated into the SK detector and the FSI+SI+PN contributions.

The SK energy scale uncertainty is applied independently from other
parameters. The energy scale uncertainty is applied as a normalization of
$E_{\nu}^{rec}$ for each event, which may vary the total event rate by
shifting the events into the cut regions of the visible energy and
reconstructed neutrino energy selection criteria. The SK energy scale
uncertainty is estimated to be 2.4%.

The effect of the systematic uncertainties on the predicted event rates of the
$\nu$\- and $\overline{\nu}$-mode samples are summarized in Tab. 19 and Tab.
20 respectively. The $1\sigma$ uncertainties are obtained by throwing large
sets of toy experiments, varying only the selected systematic parameters for
each experiment, and calculating the relative uncertainties from the
distributions of the event rates. The anti-correlations between flux and cross
section parameters reduce the systematic uncertainties when both these sources
are taken into account.

Table 19: Effect of 1$\sigma$ variation of the systematic uncertainties on the predicted event rates of the $\nu$-mode samples. Source of uncertainty |  $\nu_{e}$ CCQE-like | $\nu_{\mu}$ |  $\nu_{e}$ CC1$\pi^{+}$
---|---|---|---
| $\delta N/N$ | $\delta N/N$ | $\delta N/N$
Flux | 3.7% | 3.6% | 3.6%
(w/ ND280 constraint) |  |  |
Cross section | 5.1% | 4.0% | 4.9%
(w/ ND280 constraint) |  |  |
Flux+cross-section |  |  |
(w/o ND280 constraint) | 11.3% | 10.8% | 16.4%
(w/ ND280 constraint) | 4.2% | 2.9% | 5.0%
FSI+SI+PN at SK | 2.5% | 1.5% | 10.5%
SK detector | 2.4% | 3.9% | 9.3%
All |  |  |
(w/o ND280 constraint) | 12.7% | 12.0% | 21.9%
(w/ ND280 constraint) | 5.5% | 5.1% | 14.8%
Table 20: Effect of 1$\sigma$ variation of the systematic uncertainties on the predicted event rates of the $\overline{\nu}$\- mode samples. Source of uncertainty |  $\overline{\nu}_{e}$ CCQE-like | $\overline{\nu}_{\mu}$
---|---|---
| $\delta N/N$ | $\delta N/N$
Flux | 3.8% | 3.8%
(w/ ND280 constraint) |  |
Cross section | 5.5% | 4.2%
(w/ ND280 constraint) |  |
Flux+cross-section |  |
(w/o ND280 constraint) | 12.9% | 11.3%
(w/ ND280 constraint) | 4.7% | 3.5%
FSI+SI+PN at SK | 3.0% | 2.1%
SK detector | 2.5% | 3.4%
All |  |
(w/o ND280 constraint) | 14.5% | 12.5%
(w/ ND280 constraint) | 6.5% | 5.3%

##  IX Sensitivity of oscillation parameters to neutrino interaction modeling

The neutrino interaction uncertainties are one of the main contributors to the
systematic uncertainty on all oscillation measurements and there is a global
effort underway to improve the understanding of neutrino cross sections. This
has lead to the creation of a number of interaction models which can, at least
partially, describe existing cross-section data but which also produce
different predictions for the oscillated event rates and spectra at SK. This
analysis uses the high statistics near detector data to constrain both the
flux and cross-section model uncertainties, improving the prediction of the
far detector event rate, and reducing the uncertainty on that prediction.
However, the near and far detectors observe different neutrino energy spectra
(mainly due to oscillations) so the underlying neutrino interactions they
sample will be different. The different design of the two detectors also means
that they are sensitive to different regions of the lepton kinematic phase
space. The near detector fit therefore tunes the flux and cross-section models
to a set of neutrino interactions and a lepton kinematic phase space that is
not the same as that observed at the far detector. These differences could
then be incorrectly attributed to neutrino oscillation effects in the
oscillation analysis, which would result in biased oscillation measurements.

This has been studied in a phenomenological context using a simplified
oscillation analysis, summarized in Ref. Ankowski:2016jdd . The results of the
studies show that long-baseline experiments may be biased by cross-section
model choices, as might be expected from the qualitative arguments above.
Given the potential impact of these model choices it is important to
investigate these issues using the full T2K oscillation analysis framework.
This has been performed for a range of neutrino interaction model variations,
which were discussed in Sec. III.2.

####  IX.0.1 Production and analysis of simulated datasets

The analysis is performed using replica datasets created by changing one part
of the nominal MC. Simulated data are created by applying the event selection
described in Sec. IV.2 and Sec. VI to the nominal T2K MC. A weight is then
applied to each selected event, calculated as the ratio of the altered
interaction model to the nominal cross section model. For the SK simulated
data, the relevant oscillation probability is also applied. This produces
event samples corresponding to the alternative interaction model. An example
of this is shown in Fig. 30 for the ND280 samples, where the left plot shows
the ratio of the alternative 2p2h simulated data to the nominal MC for the
FGD1 CC-0$\pi$ sample while the right shows the ratio for the alternative 1p1h
model.

Figure 30: Ratio of the alternative 2p2h simulated data (left) and alternative
1p1h simulated data (right) to the nominal MC, shown as a function of
reconstructed lepton kinematics for the FGD1 CC-0$\pi$ sample.

For each interaction model variation, simulated data are generated at both the
ND280 and SK. The nominal flux and interaction models are then fit to the
ND280 simulated data using the procedure described in Sec. V. The result is
then used as the ND280 input to the oscillation analysis, described in Sec.
VIII.1. The SK simulated data are fit to produce a set of likelihood contours
for the oscillation parameters. These likelihood contours are then compared to
the expected likelihood contour for fits where the nominal MC is used as data.
For each oscillation parameter a bias is calculated, defined as the difference
in the parameter best-fit point between the simulated data fit and the nominal
data fit, divided by the $1\sigma$ uncertainty on that parameter from the
nominal fit.

####  IX.0.2 Results

Results from the ND280 fit are used for central values and uncertainties for
the neutrino flux and cross-section model parameters. When fitting the nominal
MC, the parameter estimators are found to be unbiased while their
uncertainties are reduced. For the fits to simulated data generated from
alternative models the values of both the flux and cross-section parameters
change, as shown in Fig. 31. These show the results obtained by fitting the
data generated using the alternative 2p2h model, rather than the nominal
model.

Figure 31: The nominal values (red) and postfit values (blue) of the ND280
$\nu$-mode $\nu_{\mu}$ flux (left) and the cross-section (right) parameters,
shown as a fraction of their nominal value, following the near detector fit to
the alternative 2p2h ND280 simulated data.

Figure 31 demonstrates how differences between the true and model interaction
cross sections can be absorbed by both the flux and cross-section
parameters.The ND280 fit result is used to predict the unoscillated event
spectra at SK, producing the $\nu$-mode $\nu_{\mu}$ and $\nu_{e}$ event
spectra shown in Fig. 32. The predicted spectra falls between the nominal and
2p2h models and the uncertainty typically covers much of the difference
between the models.

Figure 32: Single-ring expected event spectrum at SK as a function of $E_{\nu}^{rec}$, assuming CCQE kinematics, for the selected $\mu$-like (left) and e-like (right) events in $\nu$-mode. The nominal SK prediction is shown by the dashed blue line, the SK simulated data by the solid black line and the prediction from the ND280 simulated data fit by the hashed red region. Table 21: Oscillation parameter values used as inputs for the studies of simulated data. Parameter | Maximal | Non-maximal
---|---|---
| oscillation | oscillation
$\sin^{2}(\theta_{23})$ | 0.523 | 0.6
$\delta_{CP}$ (rad) | -1.601 | -1.601
$\sin^{2}(\theta_{13})$ | 0.0217
$\sin^{2}(\theta_{12})$ | 0.304
$\Delta m^{2}_{21}$ $(10^{-5}$ eV${}^{2})$ | 7.53
$\Delta m^{2}_{32}$ $(10^{-3}$ eV${}^{2})$ | 2.51

To evaluate possible bias in the oscillation parameter estimators, a fit is
performed with SK simulated data. The likelihood surfaces for the oscillation
parameters of interest obtained by this fit are shown in Fig. 33. These
results use the true parameter values listed in the “maximal oscillation”
column of Tab. 21.

Figure 33: The likelihood surfaces (red) from the oscillation fit of the
alternative 2p2h simulated data at SK, using the result of the fit to the
ND280 simulated data as an input. The oscillation parameters of Tab. 21
(maximal mixing) were used. The likelihood surfaces from the nominal fit are
shown by the black line. The contours for $\sin^{2}\theta_{23}$ (top left),
$\Delta m^{2}_{32}$ (top right), $\delta_{CP}$ (bottom left) and $\delta_{CP}$
vs $\sin^{2}\theta_{13}$ (bottom right) are shown. For the one-dimensional
likelihoods the point of minimum $\Delta\chi^{2}$ is given in each legend.

Fig. 33 shows that if the alternative 2p2h model is correct, using the Nieves
model in the nominal MC produces small ($\sim$2%) biases in the measured
values of $\sin^{2}\theta_{23}$ and $\Delta m^{2}_{32}$ ($<1\%$). The study
also shows a change in the $\delta_{CP}$ likelihood contour, increasing the
exclusion around $\delta_{CP}=\pi/2$ by 1.2 units of $\Delta\chi^{2}$.
Overall, as seen in the bottom right of Fig. 33, the effect on the appearance
analysis is small, compared to the statistical uncertainty.

Figure 34 shows fits to the same alternative 2p2h simulated data as in Fig.
33, but using the non-maximal oscillation parameters from Tab. 21.

Figure 34: The likelihood surfaces (red) from the oscillation fit of the
alternative 2p2h simulated data at SK, using the result of the fit to the
ND280 simulated data as an input. The oscillation parameters of Tab. 21 (non-
maximal mixing) were used. The likelihood surfaces from the nominal fit are
shown by the black line. The contours for $\sin^{2}\theta_{23}$ (top left),
$\Delta m^{2}_{32}$ (top right), $\delta_{CP}$ (bottom left) and $\delta_{CP}$
vs $\sin^{2}\theta_{13}$ (bottom right) are shown. For the one-dimensional
likelihoods the point of minimum $\Delta\chi^{2}$ is given in each legend.

For non-maximal disappearance, a larger change is observed in the
$\sin^{2}\theta_{23}$ likelihood contour. The effect on the point estimate and
68% CL interval is relatively small.

Fig. 35 shows the results of the fits to alternative 1p1h simulated dataset
described in Sec. III.2. The alternative 1p1h model changes the T2K MC as a
function of both neutrino energy and the angle at which the lepton is produced
relative to the neutrino direction. Similarly to the alternative 2p2h case,
this model variation introduces a small bias in the measurement of $\Delta
m^{2}_{32}$ and has a small effect on the $\delta_{CP}$ likelihood contour.

Figure 35: The likelihood surfaces (red) from the oscillation fit of the
alternative 1p1h simulated data at SK, using the result of the fit to the
ND280 simulated data as an input. The oscillation parameters of Tab. 21
(maximal mixing) were used. The likelihood surfaces from the nominal fit are
shown by the black line. The contours for $\sin^{2}\theta_{23}$ (top left),
$\Delta m^{2}_{32}$ (top right), $\delta_{CP}$ (bottom left) and $\delta_{CP}$
vs $\sin^{2}\theta_{13}$ (bottom right) are shown. For the one-dimensional
likelihoods the point of minimum $\Delta\chi^{2}$ is given in each legend.

####  IX.0.3 Summary

Detailed studies have been performed to evaluate the sensitivity of the T2K
oscillation analysis to neutrino interaction modeling. By using the near
detector data to estimate flux and interaction model parameters, the
oscillation parameter estimates are found to be essentially unaffected by
changes to the interaction model. The largest observed effect on the
oscillation parameter likelihood contours is shown in Fig. 34, while Fig. 35
is more representative of the majority of the model variations discussed in
Sec. III.2.

A summary of the maximum bias observed for all of the alternative models
studied is shown in Tab. 22. The bias is presented as a fraction of the
expected 1$\sigma$ uncertainty on each oscillation parameter, and is the
maximum bias seen from all true oscillation parameter values tested. T2K does
not expect a significant constraint on $\delta_{CP}$ given the integrated POT
available for this analysis. As a result, it is difficult to quantify the
effect on $\delta_{CP}$ in a single number, and so this parameter is not
included in the summary table. For reference, the largest effect observed for
$\delta_{CP}$ is shown in Fig. 34.

Table 22: Summary of the maximum bias seen in oscillation parameters when fitting simulated datasets (defined in Sect. III.2), presented as a fraction of the expected 1$\sigma$ uncertainty on each parameter. Fits were performed for a number of true oscillation parameter assumptions. The numbers shown here are the maximum bias found for a given parameter and alternative model across all true oscillation parameter values tested. Alternative model | Maximum bias on parameter ($\sigma$)
---|---
|  $\Delta m^{2}_{32}$ |  $\sin^{2}\theta_{23}$ |  $\sin^{2}\theta_{13}$
SF |  0.09 |  0.17 |  0.17
Effective RPA |  0.00 |  0.06 |  0.00
Alternative 1p1h |  0.20 |  0.09 |  0.07
Alternative 2p2h |  0.20 |  0.21 |  0.18
Delta-enhanced 2p2h |  0.10 |  0.10 |  0.00
Not-Delta 2p2h |  0.11 |  0.07 |  0.05

##  X $\overline{\nu}_{e}$ appearance analysis

T2K has already observed $\nu_{e}$-appearance in a $\nu_{\mu}$ beam
Abe:2013hdq , but no direct evidence of $\overline{\nu}_{e}$-appearance has
been reported so far. A search for this process, with the data collected by
T2K in $\overline{\nu}$-mode, has been performed using the analysis tools
described in Sec. VIII.

In order to look for $\overline{\nu}_{e}$-appearance independently from
$\nu_{e}$-appearance, a parameter $\beta$ is introduced which multiplies the
$\overline{\nu}_{e}$-appearance probability

$P(\overline{\nu}_{\mu}\rightarrow\overline{\nu}_{e})\rightarrow\beta\times
P(\overline{\nu}_{\mu}\rightarrow\overline{\nu}_{e})$ (23)

In this analysis, $\beta$ can have two values: $\beta=0$, which corresponds to
no $\overline{\nu}_{\mu}\rightarrow\overline{\nu}_{e}$ oscillations; and
$\beta=1$, which corresponds to the appearance probability as predicted by the
PMNS framework.

The two hypotheses are tested using rate-only and rate-plus-shape analyses.
The test statistic is different in the two cases. The test statistic used in
the rate-only analysis is the number of e-like candidates observed at SK in
$\overline{\nu}$-mode, while for the rate-plus-shape analysis, the likelihood
ratio for the two hypotheses, assumining PMNS oscillations and no
$\overline{\nu}_{\mu}\rightarrow\overline{\nu}_{e}$ is used:

$\displaystyle\Delta\chi^{2}$
$\displaystyle=\chi^{2}(\beta=0)-\chi^{2}(\beta=1)$
$\displaystyle=2\Big{[}-\;\textrm{ln}\;\mathcal{L}(\beta=0;x)+\;\textrm{ln}\;\mathcal{L}(\beta=1;x)\Big{]},$
(24)

A toy ensemble is generated according to the prior knowledge on the
oscillation parameters defined in Sec. VII.2 for both values of $\beta$. The
T2K information from the $\nu_{\mu}$-disappearance,
$\overline{\nu}_{\mu}$-disappearance and $\nu_{e}$-appearance 111Note that the
e-like CC1$\pi^{+}$ sample was not included in this analysis, only CCQE-like
samples were used. channels are taken into account by using a posterior
predictive method demortier2008p . In this method, the data from these
channels are used to constrain the prior parameter space in generating the toy
ensemble for the $\overline{\nu}_{e}$ sample while preserving their
correlations.

The sensitivity of the analysis is computed by producing a simulated sample
without statistical fluctuations using the values of the oscillation
parameters defined in Tab. 13. The test statistics for this particular dataset
are 6.28 expected $\overline{\nu}_{e}$ candidates for the rate-only analysis,
corresponding to the total number of expected events in $\overline{\nu}$-mode
in Tab. 14, and $\Delta\chi^{2}=2.54$ for the rate-plus-shape analysis. When
compared to the generated toy ensembles for the no-$\overline{\nu}_{e}$
appearance hypothesis the p-values are 0.047 for both, the rate-only and rate-
plus-shape analyses. When compared with the toy ensemble generated with
$\beta=1$, p-values of 0.52 and 0.41 are found for the rate-only and rate-
plus-shape analyses respectively.

The same analysis is then applied to the $e$-like data selected in
$\overline{\nu}$-mode. The test statistics in data for the rate-only and rate-
plus-shape analyses are $N^{obs}=4$ $\overline{\nu}_{e}$ candidates and
$\Delta\chi^{2}=-2.51$, obtained from the {$E_{\nu}^{rec}$, $\theta_{lep}$}
distribution of the $e$-like candidates in $\overline{\nu}$-mode shown in Fig.
28. As shown in Tab. 23 this test statistics return a p-value of 0.41 (0.46)
for the $\beta=0$ hypothesis for the rate-only (rate-plus-shape) analysis,
while the p-value for $\beta=1$ is 0.21 (0.07). This analysis shows that, with
the available data, the rate-plus-shape analysis excludes $\beta=1$ at 90% CL.

The observed p-values for the $\beta=0$ hypothesis are larger than expectation
due to the lower than expected number of observed events in the e-like
$\overline{\nu}$-mode sample. The lower p-value for the rate $\beta=1$
hypothesis is driven by the larger discrepancy between the selected $\nu_{e}$
and $\overline{\nu}_{e}$ rate (32 and 4 events respectively) than predicted by
the simulated dataset generated, even under the assumption of a maximal CP-
violation $\delta=-\pi/2$ (28.55 and 6.28 events respectively). The p-value in
the rate-plus-shape analysis for the $\beta=1$ hypothesis is reduced further,
due to the distribution of the $\overline{\nu}_{e}$-appearance candidates in
electron momentum and angle. Their distribution is more compatible with the
background expected distribution, as shown in Fig. 19.

Table 23: Expected and observed (data) p-values for the $\beta=0$ and $\beta=1$ hypotheses, for both rate-only and rate-plus-shape analyses. The expected p-values are estimated for the simulated dataset defined in Tab. 13. | Rate-only | Rate-plus-shape
---|---|---
Simulated sample |  |
$\beta=0$ | 0.047 | 0.047
$\beta=1$ | 0.52 | 0.41
Data |  |
$\beta=0$ | 0.41 | 0.46
$\beta=1$ | 0.21 | 0.07

##  XI Joint neutrino and antineutrino oscillation analysis results

In this section, joint three-flavor oscillation analyses performed with both
the frequentist and the Bayesian approaches are presented. The five SK samples
introduced in Sec. VI are used, which allows the simultaneous study of the
$\nu_{e}$ and $\overline{\nu}_{e}$ appearance channels and the $\nu_{\mu}$ and
$\overline{\nu}_{\mu}$ disappearance channels. The oscillation parameters
$|\Delta m^{2}_{32}|$, $\textrm{sin}^{2}\theta_{23}$,
$\textrm{sin}^{2}\theta_{13}$, $\delta_{\mathrm{CP}}$, and the mass ordering
are determined with and without using the measurement of
$\textrm{sin}^{2}\theta_{13}$ from reactor experiments as a constraint.

###  XI.1 Frequentist analysis

Although two frequentist analyses were introduced in Sec. VIII.1, the results
are similar, so detailed results are only presented for the analysis which
uses {$E_{\nu}^{rec}$, $\theta_{lep}$} templates in this section. Comparisons
between both frequentist, and the Bayesian analysis are shown in Sec. XI.3.

####  XI.1.1 Results without reactor constraints

This section describes the results obtained by the frequentist analysis when
only T2K data are used to estimate the oscillation parameters. The main
parameters of interest in this case are $\delta_{\mathrm{CP}}$ and
$\sin^{2}\theta_{13}$ that can be directly compared to the reactor
measurements. The point estimates for these oscillation parameters and the
constant $\Delta\chi^{2}=1$ intervals are given for normal and inverted
ordering in Tab. 24. The $\Delta\chi^{2}$ surfaces are shown in Fig. 36. These
intervals have been produced via marginalization of all nuisance and
oscillation parameters which are not of interest, as described in Sec. VIII.1.

Table 24: Point estimates and 1$\sigma$ confidence intervals under the constant $\Delta\chi^{2}$ approximation from an analysis considering T2K oscillation data only. Parameter | Normal ordering | Inverted ordering
---|---|---
Best-fit | $\pm 1\sigma$ | Best-fit | $\pm 1\sigma$
$\delta_{CP}$ | -1.791 | [-2.789; -0.764] | -1.382 | [-2.296;-0.524]
$\sin^{2}\theta_{13}$ | 0.0271 | [0.0209; 0.0342] | 0.0299 | [0.0232; 0.0380]

Figure 36:  One-dimensional $\Delta\chi^{2}$ surfaces for oscillation
parameters $\delta_{CP}$ and $\sin^{2}\theta_{13}$ using T2K-only data. The
yellow band on the right plot corresponds to the reactor value on
$\sin^{2}\theta_{13}$ from the PDG 2015 PDG2015 .

Two-dimensional contours of constant $\Delta\chi^{2}$ for the parameters
$\delta_{CP}$ and $\sin^{2}\theta_{13}$, along with a comparison with the
constraint on $\sin^{2}\theta_{13}$ from reactor experiments, are shown in
Fig. 37. The point estimate and constant $\Delta\chi^{2}$ confidence intervals
of $\sin^{2}\theta_{13}$ from T2K data only are slightly larger than what is
found by the reactor experiments. However, the T2K-only measurement of
$\sin^{2}\theta_{13}$ is still consistent with the reactor measurement at the
68% confidence level.

Figure 37:  Two-dimensional constant $\Delta\chi^{2}$ contours for oscillation
parameters $\delta_{CP}$ and $\sin^{2}\theta_{13}$ using T2K data only. The
yellow band corresponds to the reactor value on $\sin^{2}\theta_{13}$ from the
PDG 2015 PDG2015 .

As mentioned above, in this analysis a fifth sample selecting $\nu_{e}$
candidates at SK with one delayed Michel electron in the final state has been
added for the first time. A comparison of the $\Delta\chi^{2}$ surface for
$\delta_{\mathrm{CP}}$ only including the four single-ring samples used in
previous analyses and the results obtained with the inclusion of the fifth
sample is shown in Fig. 38. As expected, a small improvement is observed when
the new sample is included.

Figure 38:  A comparison of one-dimensional constant $\Delta\chi^{2}$ contours
for normal ordering for $\delta_{CP}$ using T2K-only data for the four- and
five-sample fits.

Fig. 39 shows a comparison of the constraints in the
$\delta_{CP}$–$\sin^{2}\theta_{13}$ plane when appearance channels taken in
$\nu$-mode and in $\overline{\nu}$-mode are considered independently. Both
$\nu$\- and $\overline{\nu}$-mode disappearance channels are used in both
fits. The $\nu$\- and $\overline{\nu}$-mode datasets alone prefer different
values of $\sin^{2}\theta_{13}$, which is driven by the absolute appearance
rate. It is clear that the $\overline{\nu}$-mode appearance sample does not
have the power to exclude a zero value of $\sin^{2}\theta_{13}$ by itself. In
either case, the reactor measurement of $\sin^{2}\theta_{13}$ sits near the
upper and lower 68% confidence contours for the $\overline{\nu}$-mode and
$\nu$-mode samples respectively.

Figure 39:  Contours in the $\sin^{2}\theta_{13}$–$\delta_{CP}$ plane using
T2K-only data, obtained by analysing either the $\nu$\- or
$\overline{\nu}$-mode appearance datasets are compared for both orderings.
Both $\nu$\- and $\bar{\nu}$-mode disappearance datasets were used in all
fits. The yellow band corresponds to the reactor value on
$\sin^{2}\theta_{13}$ from the PDG 2015 PDG2015 .

####  XI.1.2 Results with reactor constraints

Here, the oscillation parameters obtained by the T2K data fit where
$\sin^{2}\theta_{13}$ is marginalized using the reactor constraint given in
Tab. 18. The best fit values of the T2K data with the reactor constraint are
summarized in Tab. 25.

Fig. 40 shows the 90% constant $\Delta\chi^{2}$ surface in the
$\sin^{2}\theta_{23}$–$\Delta m_{32}^{2}$ plane, assuming normal mass
ordering. The interval is compared with other experiments, showing good
agreement with IceCube and SK and some tension with MINOS and NO$\nu$A.

Figure 40: Allowed region at 90% confidence level for oscillation parameters
$\sin^{2}\theta_{23}$ and $\Delta m_{32}^{2}$ using T2K data with the reactor
constraint ($\sin^{2}(2\theta_{13})=0.085\pm 0.005$). The normal mass ordering
is assumed and the T2K results are compared with NO$\nu$A Adamson:2017qqn ,
MINOS PhysRevLett.110.251801 , Super-K Wendell:2015onk , and IceCube
Aartsen:2016psd .

The NO$\nu$A collaboration published $\nu_{\mu}$-disappearance results
disfavoring maximal mixing for $\textrm{sin}^{2}\theta_{23}$ at 2.6$\sigma$
Adamson:2017qqn . The T2K data in the $\nu_{\mu}$\- and
$\overline{\nu}_{\mu}$-disappearance channels, together with the T2K best fit
and the expected spectrum produced using the NO$\nu$A best fit value for
$\textrm{sin}^{2}\theta_{23}$ (higher octant) and $\Delta m^{2}_{32}$, are
shown in Fig. 41.

Figure 41: Comparison of the T2K data in $\nu_{\mu}$ (left) and
$\overline{\nu}_{\mu}$ (right) disappearance channels with the expected
spectra obtained with the T2K most probable values of the oscillation
parameters and using the NO$\nu$A most probable values for
$\textrm{sin}^{2}\theta_{23}$ (higher octant) and $\Delta m^{2}_{32}$ taken
from Ref. Adamson:2017qqn .

Fig. 42 shows the 2D $\sin^{2}\theta_{13}$–$\delta_{CP}$ confidence level
contours for the data fit including the reactor constraint. The comparisons
with the four-sample joint fit are also shown to demonstrate the effect of the
inclusion of the CC1$\pi^{+}$ e-like sample in the appearance analysis.
Compared to the best-fit results obtained with the T2K-only data fit in Sec.
XI.1.1, the inclusion of the CC1$\pi^{+}$ e-like sample in the data fit with
the reactor constraint results in a shift of best-fit value for the
$\delta_{CP}$ phase towards the maximally violating phase of $-\pi/2$.

Table 25: Best-fit results and the 1$\sigma$ confidence interval of the T2K data fit with the reactor constraint with normal and inverted hypotheses. Parameter | Normal ordering | Inverted ordering
---|---|---
Best-fit | $\pm 1\sigma$ | Best-fit | $\pm 1\sigma$
$\delta_{CP}$ | -1.728 | [-2.538;-0.877] | -1.445 | [-2.170;-0.768]
$\sin^{2}\theta_{23}$ | 0.550 | [0.465;0.601] | 0.5525 | [0.470;0.601]
$\Delta m^{2}_{32}$ | 2.54 | [2.460;2.621] | 2.51 | [2.429;2.588]
($10^{-3}$ eV${}^{2}/c^{4}$)

Figure 42:  A comparison of two-dimensional constant $\Delta\chi^{2}$ contours
in the $\delta_{CP}$–$\sin^{2}\theta_{13}$ plane using T2K data with the
reactor constraint, for both four-sample (red) and five-sample (black)
analyses with normal (left) and inverted (right) mass ordering hypotheses. The
contours are produced by marginalizing the likelihood with respect to all
parameters other than the parameters of interest.

Since there is a physical boundary at $\delta_{CP}=\pm\frac{\pi}{2}$,
calculating the coverage near the boundary using a Gaussian approximation may
not produce accurate results. To solve this problem, the coverage of the 1D
$\Delta\chi^{2}$ distribution as a function of $\delta_{CP}$ is computed using
the Feldman-Cousins approach, described in Sec. VIII. In order to perform the
study, 10,000 toy MC experiments were generated for different values of
$\delta_{CP}$ and the mass ordering. The 1D $\Delta\chi^{2}$ surface obtained
with the Feldman-Cousins approach is used to evaluate the 90% confidence
intervals for $\delta_{CP}$ in both ordering cases, as shown in Fig. 43. In
this analysis, CP-conserving values of $\delta_{\mathrm{CP}}~{}=~{}0,~{}\pi$
are excluded at 90% and 2$\sigma$ confidence levels respectively. Values of
$\delta_{\mathrm{CP}}$ in the intervals [-2.95,-0.44] ([-1.47, -1.27]) are
allowed at 90% confidence for normal (inverted) ordering.

Figure 43: One dimensional $\Delta\chi^{2}$ surfaces for oscillation parameter
$\delta_{CP}$ using T2K data with the reactor constraint. The critical
$\Delta\chi^{2}$ values obtained with the Feldman-Cousins method are used to
evaluate the 90% confidence level with the proper coverage.

A useful way to visualize the results is to compare the observed number of
events in the $\nu$-mode (in both CCQE-like and CC1$\pi^{+}$-like samples) and
$\overline{\nu}$-mode e-like samples with the expected events for different
values of $\delta_{\mathrm{CP}}$, $\sin^{2}\theta_{23}$, and mass ordering. As
it is shown in Fig. 44 the T2K data falls outside the physically allowed
region.

Figure 44: Total predicted $\nu^{\hbox{\kern-8.0pt\raise
1.0pt\hbox{{\tiny(}{\lower 1.4pt\hbox{\bf--}}{\tiny)}}}}_{e}$-appearance event
rates in the $\nu$-mode samples and in the $\bar{\nu}$-mode sample as a
function of $\delta_{CP}$ for different values of $\sin^{2}\theta_{23}$ and
both mass orderings, compared to T2K data. The dashed line distinguishes the
two solutions for the octant of $\theta_{23}$.

In order to quantify whether the T2K dataset is consistent with the PMNS
framework in terms of significance, an additional toy MC study was performed.
An ensemble of 10,000 simulated datasets was obtained in the same way as
described in Sec. VIII for the Feldman-Cousins method, with
$\delta_{CP}=-\pi/2$ and normal mass ordering. The values of $-2\Delta\ln L$
that contain 68.3% and 95.5% of the MC toys were computed and compared to the
distribution obtained with the fit of the T2K dataset. As shown in Fig. 45,
the T2K data $-2\Delta\ln L$ distribution shows a less extreme fluctuation
than at least 5% of the toys MC for all the values of $\delta_{CP}$ and mass
ordering, i.e. if the experiment is repeated many times and the true value is
$\delta_{CP}=-\pi/2$ with normal ordering, more than 5% of the experiments are
expected to show a more extreme statistical fluctuation than the current T2K
dataset over the whole range of $\delta_{CP}$ and mass ordering. From Fig. 45,
the fraction of experiments that would exclude $\delta_{CP}=0,\pi$ at 90% or
$2\sigma$ confidence level can be estimated. Assuming a true value of
$\delta_{\mathrm{CP}}$ of -$\pi$/2 and normal ordering, 24.3% (21.3%) of toy
MC experiments exclude $\delta_{\mathrm{CP}}=0$ ($\pi$) at 90% CL. The same
can be repeated for different values of $\delta_{CP}$ and mass ordering as
shown in Tab. 26.

Figure 45: One-dimensional marginal $\Delta\chi^{2}$ surfaces for oscillation
parameters $\delta_{CP}$ and $\sin^{2}\theta_{13}$ using T2K data with the
reactor constraint. The contour is produced by marginalizing the likelihood
with respect to all parameters other than the parameter of interest. The red
line shows the critical $\Delta\chi^{2}$ values obtained with the Feldman-
Cousins method, used to evaluate the 90% confidence level with the proper
coverage. The green line show the $\Delta\chi^{2}$ obtained with the fit to
the T2K data. Table 26:  The fraction of toy experiments for which
$\delta_{CP}=0,\pi$ and normal and inverted ordering are excluded at 90% and
$2\sigma$ confidence is shown for different true values of $\delta_{CP}$ and
mass ordering. 10,000 toy experiments are used for each set of values.  True:
$\delta_{CP}=-\pi/2$ — normal ordering
---
$\delta_{CP}$ | Ordering | 90% CL |  $2\sigma$ CL
0 | Normal | 0.243 | 0.131
$\pi$ | Normal | 0.216 | 0.105
0 | Inverted | 0.542 | 0.425
$\pi$ | Inverted | 0.559 | 0.436
True: $\delta_{CP}=0$ — normal ordering
$\delta_{CP}$ | Ordering | 90% CL |  $2\sigma$ CL
0 | Normal | 0.104 | 0.0490
$\pi$ | Normal | 0.130 | 0.0591
0 | Inverted | 0.229 | 0.137
$\pi$ | Inverted | 0.205 | 0.122
True: $\delta_{CP}=-\pi/2$ — inverted ordering
$\delta_{CP}$ | Ordering | 90% CL |  $2\sigma$ CL
0 | Normal | 0.124 | 0.0515
$\pi$ | Normal | 0.102 | 0.0413
0 | Inverted | 0.290 | 0.194
$\pi$ | Inverted | 0.308 | 0.207

###  XI.2 Bayesian analysis

####  XI.2.1 Results without reactor constraints

This section describes the results obtained by the Bayesian analysis when
using only T2K data to estimate the parameters $\sin^{2}\theta_{23}$, $\Delta
m_{32}^{2}$, $\sin^{2}\theta_{13}$ and $\delta_{CP}$ with the MCMC method
described in Sec. VIII.2. In contrast with the frequentist analysis presented
in Sec. XI.1, the Markov chain walks in a parameter space where the sign of
$\Delta m_{32}^{2}$ can flip, and results are presented for both mass
orderings. The best-fit point and $\pm 1\sigma$ credible interval for each
parameter, obtained with the KDE method, are summarized in Tab. 27. The best
fit point is the mode of the four-dimensional histogram where the axes are the
oscillation parameters.

Table 27: Best-fit results and the 1$\sigma$ credible interval of the T2K data fit without the reactor constraint with the MCMC analyses including both mass orderings. Parameter | Best-fit | $\pm 1\sigma$
---|---|---
$\delta_{CP}$ | -1.815 | [-2.275; -0.628]
$\sin^{2}\theta_{13}$ | 0.0254 | [0.0210; 0.0350]
$\sin^{2}\theta_{23}$ | 0.513 | [0.460 ; 0.550]
$\Delta m_{32}^{2}$ | $2.539\times 10^{-3}{eV^{2}/c^{4}}$ | $[-2.628;-2.544]\times 10^{-3}{eV^{2}/c^{4}}$
$[2.436;2.652]\times 10^{-3}{eV^{2}/c^{4}}$

The $\pm 1\sigma$ credible intervals, which have a 68.3% probability of
containing the true value, are computed, for each parameter, from the
posterior probability density marginalized over all the other parameters as
shown in Fig. 46. Fig. 46 also shows the correlations between the oscillation
parameters with the map of the marginal posterior density probability and the
credible intervals in the space formed by two parameters.

Figure 46: The two-dimensional histograms represent the marginal posterior
probability in the two-parameter space as a blue gradient. The white solid
(dashed) line is the 90% ($1\sigma$) credible interval. The one-dimensional
histograms represent the posterior probability density (post. proba.) of the
oscillation parameter in the x-axis of the column marginalized over all other
parameters. The blue areas are respectively the $1\sigma$ (dark), 90%
(medium), and 95% (light) credible interval.

The proportion of the MCMC points with $\sin^{2}\theta_{23}>0.5$ or $<0.5$
gives the posterior probability of the octant. Similarly, the relative
proportion of steps with $\Delta m_{32}^{2}>$ or $<0$ gives the posterior
probability of each mass ordering. They are shown in Tab. 28. A Bayes factor
can be computed as a ratio of the posterior probabilities jeffreys1998 . The
Bayes factor for normal ordering is $B(\mathrm{NH}/\mathrm{IH})=2.28$; the
Bayes factor for the upper octant is
$B(\sin^{2}\theta_{23}>0.5/\sin^{2}\theta_{23}<0.5)=1.32$. Neither can be
considered decisive.

Table 28: Posterior probabilities for the mass orderings and $\sin^{2}\theta_{23}$ when fitting T2K data only with an MCMC method. | $\sin^{2}\theta_{23}<0.5$ | $\sin^{2}\theta_{23}>0.5$ | Line Total
---|---|---|---
Inverted ordering | 0.137 | 0.168 | 0.305
Normal ordering | 0.294 | 0.401 | 0.695
Column total | 0.431 | 0.569 | 1

####  XI.2.2 Results with reactor constraints

This section presents the results obtained with the MCMC analysis when adding
a Gaussian prior on $\sin^{2}\theta_{13}$ with the value given in Tab. 18. The
posterior mode marginalized over the nuisance parameters is given in Tab. 29.
Including the reactor prior on $\sin^{2}\theta_{13}$, the best-fit is closer
to that obtained by the reactor experiments compared to the T2K-only results.
The $\delta_{CP}$ best-fit is closer to the maximum violating value of
$-\pi/2$ due to the correlations with $\sin^{2}\theta_{13}$ shown in Fig. 46.

Table 29: Best-fit results and the 1$\sigma$ credible interval of the T2K data fit with the reactor constraint with the MCMC analyses including both mass orderings. Parameter | Best-fit | $\pm 1\sigma$
---|---|---
$\delta_{CP}$ | -1.789 | [-2.450; -0.880]
$\sin^{2}\theta_{13}$ | 0.0219 | [0.0208; 0.0233]
$\sin^{2}\theta_{23}$ | 0.534 | [0.490 ; 0.580]
$\Delta m_{32}^{2}$ | $2.539\times 10^{-3}\,\mathrm{eV^{2}/c^{4}}$ | [-3.000; -2.952] $\times 10^{-3}{eV^{2}/c^{4}}$
[2.424; 2.664]$\times 10^{-3}{eV^{2}/c^{4}}$

The MCMC algorithm uses a flat prior on $\delta_{CP}$, but its dependence on
this choice of prior has been tested by computing the credible intervals with
a flat prior on $\sin\delta_{CP}$. The two sets of intervals are in reasonable
agreement as shown in Fig. 47.

Figure 47: $\delta_{\mathrm{CP}}$ marginal posterior probability as obtained
with the MCMC method. The credible intervals for $\pm 1\sigma$, $\pm 90\%$ and
$\pm 95\%$ are shown when using a flat prior in $\delta_{CP}$ (usual fit), and
when converting to a flat prior in $\sin\delta_{\mathrm{CP}}$.

The Bayes factor for the mass ordering and the $\theta_{23}$ octant can be
computed with the method described in Sec. XI.2.1. Using the values from Tab.
30, they are found to be $B(\mathrm{NH}/\mathrm{IH})=3.71$ and
$B(\sin^{2}\theta_{23}>0.5/\sin^{2}\theta_{23}<0.5)=2.39$ respectively. Also
in this case, these cannot be considered decisive.

Table 30: Posterior probabilities for the mass orderings and $\sin^{2}\theta_{23}$ with an MCMC method when fitting T2K data with the reactor constraint. | $\sin^{2}\theta_{23}<0.5$ | $\sin^{2}\theta_{23}>0.5$ | Line Total
---|---|---|---
Inverted ordering | 0.060 | 0.152 | 0.212
Normal ordering | 0.235 | 0.553 | 0.788
Column total | 0.295 | 0.705 | 1

###  XI.3 Comparison among the oscillation analyses

The frequentist likelihood multiplied by Gaussian penalty terms for the
nuisance parameters and uniform priors for the oscillation parameters is
equivalent to the Bayesian posterior density, for the same priors. In order to
compare the analyses, the posterior probability densities sampled by the
Bayesian analyses are converted into a $\Delta\chi^{2}$ function and the
intervals are recalculated to extract confidence intervals that are compared
with the frequentist analyses. Fig. 48 shows the constant $\Delta\chi^{2}$ 68%
and 90% intervals for all three oscillation analyses in the
$\sin^{2}\theta_{13}$–$\delta_{CP}$ plane, assuming normal ordering and only
using T2K data. Differences exist among the three methods as the 2D templates
fitted in the appearance samples are different and the Bayesian analyses does
a combined fit of near and far detector samples but no major differences are
found between the contours.

Figure 48: Constant $\Delta\chi^{2}$ 68% and 90% intervals in the
$\sin^{2}\theta_{13}$–$\delta_{CP}$ plane for both frequentist analyses and
the Bayesian fit.

Fig. 49 shows the constant $\Delta\chi^{2}$ 68% and 90% intervals in the
$\sin^{2}\theta_{23}$–$\Delta m_{32}^{2}$ plane for both frequentist and
Bayesian fits. Both distributions and intervals agree between fitters,
validating the extrapolation of the constraints on the nuisance parameters
obtained in the near detector fit to SK.

Figure 49: Constant $\Delta\chi^{2}$ 68% and 90% intervals in the
$\sin^{2}\theta_{23}$–$\Delta m_{32}^{2}$ plane for the frequentist fit which
uses 2D {$E_{\nu}^{rec}$, $\theta_{lep}$} templates and the Bayesian fit,
assuming normal ordering.

###  XI.4 Best fit spectra

An estimate of the oscillation parameters $\Delta m_{32}^{2}$,
$\sin^{2}\theta_{13}$ and $\delta_{CP}$ have been obtained with both
frequentist and Bayesian analyses. The agreement between the fit results and
the data has been evaluated by comparing the expected spectra after the
oscillation fit with the data points as shown in Fig. 50. The best-fit spectra
are obtained by sampling 2000 points from the MCMC including the reactor
constraint and fitting a Gaussian distribution to calculate the most probable
value for the predicted number of events in each energy bin.

In order to extract the ratio of oscillated to unoscillated spectra, the
expected spectra are also tuned to the no oscillation case. A coarser binning
than the one used in the fit has been used for readability.

Figure 50: Comparison of the best-fit oscillated MC energy spectra,
unoscillated spectra and T2K data for the five samples used in the fit:
$\mu$-like sample in $\nu$-mode and $\overline{\nu}$-mode (top left and
right), single ring e-like appearance sample in $\nu$-mode and
$\overline{\nu}$-mode (middle left and right), CC1$\pi^{+}$ e-like appearance
sample in $\nu$-mode (bottom). The larger unoscillated spectra in the
CC1$\pi^{+}$ e-like sample compared to the single ring sample is due to the
relatively large background of $\nu_{\mu}$ in the CC1$\pi^{+}$ sample, which
does not disappear in the no-oscillation case. The ratio of the best fit to
unoscillated spectra are also shown.

##  XII Conclusions

All data collected by the T2K experiment between 2010 and 2016 have been
analyzed to estimate the oscillation parameters $|\Delta m^{2}_{32}|$,
$\textrm{sin}^{2}\theta_{23}$, $\textrm{sin}^{2}\theta_{13}$,
$\delta_{\mathrm{CP}}$ and the mass ordering. These parameters are estimated
by doing a joint fit of the $\nu_{\mu}$ and $\overline{\nu}_{\mu}$
disappearance channels and $\nu_{e}$ and $\overline{\nu}_{e}$ appearance
channels by using five samples selected at the far detector, thus, including
the new additional CC1$\pi^{+}$ sample. The data related to this measurement
can be found in datarelease . A comprehensive study has been performed to
evaluate the sensitivity of oscillation parameter estimates to neutrino
interaction modeling, showing that the impact of these uncertainties is small
compared to the total uncertainties on the measurement of all the oscillation
parameters.

The general approach followed in this paper, which combines separate analyses
of beamline, neutrino interactions, near and far detector selections through
sets of systematic parameters and their covariances will be extended with
additional data which will be collected by T2K in the coming years in both
$\nu$\- and $\overline{\nu}$-modes, and improved near and far detector
samples. This is expected to greatly enhance the sensitivity of the T2K
experiment to the measurement of the CP-violation phase $\delta_{\mathrm{CP}}$
as well as more precise measurements of the atmospheric parameters $|\Delta
m^{2}_{32}|$ and $\textrm{sin}^{2}\theta_{23}$.

###### 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, Japan; NSERC (Grant No. SAPPJ-2014-00031),
NRC and CFI, Canada; CEA and CNRS/IN2P3, France; DFG, Germany; INFN, Italy;
National Science Centre (NCN) and Ministry of Science and Higher Education,
Poland; RSF, RFBR, and MES, Russia; MINECO and ERDF funds, Spain; SNSF and
SERI, Switzerland; STFC, UK; and DOE, USA. We also thank CERN for the
UA1/NOMAD magnet, DESY for the HERA-B magnet mover system, NII for SINET4, the
WestGrid and SciNet consortia in Compute Canada, and GridPP and the Emerald
High Performance Computing facility in the United Kingdom. In addition,
participation of individual researchers and institutions has been further
supported by funds from ERC (FP7), H2020 Grant No. RISE-GA644294-JENNIFER, EU;
JSPS, Japan; Royal Society, UK; the Alfred P. Sloan Foundation and the DOE
Early Career program, USA.

## References

  * [1] Y. Fukuda et al.  Evidence for oscillation of atmospheric neutrinos.  Phys. Rev. Lett., 81:1562–1567, 1998.
  * [2] Q. R. Ahmad et al.  Direct evidence for neutrino flavor transformation from neutral current interactions in the Sudbury Neutrino Observatory.  Phys. Rev. Lett., 89:011301, 2002.
  * [3] M. H. Ahn et al.  Measurement of neutrino oscillation by the K2K experiment.  Phys. Rev., D74:072003, 2006.
  * [4] P. Adamson et al.  Combined analysis of $\nu_{\mu}$ disappearance and $\nu_{\mu}\rightarrow\nu_{e}$ appearance in MINOS using accelerator and atmospheric neutrinos.  Phys.Rev.Lett., 112:191801, 2014.
  * [5] K. Eguchi et al.  First results from KamLAND: Evidence for reactor anti-neutrino disappearance.  Phys. Rev. Lett., 90:021802, 2003.
  * [6] K. Abe et al.  Indication of Electron Neutrino Appearance from an Accelerator-produced Off-axis Muon Neutrino Beam.  Phys. Rev. Lett., 107:041801, 2011.
  * [7] K. Abe et al.  Observation of Electron Neutrino Appearance in a Muon Neutrino Beam.  Phys.Rev.Lett., 112:061802, 2014.
  * [8] F. P. An et al.  Observation of electron-antineutrino disappearance at Daya Bay.  Phys. Rev. Lett., 108:171803, 2012.
  * [9] Soo-Bong Kim et al.  Observation of Reactor Electron Antineutrino Disappearance in the RENO Experiment.  Phys. Rev. Lett., 108:191802, 2012.
  * [10] Y. Abe et al.  Indication for the disappearance of reactor electron antineutrinos in the Double Chooz experiment.  Phys. Rev. Lett., 108:131801, 2012.
  * [11] P. Adamson et al.  First measurement of electron neutrino appearance in NOvA.  Phys. Rev. Lett., 116(15):151806, 2016.
  * [12] K. Abe et al.  Combined Analysis of Neutrino and Antineutrino Oscillations at T2K.  Phys. Rev. Lett., 118(15):151801, 2017.
  * [13] T. Sekiguchi et al.  Development and operational experience of magnetic horn system for T2K experiment.  Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 789:57 – 80, 2015.
  * [14] K. Suzuki et al.  Measurement of the muon beam direction and muon flux for the t2k neutrino experiment.  Progress of Theoretical and Experimental Physics, 2015(5):053C01, 2015.
  * [15] K. Abe et al.  The T2K Experiment.  Nucl. Instrum. Methods Phys. Res., Sect. A, 659:106–135, 2011.
  * [16] K. Abe et al.  T2K neutrino flux prediction.  Phys. Rev. D, 87:012001, Jan 2013.
  * [17] A. Ferrari, P. R. Sala, A. Fasso, and J. Ranft.  FLUKA: A multi-particle transport code.  CERN-2005-010, SLAC-R-773, INFN-TC-05-11.
  * [18] T.T. Bohlen, F. Cerutti, M.P.W. Chin, A. Fasso, A. Ferrari, P.G. Ortega, A. Mairani, P.R. Sala, G. Smirnov, and V. Vlachoudis.  The FLUKA Code: Developments and Challenges for High Energy and Medical Applications.  Nuclear Data Sheets, 120:211–214, 2014.
  * [19] N Abgrall et al.  Measurements of Cross Sections and Charged Pion Spectra in Proton-Carbon Interactions at 31 GeV/c.  Phys. Rev. C, 84:034604, 2011.
  * [20] N. Abgrall et al.  Measurement of Production Properties of Positively Charged Kaons in Proton-Carbon Interactions at 31 GeV/c.  Phys. Rev. C, 85:035210, 2012.
  * [21] N. Abgrall et al.  Measurements of $\pi^{\pm}$, $K^{\pm}$, $K^{0}_{S}$, $\Lambda$ and proton production in proton–carbon interactions at 31 GeV/c with the NA61/SHINE spectrometer at the CERN SPS.  Eur. Phys. J. C, 76:84, 2016.
  * [22] T. Eichten et al.  Particle production in proton interactions in nuclei at 24 GeV/c.  Nucl. Phys. B, 44, 1972.
  * [23] J. V. Allaby et al.  Technical Report 70-12, CERN, 1970.
  * [24] I. Chemakin et al.  Pion production by protons on a thin beryllium target at 6.4-Ge, 12.3-GeV/c, and 17.5-GeV/c incident proton momenta.  Phys. Rev. C, 77:015209, 2008.
  * [25] R. Brun, F. Carminati, and S. Giani.  CERN-W5013, 1994.
  * [26] C. Zeitnitz and T. A. Gabriel.  Proc. of International Conference on Calorimetry in High Energy Physics, 1993.
  * [27] K. Abe et al.  Measurements of neutrino oscillation in appearance and disappearance channels by the T2K experiment with 6.6E20 protons on target.  Phys.Rev., D91:072010, 2015.
  * [28] Y. Hayato.  A Neutrino Interaction Simulation Program Library NEUT.  Acta Physica Polonica B, 40(9):2477, Sep 2009.
  * [29] A. M. Ankowski and C. Mariani.  Systematic uncertainties in long-baseline neutrino-oscillation experiments.  J. Phys., G44(5):054001, 2017.
  * [30] M. Martini, M. Ericson, and G. Chanfray.  Neutrino energy reconstruction problems and neutrino oscillations.  Phys. Rev., D85:093012, 2012.
  * [31] J. Nieves, F. Sanchez, I. Ruiz Simo, and M. J. Vicente Vacas.  Neutrino Energy Reconstruction and the Shape of the CCQE-like Total Cross Section.  Phys. Rev., D85:113008, 2012.
  * [32] P. Coloma and P. Huber.  Impact of nuclear effects on the extraction of neutrino oscillation parameters.  Phys. Rev. Lett., 111(22):221802, 2013.
  * [33] A. M. Ankowski et al.  Comparison of the calorimetric and kinematic methods of neutrino energy reconstruction in disappearance experiments.  Phys. Rev., D92(7):073014, 2015.
  * [34] O. Lalakulich, U. Mosel, and K. Gallmeister.  Energy reconstruction in quasielastic scattering in the MiniBooNE and T2K experiments.  Phys. Rev., C86:054606, 2012.
  * [35] C. Wilkinson et al.  Testing charged current quasi-elastic and multinucleon interaction models in the NEUT neutrino interaction generator with published datasets from the MiniBooNE and MINER$\nu$A experiments.  Phys. Rev., D93(7):072010, 2016.
  * [36] C. H. Llewellyn Smith.  Neutrino reactions at accelerator energies.  Phys. Rep., 3:261–379, 1972.
  * [37] R. Bradford, A. Bodek, H. Budd, and J. Arrington.  A new parameterization of the nucleon elastic form factors.  Nuclear Physics B - Proceedings Supplements, 159(0):127 – 132, 2006\.
  * [38] E. J. Moniz and R. A. Smith.  Neutrino reactions on nuclear targets.  Nucl. Phys., B101:605–622, 1972.
  * [39] O. Benhar and A. Fabrocini.  Two nucleon spectral function in infinite nuclear matter.  Phys. Rev., C62:034304, 2000.
  * [40] J. Nieves, I. Ruiz Simo, and M. J. Vicente Vacas.  Inclusive Charged–Current Neutrino–Nucleus Reactions.  Phys. Rev., C83:045501, 2011.
  * [41] R. Gran, J. Nieves, F. Sanchez, and M. J. Vicente Vacas.  Neutrino-nucleus quasi-elastic and 2p2h interactions up to 10 GeV.  Phys. Rev., D88(11):113007, 2013.
  * [42] A.A. Aguilar-Arevalo et al.  First measurement of the muon neutrino charged current quasielastic double-differential cross section.  Phys. Rev., D81:092005, 2010.
  * [43] A.A. Aguilar-Arevalo et al.  First measurement of the muon antineutrino double-differential charged current quasi-elastic cross section.  Phys. Rev., D88:032001, 2013.
  * [44] G.A. Fiorentini et al.  Measurement of Muon Neutrino Quasi-Elastic Scattering on a Hydrocarbon Target at $E_{\nu}\sim 3.5$ GeV.  Phys. Rev. Lett., 111:022502, 2013.
  * [45] L. Fields et al.  Measurement of Muon Antineutrino Quasi-Elastic Scattering on a Hydrocarbon Target at $E_{\nu}\sim 3.5$ GeV.  Phys. Rev. Lett., 111:022501, 2013.
  * [46] D. Rein and Sehgal L. M.  Neutrino-excitation of baryon resonances and single pion production.  Annals of Physics, 133(1):79–153, 1981.
  * [47] K. M. Graczyk and J. T. Sobczyk.  Form factors in the quark resonance model.  Phys. Rev. D, 77:053001, Mar 2008.
  * [48] G.M. Radecky, V.E. Barnes, D.D. Carmony, A.F. Garfinkel, M. Derrick, et al.  Study of Single Pion Production by Weak Charged Currents in Low-energy Neutrino $d$ Interactions.  Phys. Rev., D25:1161–1173, 1982.
  * [49] T. Kitagaki, H. Yuta, S. Tanaka, A. Yamaguchi, K. Abe, et al.  Charged Current Exclusive Pion Production in Neutrino Deuterium Interactions.  Phys. Rev., D34:2554–2565, 1986.
  * [50] C. Wilkinson, P. Rodrigues, S. Cartwright, L. Thompson, and K. S. McFarland.  Reanalysis of bubble chamber measurements of muon-neutrino induced single pion production.  Phys. Rev., D90(11):112017, 2014.
  * [51] A. A. Aguilar-Arevalo et al.  Measurement of $\nu_{\mu}$ and $\bar{\nu}_{\mu}$ induced neutral current single $\pi^{0}$ production cross sections on mineral oil at $E_{\nu}\sim\text{O(1 GeV)}$.  Phys. Rev., D81:013005, 2010.
  * [52] A.A.Aguilar-Arevalo et al.  Measurement of $\nu_{\mu}$-induced charged-current neutral pion production cross sections on mineral oil at $E_{\nu}\in 0.5-2.0\,\text{GeV}$.  Phys. Rev., D83:052009, 2011.
  * [53] A. A. Aguilar-Arevalo et al.  Measurement of Neutrino-Induced Charged-Current Charged Pion Production Cross Sections on Mineral Oil at $E_{\nu}\sim 1~{}\text{GeV}$.  Phys.Rev., D83:052007, 2011.
  * [54] B. Eberly et al.  Charged Pion Production in $\nu_{\mu}$ Interactions on Hydrocarbon at $\langle E_{\nu}\rangle$= 4.0 GeV.  Phys. Rev., D92:092008, 2015.
  * [55] Y. Kurimoto et al.  Measurement of Inclusive Neutral Current Neutral Pion Production on Carbon in a Few-GeV Neutrino Beam.  Phys. Rev., D81:033004, 2010.
  * [56] S. Nakayama.  Measurement of single $\pi^{0}$ production in neutral current neutrino interactions with water by a 1.3-GeV wide band muon neutrino beam.  Phys. Lett., B619:255, 2005.
  * [57] D. Rein and L. M. Sehgal.  Coherent $\pi^{0}$ production in neutrino reactions.  Nucl. Phys., B223:29, 1983.
  * [58] D. Rein and L.M. Sehgal.  PCAC and the deficit of forward muons in $\pi^{+}$ production by neutrinos.  Phys. Lett., B657:207–209, 2007.
  * [59] A. Higuera et al.  Measurement of coherent production of $\pi^{\pm}$ in neutrino and antineutrino beams on carbon from $E_{\nu}$ of $1.5$ to $20$ GeV.  Phys. Rev. Lett., 113(26):261802, 2014.
  * [60] Ch. Berger and L. M. Sehgal.  PCAC and coherent pion production by low energy neutrinos.  Phys. Rev., D79:053003, 2009.
  * [61] M. Glück, E. Reya, and A. Vogt.  Dynamical parton distributions revisited.  The European Physical Journal C - Particles and Fields, 5(3):461–470, 1998.
  * [62] A.Bodek and U.K. Yang.  Modeling neutrino and electron scattering cross sections in the few-GeV region with effective LO PDFs.  AIP Conf. Proc., 670:110–117, 2003.
  * [63] T. Sjöstrand.  High-energy-physics event generation with PYTHIA 5.7 and JETSET 7.4.  Computer Physics Communications, 82(1):74 – 89, 1994.
  * [64] P. Adamson et al.  Neutrino and Antineutrino Inclusive Charged-current Cross Section Measurements with the MINOS Near Detector.  Phys. Rev., D81:072002, 2010.
  * [65] E. Wang, L. Alvarez-Ruso, Y. Hayato, K. Mahn, and J. Nieves.  Photon emission in neutral current interactions at the T2K experiment.  Phys. Rev., D92(5):053005, 2015.
  * [66] E. Wang, L. Alvarez-Ruso, and J. Nieves.  Photon emission in neutral current interactions at intermediate energies.  Phys. Rev., C89(1):015503, 2014.
  * [67] R. D. Woods and D. S. Saxon.  Diffuse surface optical model for nucleon-nuclei scattering.  Phys. Rev., 95:577–578, Jul 1954.
  * [68] D.S. Baranov, A.P. Bugorsky, A.A. Ivanilov, V.I. Kochetkov, V.I. Konyushko, et al.  Measurements of the $\nu_{\mu}N$ total cross section at 2-GeV - 30-GeV in SKAT neutrino experiment.  Phys.Lett., B81:255, 1979.
  * [69] M. Day and K. S. McFarland.  Differences in quasi-elastic cross sections of muon and electron neutrinos.  Phys. Rev., D86:053003, 2012.
  * [70] M. Martini, M. Ericson, and G. Chanfray.  Neutrino quasielastic interaction and nuclear dynamics.  Phys. Rev., C84:055502, 2011.
  * [71] K. Abe et al.  Measurements of the T2K neutrino beam properties using the INGRID on-axis near detector.  Nucl. Instrum. Meth., A694:211–223, 2012.
  * [72] N. Abgrall et al.  Time Projection Chambers for the T2K Near Detectors.  Nucl. Instrum. Meth., A637:25–46, 2011.
  * [73] P. A. Amaudruz et al.  The T2K Fine-Grained Detectors.  Nucl. Instrum. Meth., A696:1–31, 2012.
  * [74] S. Fukuda et al.  The Super-Kamiokande detector.  Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 501(2):418 – 462, 2003.
  * [75] K. A. Olive et al.  Review of Particle Physics.  Chin. Phys., C38:090001, (2014) and 2015 update.
  * [76] Z. Maki, M. Nakagawa, and S. Sakata.  Remarks on the unified model of elementary particles.  Prog. Theor. Phys., 28:870–880, 1962.
  * [77] B. Pontecorvo.  Neutrino Experiments and the Problem of Conservation of Leptonic Charge.  Sov. Phys. JETP, 26:984–988, 1968.  [Zh. Eksp. Teor. Fiz.53,1717(1967)].
  * [78] Jiro Arafune, Masafumi Koike, and Joe Sato.  CP violation and matter effect in long baseline neutrino oscillation experiments.  Phys. Rev., D56:3093–3099, 1997.  [Erratum: Phys. Rev.D60,119905(1999)].
  * [79] R. D. Cousins and V. L. Highland.  Incorporating systematic uncertainties into an upper limit.  Nucl. Instrum. Meth., A320:331–335, 1992.
  * [80] G. J. Feldman and R. D. Cousins.  Unified approach to the classical statistical analysis of small signals.  Phys. Rev., D57:3873–3889, Apr 1998.
  * [81] W. K. Hastings.  Monte Carlo sampling methods using Markov chains and their applications.  Biometrika, Vol. 57,:090001, 1970.
  * [82] P. Mills.  Efficient statistical classification of satellite measurements.  International Journal of Remote Sensing, 32:6109–6132, 2011.
  * [83] K. Cranmer.  Kernel estimation in high-energy physics.  Computer Physics Communications, 136:198–207, 2001.
  * [84] F. James and M. Roos.  Minuit - a system for function minimization and analysis of the parameter errors and correlations.  Computer Physics Communications, 10:343–367, 1975.
  * [85] L. Demortier.  P values and nuisance parameters.  2008\.
  * [86] P. Adamson et al.  Measurement of the neutrino mixing angle $\theta_{23}$ in NOvA.  Submitted to: Phys. Rev. Lett., 2017.
  * [87] P. Adamson et al.  Measurement of neutrino and antineutrino oscillations using beam and atmospheric data in minos.  Phys. Rev. Lett., 110:251801, Jun 2013.
  * [88] R. Wendell.  Atmospheric Neutrino Oscillations at Super-Kamiokande.  PoS, ICRC2015:1062, 2016.
  * [89] M. G. Aartsen et al.  Neutrino oscillation studies with IceCube-DeepCore.  Nucl. Phys., B908:161–177, 2016.
  * [90] H. Jeffreys.  The Theory of Probability.  OUP Oxford, 1998.
  * [91] Abe, K. et al.  T2K public data.  http://t2k-experiment.org/.