Search for CP violation in Neutrino and Antineutrino Oscillations by the T2K experiment with $2.2\times10^{21}$ protons on target

The T2K experiment measures muon neutrino disappearance and electron neutrino appearance in accelerator-produced neutrino and antineutrino beams. With an exposure of $14.7(7.6)\times 10^{20}$ protons on target in neutrino (antineutrino) mode, 89 $\nu_e$ candidates and 7 anti-$\nu_e$ candidates were observed while 67.5 and 9.0 are expected for $\delta_{CP}=0$ and normal mass ordering. The obtained $2\sigma$ confidence interval for the $CP$ violating phase, $\delta_{CP}$, does not include the $CP$-conserving cases ($\delta_{CP}=0,\pi$). The best-fit values of other parameters are $\sin^2\theta_{23} = 0.526^{+0.032}_{-0.036}$ and $\Delta m^2_{32}=2.463^{+0.071}_{-0.070}\times10^{-3} \mathrm{eV}^2/c^4$.

Introduction. -Observation of neutrino oscillations has established that each of the three flavor states of neutrinos is a superposition of at least three mass eigenstates, m 1 , m 2 and m 3 [1][2][3][4]. As a consequence of threegeneration mixing, the flavor-mass mixing matrix, the PMNS matrix [5,6], can have an irreducible imaginary component, and CP symmetry can be violated in neutrino oscillations, analogous to the case of the quark sector. The PMNS matrix is parametrized by three mixing angles θ 12 , θ 13 and θ 23 , and one CP violation phase, δ CP , which gives rise to asymmetries between neutrino oscillations and antineutrino oscillations if sin δ CP = 0. The magnitude of CP violation is determined by the invariant J CP = 1 8 cos θ 13 sin 2θ 12 sin 2θ 23 sin 2θ 13 sin δ CP ≈ 0.033 sin δ CP [7,8] and could be large compared to the quark sector value (J CP ≈ 3 × 10 −5 ). The most feasible way to probe δ CP is by measuring the appearance of electron (anti)neutrinos ( ( ν ) e ) by using accelerator-produced muon (anti)neutrino ( ( ν ) µ ) beams. T2K has reported that the CP conservation hypothesis (δ CP = 0, π) is excluded at 90% confidence level (C.L.) using the data collected up to May 2016 [9,10]. Since then, the neutrino mode data set has doubled, and the electron neutrino and antineutrino event selection efficiencies have increased by 30% and 20%, respectively. In this Letter, we report new results on δ CP , sin 2 θ 23 and ∆m 2 (∆m 2 32 ≡ m 2 3 − m 2 2 for normal or ∆m 2 13 ≡ m 2 1 − m 2 3 for inverted mass ordering) obtained by analyzing both muon (anti)neutrino disappearance and electron (anti)neutrino appearance data collected up to May 2017 using a new event selection method.
The T2K experiment [11]. -The 30 GeV proton beam from the J-PARC accelerator strikes a graphite target to produce charged pions and kaons which are focused or defocused by a system of three magnetic horns. The focused charge is defined by the horn current direction, producing either a muon neutrino or antineutrino beam from the focused secondaries decaying in the 96 m long decay volume. An on-axis near detector (INGRID) and a detector 2.5 • off the beam axis (ND280) sample the unoscillated neutrino beam 280 m downstream from the target station and monitor the beam direction, composition, and intensity. The off-axis energy spectrum peaks at 0.6 GeV and has significantly less ( ν ) e contamination at the peak energy and less high energy neutrino flux than on axis. The Super-Kamiokande (SK) 50 kt water-Cherenkov detector [12], as a far detector, samples the oscillated neutrino beam 2.5 • off axis and 295 km from the production point.
Data set. -The results presented here are based on data collected from Jan 2010 to May 2017. The data sets include a beam exposure of 14.7 × 10 20 protons on target (POT) in neutrino mode and 7.6 × 10 20 POT in antineutrino mode for the far-detector (SK) analysis, and an exposure of 5.8×10 20 POT in neutrino mode and 3.9×10 20 POT in antineutrino mode for the near-detector (ND280) analysis.
Analysis strategy. -Oscillation parameters are determined by comparing model predictions with observations at the near and far detectors. The neutrino flux is modeled based on a data-driven simulation. The neutrinonucleus interactions are simulated based on theoretical models with uncertainties estimated from data and models. The flux and interaction models are refined by the observation of the rate and spectrum of charged-current (CC) neutrino interactions by ND280. Since ND280 is magnetized, wrong-sign contamination in the beam can be estimated from charge-selected near-detector samples. The prediction of the refined model is compared with the observation at SK to estimate the oscillation parameters. Overall analysis method is the same as in previous T2K results [10], but this analysis uses improved theoretical models to describe neutrino interactions and a new reconstruction algorithm at SK, which improves signal/background discrimination and allows an expanded fiducial volume.
Neutrino flux model. -A data-driven simulation is used to calculate the neutrino and antineutrino fluxes and their uncertainties at each detector, including correlations [10,13]. The interactions of hadrons in the target and other beamline materials are tuned using external thin-target hadron-production data, mainly measurements of 30 GeV protons on a graphite target by the NA61/SHINE experiment [14]. The simulation reflects the proton beam condition, horn current and neutrino beam-axis direction as measured by monitors. Near the peak energy, and in the absence of oscillations, 97.2% (96.2%) of the (anti)neutrino mode beam is ( ν ) µ . The remaining components are mostlyν µ (ν µ ); contamination of ( ν ) e is only 0.42% (0.46%). The dominant source of systematic error in the flux model is the uncertainty of the hadron-production data. Some of the beamline conditions are different depending on time. Stability of the neutrino flux has been monitored by INGRID throughout the whole data taking period. The flux covariance matrix was constructed by removing the near-far correlations for time-dependent systematics for the period during which ND280 data was not used in this analysis. While the flux uncertainty is approximately 9% at the peak energy, its impact on oscillation parameter uncertainties, given that the near and far detector measurements sample nearly the same flux, is significantly smaller.
Neutrino interaction model. -Events are simulated with the NEUT [15] neutrino interaction generator. The dominant charged-current quasi-elastic (CCQE)-like interaction (defined as those with a charged lepton, and no pions in the final state) is modeled with a relativistic Fermi gas (RFG) nuclear model including longrange correlations using the random phase approximation (RPA) [16]. The 2p-2h model of Nieves et al. [17,18] predicts multi-nucleon contributions to CCQE-like processes. These can be divided into meson exchange current (∆-like) contributions, which include both diagrams with an intermediate ∆ and contributions from pion inflight and pion contact terms (see Ref. [17] for details), and contributions from interactions with correlated nucleonnucleon pairs (non-∆-like), which introduce different biases in the reconstructed neutrino energy, E rec , calculated assuming QE scattering [10] 1 . New parameters are introduced to vary the relative contribution of ∆-like and non-∆-like terms 2 for 12 C and 16 O, with a 30% correlation between the two nuclei. The total 2p-2h normalization is varied separately for ν andν with flat priors. There is an additional uncertainty on the ratio of 12 C to 16 O 2p-2h normalizations, with a 20% uncertainty. The Q 2 dependence of the RPA correction is allowed to vary by the addition of four variable parameters designed to span the total theoretical uncertainty in the Q 2 dependence [19,20]. Processes producing a single pion and one or more nucleons in the final state are described by the Rein-Sehgal model [21]. Parameters describing the ∆ axial form factor and single pion production not through baryon resonances are tuned to match D 2 measurements [22][23][24] in a method similar to [25]. Production of pions in coherent inelastic scattering is described by a tuned model of Rein-Sehgal [26] which agrees with recent measurements [27,28]. As in [10], differences between muon-and electron-neutrino interactions occur because of final-state lepton mass and radiative corrections, and are largest at low energies. To account for this, we add a 2% uncorrelated uncertainty for each of electron neutrino and antineutrino cross sections relative to those of muons (σ CC (ν e )/σ CC (ν µ ) and σ CC (ν e )/σ CC (ν µ )), and another 2% uncertainty anticorrelated between the two ratios [29]. The cross-section parametrization is otherwise as described in [10], with the exception of variations of the nucleon removal energy, E b , by 25(27) +18 −9 MeV for 12 C( 16 O) [30]. Some systematic uncertainties are not easily implemented by varying model parameters. These are the subjects of "simulated data" studies, where simulated data generated from a variant model is analyzed under the assumptions of the default model. Studies include varying E b , replacing the RFG model with a local Fermi gas model [17] or a spectral function model [31], changing the 2p-2h model to an alternate one [32] or fixing the 2p-2h model to be fully "∆-like" or "non-∆-like", varying the axial nucleon form factor to allow more realistic high Q 2 uncertainties [33,34], and using an alternative single pion production model described in [35]. Additional simulated data studies, based on an excess observed at low muon momentum (p µ ≤ 400 MeV) and moderate angle (0.6 ≤ cos θ µ ≤ 0.8) in the near detector, quantified possible biases in neutrino energy reconstruction by modeling this as an additional ad hoc interaction under hypotheses that it had 1p-1h, ∆-like 2p-2h or non-∆-like 2p-2h kinematics. Finally, a discrepancy in the pion kinematic spectrum observed at the near detector motivated a simulated data study to check the impact on the signal samples at SK.
Fits to these simulated data sets showed no significant biases in δ CP or sin 2 θ 13 ; however biases in ∆m 2 comparable to the total systematic uncertainty were seen for most data sets. This bias was accounted for by adding an additional source of uncertainty into the confidence intervals in ∆m 2 , as described later. As well as biases in ∆m 2 , fits to the varied E b simulated data sets also showed biases in sin 2 θ 23 comparable to the total systematic uncertainty. To account for this bias, an additional degree of freedom was added to the fit, which allows the model to replicate the spectra expected at the far detector when E b is varied. After the addition of these additional uncertainties, fits to the simulated data sets no longer show biases that are significant compared to the total systematic error.
Fit to the near detector data. -Fitting the unoscillated spectra of CC candidate events in ND280 constrains the systematic parameters in the neutrino flux and cross-section models [11]. The CC samples are composed of reconstructed interactions in one of the two Fine-Grained Detectors (FGD) with particle tracking through Time Projection Chambers (TPC) interspersed among the FGDs. While both FGDs have active layers of segmented plastic scintillator, the second FGD (FGD2) additionally contains six water-target modules, allowing direct constraints of neutrino interactions on H 2 O, the same target as SK. The ND280 event selection is unchanged from the previous T2K publication [36]. The CC inclusive events are separated into different samples depending on the FGD in which the interaction occurred, the beam mode, the muon charge and the final-state pion multiplicity. The negative muon candidates from data taken in neutrino mode are divided into three samples per FGD based on reconstructed final state topologies: no pion candidate (CC0π), one π + candidate (CC1π), and all the other CC event candidates (CC Other), dominated respectively by the CCQE-like process, CC single pion production, and deep inelastic scattering. In antineutrino mode, positively-and negatively-charged muon tracks are used to define CC event candidates, which are distributed in two topologies: those with only a single muon track reconstructed in the TPC (CC 1-track), and those with at least one other track reconstructed in the TPC (CC N-track). All event samples are binned according to the candidate's momentum p µ and cos θ µ , where θ µ is the angle between the track direction and the detector axis. A binned likelihood fit to the data is performed assuming a Poisson-distributed number of events in each bin with an expectation computed from the flux, cross-section and ND280 detector models. The near detector systematic and flux parameters are marginalized in estimating the far detector flux and cross-section parameters and their covariances. The uncertainties on neutral current and ν e interactions cannot be constrained by the current ND280 selection, therefore the fit leaves the related parameters unconstrained. Figure 1 shows data, pre-fit and post-fit Monte-Carlo p µ distributions for the FGD2 CC0π sample. A deficit of 10%-15% in the pre-fit predicted number of events is observed, which is consistent with the previous T2K publications [36]. In this previous analysis, the simulated flux was increased to compensate the deficit. This is now resolved by the new RPA treatment, by increasing the low Q 2 part of the cross section. Good agreement is observed between the post-fit model and the data, with a p-value of 0.473, which is better agreement than in the previous T2K publication [36], partly due to the modified cross-section parametrization. The fit to the ND280 data reduces the flux and the ND280-constrained interaction model uncertainties on the predicted event rate at the far detector from 11-14% to 2.5-4% for the different samples.
Far detector event selection and data. -Events at the far detector are required to be time-coincident with the beam and to be fully contained in the SK inner detector, by requiring limited activity in the outer detector. A newly-deployed Cherenkov-ring reconstruction algorithm, previously used only for neutral current (NC) π 0 background suppression [37], is used to classify events into five analysis samples, enriched in: ( ν ) µ CCQE; ( ν ) e CCQE; and ν e CC1π + where the π + is below Cherenkov threshold. The reconstruction algorithm uses all the information in an event by simultaneously fitting the time and charge of every photosensor in the detector. This results in an improved resolution of reconstructed quantities and particle identification. The fiducial volume is defined for each sample in terms of the minimum distance between the neutrino interaction vertex and the detector wall (wall ) and the distance from the vertex to the wall in the direction of propagation (towall ). These criteria are optimized taking into account both statistical and systematic uncertainties, with the systematic parameters related to ring-counting and e/µ, e/π 0 and µ/π + separation being constrained in a fit to SK atmospheric data. Other systematic uncertainties related to the modeling of the far detector are estimated using non-neutrino control samples. Detector systematic error covariances between samples and bins for the oscillation analysis are constructed in the same way as was described in previous T2K publications [37].
The π 0 and π + NC suppression cuts are optimized by running a simplified oscillation analysis [38] on a simulated data set and choosing the criteria that minimize the uncertainty on the oscillation parameters.
All selected events are required to have only one Cherenkov ring. For the ( ν ) µ CCQE-enriched samples, the single-ring events are further required to have wall > 50 cm and towall > 250 cm; be classified as µ-like by the µ/e separation cut; have a reconstructed momentum greater than 200 MeV/c; have up to one decay-electron candidate; and satisfy the π + rejection criterion. After these selection cuts are applied, 240 events are found in the neutrino-mode data and 68 in antineutrino-mode data, with an expectation of 261.6 and 62.0, respectively, for sin 2 θ 23 = 0.528 and ∆m 2 32 = 2.509 × 10 −3 eV 2 /c 4 . The E rec distributions for the data and best-fit Monte Carlo are shown in Figure 2.
The ( ν ) e CCQE-enriched samples contain e-like events with no decay electron candidates, that pass the π 0 rejection cut, have wall > 80 cm, towall > 170 cm, momentum > 100 MeV/c, and a reconstructed neutrino energy (E rec ) lower than 1250 MeV. E rec is calculated from the lepton momentum and angle assuming CCQE kinematics. The ν e CC1π + -enriched sample has the same selection criteria with the exception of the fiducial volume criteria, which are wall > 50 cm and towall > 270 cm, and the requirement of one decay electron candidate in the event, from which the presence of a π + is inferred. Like in the case of the CCQE-enriched samples, E rec for the ν e CC1π + sample is calculated from the outgoing electron kinematics, except in this case the ∆ ++ mass is assumed for the outgoing nucleon. Event yields for these samples are compared to Monte-Carlo predictions in Table II and their E rec distributions are shown in Figure 3. Compared to previous T2K publications, the optimized event selection criteria are expected to increase the acceptance for ( ν ) µ CCQE events by 15% with a 50% reduction of the NC1π + background; to increase the ( ν ) e CC events acceptance by 20% with similar purity to previous analyses; and to increase the ν e CC1π + acceptance by 33% with a 70% reduction in background caused by particle misidentification. A summary of the systematic uncertainties on the predicted event rates at SK is given in Table I. Oscillation analysis.-A joint maximum-likelihood fit to five far-detector samples constrains the oscillation parameters sin 2 θ 23 , ∆m 2 , sin 2 θ 13 and δ CP . Oscillation probabilities are calculated using the full three-flavor oscillation formulas [39] including matter effects, with a crust density of ρ = 2.6 g/cm 3 [40].
Priors for the flux and interaction cross-section parameters are obtained using results from a fit to the neardetector data. Flat priors are chosen for sin 2 θ 23 , |∆m 2 | and δ CP . The two mass orderings are each given a probability of 50%. In some fits a flat prior is also chosen for sin 2 2θ 13 ; whereas, in fits that use reactor neutrino measurements, we use a Gaussian prior of sin 2 2θ 13 = 0.0857 ± 0.0046 [41]. The θ 12 and ∆m 2 21 parameters have negligible effects and are constrained by Gaussian priors from the PDG [41].
Using the same procedure as [10], we integrate the product of the likelihood and the nuisance priors to obtain the marginal likelihood, which does not depend on the nuisance parameters. We define the marginal likelihood ratio as −2∆lnL = −2 ln(L/L max ), where L max is the maximum marginal likelihood.
Using this statistic, three independent analyses have been developed. The first and second analyses provide confidence intervals using a hybrid Bayesian-frequentist approach [42]. The third analysis provides credible intervals using the posterior probability distributions calculated with a fully Bayesian Markov chain Monte Carlo method [43]. This analysis also simultaneously fits both near-and far-detector data, which validates the extrapolation of nuisance parameters from the near to far detector. For all three analyses, the ( ν ) µ samples are binned by E rec . The first and third analyses bin the three ( ν ) e samples in E rec and lepton angle, θ, relative to the beam, while the second analysis uses lepton momentum, p, and θ. All three analyses give consistent results.
Expected event rates for various values of δ CP and mass ordering are shown in Table II. An indication of the sensitivity to δ CP can be seen from the ∼ 20% variation in the predicted total event rate between the CP conserved case (δ CP = 0, π) and when CP is maximally violated. The ( ν ) µ event rates are negligibly affected by the mass ordering, whereas the ( ν ) e rates differ by ∼ 10% between mass orderings. In the ν e CC1π + sample we see 15 events when we expected 6.9 for δ CP = −π/2 and normal ordering. The p-value to observe an upwards or downwards fluctuation of this significance in any one of the five samples used is 12%. The p-value to observe the data given the posterior expectation across all samples is greater than 35%. FIG. 5. The 68% (90%) constant −2∆lnL confidence regions in the sin 2 θ13-δCP plane using a flat prior on sin 2 (2θ13), assuming normal (black) and inverted (red) mass ordering. The 68% confidence region from reactor experiments on sin 2 θ13 is shown by the yellow vertical band. Fits to determine either one or two of the oscillation parameters are performed, while the other parameters are marginalized. The constant −2∆lnL method is then used to set confidence regions [41]. Confidence regions in the |∆m 2 |-sin 2 θ 23 plane (Fig. 4) were first computed for each mass ordering separately using the reactor measurement prior on sin 2 θ 13 . The likelihood used to generate these confidence regions is convolved with a Gaussian function in the ∆m 2 direction. The standard deviation of this Gaussian is 3.5 × 10 −5 eV 2 /c 4 , which is the quadrature sum of the biases on ∆m 2 seen in the fits to the simulated data sets.
Confidence regions in the sin 2 θ 13 -δ CP plane were calculated, without using the reactor measurement prior on sin 2 (2θ 13 ), for both the normal and inverted orderings (Fig. 5). T2K's measurement of sin 2 θ 13 agrees well with the reactor measurement.
Confidence intervals for δ CP were calculated using the Feldman-Cousins method [45], marginalized over both mass orderings simultaneously, from a fit using the reactor measurement prior. The best fit value is δ CP = −1.87(−1.43) for the normal (inverted) ordering, which is close to maximal CP violation (Fig. 6). The δ CP confidence intervals at 2σ (95.45%) are (−2.99, −0.59) for normal ordering and (−1.81, −1.01) for inverted ordering. Both intervals exclude the CP -conserving values of 0 and π. The Bayesian credible interval at 95.45% is (−3.02,−0.44), marginalizing over the mass ordering. The normal ordering is preferred with a posterior probability of 87%.
Sensitivity studies show that, if the true value of δ CP is −π/2 and the mass ordering is normal, 22% of simulated experiments exclude δ CP = 0 and π at 2σ C.L.
Conclusions. -T2K has constrained the leptonic CP violation phase (δ CP ), sin 2 θ 23 , ∆m 2 and the posterior probability for the mass orderings with additional data and with an improved event selection efficiency. The 2σ (95.45%) confidence interval for δ CP does not contain the CP -conserving values of δ CP = 0, π for either of the mass orderings. The current result is predominantly limited by statistics. T2K will accumulate 2.5 times more data, thereby improving sensitivity for the relevant oscillation parameters. The data related to the measurement and results presented in this Letter can be found in [46].