Search for Electron Antineutrino Appearance in a Long-Baseline Muon Antineutrino Beam

Electron antineutrino appearance is measured by the T2K experiment in an accelerator-produced antineutrino beam, using additional neutrino beam operation to constrain parameters of the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) mixing matrix. T2K observes 15 candidate electron antineutrino events with a background expectation of 9.3 events. Including information from the kinematic distribution of observed events, the hypothesis of no electron antineutrino appearance is disfavored with a significance of 2 . 40 σ and no discrepancy between data and PMNS predictions is found. A complementary analysis that introduces an additional free parameter which allows non-PMNS values of electron neutrino and antineutrino appearance also finds no discrepancy between data and PMNS predictions.

Electron antineutrino appearance is measured by the T2K experiment in an accelerator-produced antineutrino beam, using additional neutrino beam operation to constrain parameters of the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) mixing matrix. T2K observes 15 candidate electron antineutrino events with a background expectation of 9.3 events. Including information from the kinematic distribution of observed events, the hypothesis of no electron antineutrino appearance is disfavored with a significance of 2.40σ and no discrepancy between data and PMNS predictions is found. A complementary analysis that introduces an additional free parameter which allows non-PMNS values of electron neutrino and antineutrino appearance also finds no discrepancy between data and PMNS predictions. Introduction.-The observation of neutrino oscillations has established that each neutrino flavor state (e, μ, τ) is a superposition of at least three mass eigenstates (m 1 , m 2 , m 3 ) [1 -4]. The phenomenon of oscillation is modeled by a three-generation flavor-mass mixing matrix, called the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix [5,6]. With the discovery of nonzero θ 13 and the explicit observation of ν μ to ν e appearance oscillation [7], it is now crucial to test the PMNS framework and establish if it is sufficient to explain all neutrino and antineutrino oscillation observations. One such test is to search for the CP-reversed appearance oscillation ofν μ toν e . A search for this process in the Tokai-to-Kamioka (T2K) experiment was reported in Ref. [8], and recent results from the NOvA experiment show a significance of 4.4σ [9]. In this Letter, we report a search for electron antineutrino appearance at the T2K experiment with an improved event selection and a dataset more than a factor of 2 larger than previous T2K results.
The T2K experiment.-The T2K experiment [10] begins with a 30 GeV proton beam from the J-PARC main ring striking a graphite target, producing pions and kaons. These charged hadrons are focused by a system of three magnetic horns to decay in a 96 m decay volume. Positively charged hadrons are focused to produce a beam of predominantly neutrinos ("neutrino mode"); negatively charged hadrons are focused for a beam of predominantly antineutrinos ("antineutrino mode").
An unmagnetized on-axis near detector (INGRID) and a magnetized off-axis (2.5°) near detector (ND280) sample the unoscillated neutrino beam 280 m downstream from the target station and monitor the beam direction, composition, and intensity and constrain neutrino interaction properties. The unmagnetized Super-Kamiokande (SK) 50 kt water-Cherenkov detector is the T2K far detector, and samples the oscillated neutrino beam 2.5°off axis and 295 km from the production point.
The analysis presented here uses data collected from January 2010 to June 2018. The dataset has an exposure at SK of 1.63 × 10 21 protons on target (POT) in antineutrino mode, with an additional dataset of 1.49 × 10 21 POT in neutrino mode used to constrain PMNS oscillation parameters acting as systematic uncertainties in the analysis. The ND280 detector uses an exposure of 0.58 × 10 21 POT in neutrino mode and 0.39 × 10 21 POT in antineutrino mode.
Analysis strategy.-The significance ofν e appearance is evaluated by introducing the parameter β, which multiplies the PMNS oscillation probability Pðν μ →ν e Þ: The analysis is performed allowing both β ¼ 0 and β ¼ 1 to be the null hypothesis, where both hypotheses fully Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI. Funded by SCOAP 3 . account for uncertainties in the values of the oscillation and systematic parameters. The β ¼ 0 case determines ifν e events can be seen above background in the experiment, and the β ¼ 1 case determines if the data is consistent with PMNS. Two analyses are performed on each hypothesis to obtain corresponding p values: one uses only the number of events ("rate only"); while the other also uses information from the kinematic variables of events ("rate þ shape").
The total number of candidateν e events in the antineutrino beam mode is used as the test statistic to calculate the rate-only p value. The test statistic is used to calculate the rate þ shape p value, where the χ 2 values are calculated by marginalizing over all systematic and oscillation parameters, including the mass ordering.
In both analyses, other data samples-ν μ -like and ν e -like in neutrino beam mode andν μ -like in antineutrino beam mode-are used to constrain other PMNS oscillation parameters, as in other T2K analyses [11]. A complementary analysis allows β to be a continuous free parameter with limits between 0 and infinity. In this analysis only, in addition to β multiplying P PMNS ðν μ →ν e Þ as in Eq. (1), the probability P PMNS ðν μ → ν e Þ is multiplied by a factor 1=β. This formulation-slightly different from above-was chosen for its property of anticorrelation in shifting probability between neutrinos and antineutrinos. The extra degree of freedom allows the fit to explore areas away from the PMNS constraint to more accurately reflect the information given by the data. Credible interval contours in the Pðν μ → ν e Þ and Pðν μ →ν e Þ parameter space, the main result of the analysis, are then compared against T2K data fit with β fixed to 1 to test the compatibility between the T2K data and the PMNS model constraining the standard fit.
Neutrino beam flux.-The primary signal datasets were taken in antineutrino mode. The flux was predicted by a Monte Carlo (MC) simulation incorporating the FLUKA2011 interaction model [12] tuned to the results of recent external hadron production experiments including the NA61/SHINE experiment at CERN [13][14][15]. The INGRID detector is used to monitor the beam axis direction and total flux stability.
The resultant flux model [16][17][18] estimates unoscillated neutrino and antineutrino fluxes at all detectors as well as their uncertainties and correlations. The flux at ND280 and SK peaks at 600 MeV, where 96.2% of the beam is composed ofν μ and 0.46%ν e . The remainder of the beam is almost entirely ν μ . This wrong sign contamination is greater in antineutrino mode than neutrino mode.
Neutrino interaction model.-The NEUT (v5.3.3) neutrino interaction generator [19] is used to generate simulated neutrino events. The model used is described in Refs. [8] and [11]. The most relevant contributions for this analysis are highlighted here.
The dominant neutrino-nucleus interaction topology near 600 MeV, charged current quasielastic (CCQE)-like, is defined as an interaction with one charged lepton and zero pions in the final state. The nucleus is modeled with a relativistic Fermi gas modified by a random phase approximation (RPA) to account for long-range correlations [20]. A multinucleon component is included with the Nieves 2p-2h model [21,22], which contains both meson exchange current (Δ-like) and correlated nucleon pair (non-Δ-like) contributions. Parameters representing systematic uncertainties for the CCQE-like mode include the nucleon axial mass, M QE A ; the Fermi momentum for 12 C and 16 O; the 2p-2h normalization term for ν andν separately; four parameters controlling the RPA shape as a function of Q 2 ; and the relative contributions of the Δ-like and non-Δ-like contributions to 2p-2h in 12 C and 16 O. The RPA parameters have Gaussian priors to cover the theoretical shape uncertainty given in [23,24], and the 2p-2h shape contribution has a 30% correlation between 12 C and 16 O; all other priors are uniform. Other neutrino-nucleus processes are subdominant, and their rates are constrained via appropriate uncertainties.
Differences between muon-and electron-neutrino interactions are largest at low energies and occur because of final-state lepton mass and radiative corrections. A 2% uncorrelated uncertainty is added for each of the electron neutrino and antineutrino cross sections relative to those of muons and another 2% uncertainty anticorrelated between the two ratios [25].
Some systematic uncertainties are not easily included by varying model parameters. These are the subjects of "simulated data" studies, where simulated data generated from a variant model are analyzed under the assumptions of the default model. The model variations that produce the largest changes in theν e far detector spectra are an alternate single resonant pion model [26], and ad hoc models driven by observed discrepancies in the near detector kinematic spectra, where the discrepancy is modeled as having either 1p − 1h, 2p − 2h − Δ-like, and 2p-2h-non-Δ-like kinematics. None of the variant models studied showed differences in the sensitivity values at greater than the 0.1σ level.
Near detector data constraints.-The ND280 detector is used to fit unoscillated samples of charged current (CC) muon neutrino interaction events to constrain flux and cross section systematic uncertainties for the signal and background models of SK events. The samples-unchanged from Ref. [11]-are selected from events that begin in one of two fine-grained detectors (FGDs) and produce tracks that enter the time-projection chambers, which are interleaved with the FGDs. Both FGDs are composed of layers of bars of plastic scintillator, and the more downstream FGD additionally has panels of water interleaved between layers of scintillator.
In neutrino beam mode, in each FGD, the CC events (defined as containing negatively charged muonlike track) PHYSICAL REVIEW LETTERS 124, 161802 (2020) 161802-4 are split into three subsamples: a CC0π sample, with zero pions in the final state, enhanced in CCQE-like interactions; a CC1π þ sample, with one π þ in the final state, enhanced in resonant pion interactions; and a CC other sample, containing all other CC events. In antineutrino beam mode, in each FGD, there are selected interactions with positively charged muons (ν -like) and negatively charged muons (ν-like). The latter constrains the wrong-sign contamination, which is larger in antineutrino beam mode. Each of these selections is divided into two topologies: containing a single track and containing multiple tracks.
All samples are fit simultaneously and are binned in lepton momentum, p μ , and lepton angle, cos θ μ relative to the average beam neutrino direction. 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 fit returns central values and correlated uncertainties for systematic uncertainty parameters that are constrained by the near detector, marginalizing over near detector flux and detector systematic parameters. Some uncertainties on neutral current and ν e events cannot be constrained by these ND280 samples and those parameters are passed to the appearance analysis with their original prior.
The MC prediction before fitting underestimates the data by 10%-15%, consistent with previous T2K analyses. The agreement between the MC prediction after fitting and data is good, with a p value of 0.473. The fit to the ND280 data reduces the flux and the ND280-constrained interaction model uncertainties on the predicted electron antineutrino sample event rate at the far detector from 14.6% to 7.6%.
ν e SK selection.-Unlike in the previous analysis, SK events are reconstructed and selected using the new reconstruction algorithm described in Ref. [27]. Aν e event candidate in SK must meet the following criteria: (i) it is within the beam time window as determined from a GPS time stamp, and its Cherenkov light is fully contained in the SK inner detector, with minimal outer-detector activity; (ii) the reconstructed vertex is at least 80 cm from the innerdetector wall; (iii) only one Cherenkov ring candidate is found in the reconstruction and the ring is identified as electronlike; (iv) the distance from the vertex to the detector wall is greater than 170 cm along the track direction; (v) the visible energy in the event is greater than 100 MeV; (vi) there is no evidence of delayed activity consistent with a stopped muon decay; (vii) the reconstructed energy under a quasielastic scattering hypothesis is less than 1250 MeV; (viii) the ring is inconsistent with a π 0 decay hypothesis.
These reconstruction cuts have an efficiency of 71.5% forν e events that satisfy the fully contained and fiducial requirements. The new event selection increases the yield ofν e signal by approximately 20% compared to the previous analysis, primarily due to the new fiducial cuts, with no loss of purity. Assuming oscillation parameter values near the best fit of previous T2K analyses of sin 2 θ 23 ¼ 0.528, sin 2 θ 13 ¼ 0.0212, sin 2 θ 12 ¼ 0.304, Δm 2 32 ¼2.509×10 −3 eV 2 =c 4 , Δm 2 21 ¼ 7.53 × 10 −5 eV 2 =c 4 , δ CP ¼ −1.601, normal ordering and β ¼ 1, the total expected background is 9.3 events including 3.0 ν e interactions resulting from oscillations of ν μ in the beam. The remaining major sources of background are intrinsic ν e and ν e in the beam (4.2 events) and neutral-current interactions (2.1 events). With the oscillation parameters above, a signal yield of 7.4 events is expected, for a total prediction of 16.8 events. Figure 1 shows the 15 observed data events superimposed on a prediction generated using the above oscillation parameter values. ν e appearance.-Theν e appearance p values are calculated by considering the rate-only and rate þ shape test statistics of an ensemble of 2 × 10 4 pseudoexperiments. Each pseudoexperiment is generated by randomizing systematic parameters-including oscillation parametersand applying statistical fluctuations. Four control samples, ν mode single-ring e-like and ν e CC1π-like (single-ring e-like accompanied by electron decay) and both ν andν mode single-ring μ-like, are used to constrain the distribution of oscillation parameters of the pseudoexperiments. The four control samples of many pseudoexperiments are compared to data, and rejection sampling is used to select 2 × 10 4 that are most probable, according to data. The systematic parameters are then marginalized over using a numeric integration technique (with 2 × 10 5 samples of the systematic parameter space) when calculating the rate þ shape test statistic. Both the number of pseudoexperiments and the number of points used for the numerical integration were studied and selected to ensure p value stability.
When producing the pseudoexperiments and marginalizing over systematic uncertainties, Gaussian prior probabilities on the following oscillation parameters are used: sin 2 2θ 12 (0.846 AE 0.021), Δm 2 21 ðð7.53 AE 0.18Þ × 10 −5 eV 2 =c 4 Þ, and sin 2 2θ 13 (0.0830 AE 0.0031) [28]. The mass ordering is randomized with a probability of 0.5 for NO and 0.5 for IO. The other PMNS parameters are randomized using uniform prior probabilities with limits set based on previous experiments. Systematic parameters are randomized according to the constraints set by the near detector fit.
When predicted distributions are compared to data, a binned Poisson likelihood is used for all five SK data samples. The e-like samples use a 2D distribution in the reconstructed neutrino energy, E rec , and the reconstructed neutrino angle with respect to the average beam direction, θ. The μ-like samples use a 1D distribution in the reconstructed neutrino energy.
For the rate þ shape analysis, the likelihood for a pseudoexperiment is defined as the product of the likelihoods of theν mode single-ring e-like sample, λ¯ν e , and the control samples, λ c . The test statistic is then calculated as in Eq. (3), by averaging this likelihood over samples of the systematic parameter space, a i . When the generated distribution of the test statistic is calculated, λν e is compared to the pseudoexperiment data, E, and λ c is compared to data, D; when the test statistic for the real data is calculated, both likelihoods are compared to data, An independent, complementary analysis uses the kinematic variable of outgoing lepton momentum, p l instead of reconstructed neutrino energy, and additionally uses weighting of pseudoexperiments instead of rejection sampling. Both analyses were found to give consistent test statistic distributions and therefore p values.
The distributions of the rate-only and rate þ shape test statistics for the β ¼ 0 and β ¼ 1 hypotheses are shown in Fig. 2. These distributions are integrated from the data test statistic to obtain right(left)-tailed p values for the β ¼ 0ð1Þ hypothesis. The observed number of events in theν mode single-ring e-like sample in SK was 15, compared to a prediction of 16.8. The observed data Δχ 2 value in the rate þ shape analysis was 3.811 and the prediction was 6.3. The resulting p values are shown in Table I. Both the rate-only and rate þ shape analyses disfavor the no-ν eappearance hypothesis (β ¼ 0) more than the PMNSν e appearance hypothesis (β ¼ 1). Compared to the prediction, a slightly weaker exclusion of the noν e appearance hypothesis (β ¼ 0) is observed due to observing fewer events than expected. The rate þ shape analysis gives a stronger observed exclusion of both hypotheses than the rate-only analysis, due to the extra shape information used to discredit each hypothesis.
Continuous β.-A complementary analysis allows β to be a free parameter, which allows for a continuum of non-PMNS models, rather than only the single β ¼ 0 no-ν e -appearance case. The impact of this analysis is shown in the parameter space of Pðν μ → ν e Þ vs Pðν μ →ν e Þ, and in the ν e vsν e event rate space. Varying δ CP at a fixed energy creates an ellipse with a negatively sloping major axis in the biprobability phase space. Switching the mass ordering shifts the center of the ellipse along the Pðν μ → ν e Þ ¼ −Pðν μ →ν e Þ axis. The other oscillation parameters shift the ellipses along the identity line in the biprobability space. Two ellipses are shown in the left panel of Fig. 3 in orange and brown, with the input oscillation parameter values taken from the β ¼ 1 fit; the eccentricity of the ellipses is very large for the T2K experiment, which makes them appear like lines. In the ellipses, the bottom right corresponds to δ CP ¼ −π=2, top left to δ CP ¼ π=2, and the middle to δ CP ¼ 0; AEπ.
Credible interval contours (68% and 90%) are produced by a Bayesian Markov chain Monte Carlo for the standard,  In the biprobability fit with β fixed to 1, two lobes appear in the contours, which correspond to the two mass orderings: the upper lobe to the inverted orderings, and the lower to the normal ordering. These lobes coincide with the maximally CP-violating δ CP value regions of the two T2K expectation ovals, shown in brown (normal ordering) and orange (inverse ordering). The width of the credible intervals comes mainly from the uncertainties in sin 2 ð2θ 13 Þ and sin 2 ðθ 23 Þ, and height from δ CP and the mass ordering. This effect disappears in the bievent space after including statistical fluctuations in the contours for easier comparison against the data point.
The free β fit explores a larger area, especially in Pðν μ →ν e Þ andν e , which is expected; the lower number ofν e than ν e candidate events leads to a higher uncertainty in Pðν μ →ν e Þ, when not constrained by the PMNS model; additionally, the two probabilities are now decoupled due to the additional β parameter, giving independent results for both probabilities and both event rates. These credible intervals can be used to compare other neutrino oscillation models against the fit constrained by the PMNS model and against the free β fit that represents the information given by the T2K data with additional freedom.
The 90% and the 68% credible intervals from both continuous-β and PMNS-constrained fits significantly overlap. There is good agreement between the two fits, showing consistency between T2K data and the PMNS model. Additionally, the value of β is consistent with 1 (90% credible interval [0.3,1.06]), when marginalizing over all other oscillation parameters. The data point is well within the 68% credible interval in both fits after including the statistical fluctuations.
Conclusions.-The T2K Collaboration has searched for ν e appearance in aν μ beam using a dataset twice as large as in its previous searches. The data have been analyzed within two frameworks, and have been compared to predictions with either noν e appearance orν e appearance as expected from the PMNS model prediction. In both frameworks, the data are consistent with the presence ofν e appearance and no significant deviation from the PMNS prediction is seen. Using full rate and shape information, the no-appearance scenario is disfavored with a significance of 2.40 standard deviations. * Deceased. † Also at J-PARC, Tokai, Japan. ‡ Affiliated member at Kavli IPMU (WPI), the University of Tokyo, Japan. § Also at National Research Nuclear University "MEPhI" and Moscow Institute of Physics and Technology, Moscow, Russia. ∥ Also at the Graduate University of Science and Technology, Vietnam Academy of Science and Technology. ¶ Also at JINR, Dubna, Russia.    3. Biprobability (left) and bievent (right) credible interval comparison between the standard fit constrained by the PMNS (light blue) model and the non-PMNS fit with the free β parametrization (dark blue). The maximum posterior density point is marked as the 2D mode. The narrow T2K prediction ovals for normal and inverse mass orderings are in brown and orange, respectively. In the ellipses, the bottom right corresponds to δ CP ¼ −π=2, top left to δ CP ¼ π=2, and the middle to δ CP ¼ 0, AEπ. All probabilities are calculated at 600 MeV. The bievent credible intervals include statistical Poisson fluctuation.