Measurement of neutrino and antineutrino neutral-current quasielastic-like interactions on oxygen by detecting nuclear de-excitation γ-rays

K. Abe, R. Akutsu, A. Ali, C. Alt, C. Andreopoulos, 34 L. Anthony, M. Antonova, S. Aoki, A. Ariga, Y. Ashida, E.T. Atkin, Y. Awataguchi, S. Ban, M. Barbi, G.J. Barker, G. Barr, C. Barry, M. Batkiewicz-Kwasniak, A. Beloshapkin, F. Bench, V. Berardi, S. Berkman, 61 L. Berns, S. Bhadra, S. Bienstock, A. Blondel, ∗ S. Bolognesi, B. Bourguille, S.B. Boyd, D. Brailsford, A. Bravar, D. Bravo Berguño, C. Bronner, A. Bubak, M. Buizza Avanzini, J. Calcutt, T. Campbell, S. Cao, S.L. Cartwright, M.G. Catanesi, A. Cervera, A. Chappell, C. Checchia, D. Cherdack, N. Chikuma, G. Christodoulou, J. Coleman, G. Collazuol, L. Cook, 28 D. Coplowe, A. Cudd, A. Dabrowska, G. De Rosa, T. Dealtry, P.F. Denner, S.R. Dennis, C. Densham, F. Di Lodovico, N. Dokania, S. Dolan, O. Drapier, J. Dumarchez, P. Dunne, L. Eklund, S. Emery-Schrenk, A. Ereditato, P. Fernandez, T. Feusels, 61 A.J. Finch, G.A. Fiorentini, G. Fiorillo, C. Francois, M. Friend, † Y. Fujii, † R. Fujita, D. Fukuda, R. Fukuda, Y. Fukuda, K. Gameil, 61 C. Giganti, T. Golan, M. Gonin, A. Gorin, M. Guigue, D.R. Hadley, J.T. Haigh, P. Hamacher-Baumann, M. Hartz, 28 T. Hasegawa, † N.C. Hastings, T. Hayashino, Y. Hayato, 28 A. Hiramoto, M. Hogan, J. Holeczek, N.T. Hong Van, 27 F. Iacob, A.K. Ichikawa, M. Ikeda, T. Ishida, † T. Ishii, † M. Ishitsuka, K. Iwamoto, A. Izmaylov, 26 B. Jamieson, S.J. Jenkins, C. Jesús-Valls, M. Jiang, S. Johnson, P. Jonsson, C.K. Jung, ‡ M. Kabirnezhad, A.C. Kaboth, 53 T. Kajita, ‡ H. Kakuno, J. Kameda, D. Karlen, 61 S.P. Kasetti, Y. Kataoka, T. Katori, Y. Kato, E. Kearns, 28, ‡ M. Khabibullin, A. Khotjantsev, T. Kikawa, H. Kim, J. Kim, 61 S. King, J. Kisiel, A. Knight, A. Knox, T. Kobayashi, † L. Koch, T. Koga, A. Konaka, L.L. Kormos, Y. Koshio, ‡ K. Kowalik, H. Kubo, Y. Kudenko, § N. Kukita, S. Kuribayashi, R. Kurjata, T. Kutter, M. Kuze, L. Labarga, J. Lagoda, M. Lamoureux, M. Laveder, M. Lawe, M. Licciardi, T. Lindner, R.P. Litchfield, S.L. Liu, X. Li, A. Longhin, L. Ludovici, X. Lu, T. Lux, L.N. Machado, L. Magaletti, K. Mahn, M. Malek, S. Manly, L. Maret, A.D. Marino, J.F. Martin, T. Maruyama, † T. Matsubara, K. Matsushita, V. Matveev, K. Mavrokoridis, E. Mazzucato, M. McCarthy, N. McCauley, K.S. McFarland, C. McGrew, A. Mefodiev, C. Metelko, M. Mezzetto, A. Minamino, O. Mineev, S. Mine, M. Miura, ‡ L. Molina Bueno, S. Moriyama, ‡ J. Morrison, Th.A. Mueller, L. Munteanu, S. Murphy, Y. Nagai, T. Nakadaira, † M. Nakahata, 28 Y. Nakajima, A. Nakamura, K.G. Nakamura, K. Nakamura, 16, † S. Nakayama, 28 T. Nakaya, 28 K. Nakayoshi, † C. Nantais, T.V. Ngoc, ¶ K. Niewczas, K. Nishikawa, ‖ Y. Nishimura, T.S. Nonnenmacher, F. Nova, P. Novella, J. Nowak, J.C. Nugent, H.M. O’Keeffe, L. O’Sullivan, T. Odagawa, K. Okumura, 28 T. Okusawa, S.M. Oser, 61 R.A. Owen, Y. Oyama, † V. Palladino, J.L. Palomino, V. Paolone, W.C. Parker, P. Paudyal, M. Pavin, D. Payne, G.C. Penn, L. Pickering, C. Pidcott, E.S. Pinzon Guerra, C. Pistillo, B. Popov, ∗∗ K. Porwit, M. Posiadala-Zezula, A. Pritchard, B. Quilain, T. Radermacher, E. Radicioni, B. Radics, P.N. Ratoff, E. Reinherz-Aronis, C. Riccio, E. Rondio, S. Roth, A. Rubbia, A.C. Ruggeri, A. Rychter, K. Sakashita, † F. Sánchez, C.M. Schloesser, K. Scholberg, ‡ J. Schwehr, M. Scott, Y. Seiya, †† T. Sekiguchi, † H. Sekiya, 28, ‡ D. Sgalaberna, R. Shah, 42 A. Shaikhiev, F. Shaker, A. Shaykina, M. Shiozawa, 28 W. Shorrock, A. Shvartsman, A. Smirnov, M. Smy, J.T. Sobczyk, H. Sobel, 28 F.J.P. Soler, Y. Sonoda, J. Steinmann, S. Suvorov, 6 A. Suzuki, S.Y. Suzuki, † Y. Suzuki, A.A. Sztuc, M. Tada, † M. Tajima, A. Takeda, Y. Takeuchi, 28 H.K. Tanaka, ‡ H.A. Tanaka, 60 S. Tanaka, L.F. Thompson, W. Toki, C. Touramanis, K.M. Tsui, T. Tsukamoto, † M. Tzanov, Y. Uchida, W. Uno, M. Vagins, 5 S. Valder, Z. Vallari, D. Vargas, G. Vasseur, C. Vilela, W.G.S. Vinning, T. Vladisavljevic, 28 V.V. Volkov, T. Wachala, J. Walker, J.G. Walsh, Y. Wang, D. Wark, 42 M.O. Wascko, A. Weber, 42 R. Wendell, ‡ M.J. Wilking, C. Wilkinson, J.R. Wilson, R.J. Wilson, K. Wood, C. Wret, Y. Yamada, ‖ K. Yamamoto, †† C. Yanagisawa, ‡‡ G. Yang, T. Yano, K. Yasutome, S. Yen, N. Yershov, M. Yokoyama, ‡ T. Yoshida, M. Yu, A. Zalewska, J. Zalipska, K. Zaremba, G. Zarnecki, M. Ziembicki, E.D. Zimmerman, M. Zito, S. Zsoldos, and A. Zykova

Neutrino-and antineutrino-oxygen neutral-current quasielastic-like interactions are measured at Super-Kamiokande using nuclear de-excitation γ-rays to identify signal-like interactions in data from a 14.94 (16.35) × 10 20 protons-on-target exposure of the T2K neutrino (antineutrino) beam. The measured flux-averaged cross sections on oxygen nuclei are ⟨σν-NCQE⟩ = 1.70 ± 0.17(stat.) +0. 51 −0.38 (syst.) × 10 −38 cm 2 /oxygen with a flux-averaged energy of 0.82 GeV and ⟨σν-NCQE⟩ = 0.98 ± 0.16(stat.) +0. 26 −0.19 (syst.) × 10 −38 cm 2 /oxygen with a flux-averaged energy of 0.68 GeV, for neutrinos and antineutrinos, respectively. These results are the most precise to date, and the antineutrino result is the first cross section measurement of this channel. They are compared with various theoretical predictions. The impact on evaluation of backgrounds to searches for supernova relic neutrinos at present and future water Cherenkov detectors is also discussed.

I. INTRODUCTION
Measurements of neutrino neutral-current (NC) processes give insight into neutrino-nucleus interactions and are important for understanding the nucleon itself as well as improving the sensitivity of searches for a variety of physics phenomena. The strange quark content of the nucleon (∆s), for instance, can be probed via NC interactions (see Ref. [1] and references therein), and its measurements have been demonstrated by the BNL E734 experiment [2] and the MiniBooNE experiment [3,4]. Precision measurements of the neutrino-and antineutrino-oxygen NC interactions in the sub-GeV region, where the quasielastic process is expected to be dominant, also benefit a diverse array of searches with water Cherenkov detectors, such as Super-Kamiokande (SK) [5], its future upgrade, SK-Gd [6], and its successor, Hyper-Kamiokande [7]. In supernova relic neutrino (SRN) searches [8][9][10], the present uncertainty on these interactions induces a large error on atmospheric neutrino backgrounds, limiting the sensitivity at low energies where the SRN flux is predicted to be large. When searching for dark matter in accelerator neutrino experiments, as suggested in Refs. [11,12], the rate of NC interactions must be accurately estimated as they are indistinguishable from the signal. Another motivation arises in * now at CERN † 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 ∥ deceased * * also at JINR, Dubna, Russia † † also at Nambu Yoichiro Institute of Theoretical and Experimental Physics (NITEP) ‡ ‡ also at BMCC/CUNY, Science Department, New York, New York, U.S.A.
the search for sterile neutrinos in accelerator neutrino experiments [13][14][15]. The fact that the NC interaction cross section does not depend on the neutrino flavor makes it possible to search for a deficit of NC events, which would be interpreted as transitions from active to sterile neutrinos. NC interactions at the neutrino energies of interest here (E ν ≲ 1 GeV) are difficult to observe in water Cherenkov detectors because their final state particles are either neutral or charged but often below the Cherenkov threshold. Instead, the present work seeks to identify these interactions using Cherenkov light arising from the electromagnetic cascade produced by γ-rays emitted from the de-excitation of the recoil nucleus [16][17][18][19]. At E ν ≳ 200 MeV, the NC quasielastic nucleon knock-out (NCQE) processes, become dominant over NC inelastic processes without nucleon knock-out, ν(ν) + 16 O → ν(ν) + 16 O * [19]. The resulting excited nuclei relax to the ground state with the emission of γ-rays promptly. These γ-rays are available as a probe to study the NCQE interaction as has been demonstrated at T2K [20] and SK [21]. Previous studies at T2K measured the neutrino-oxygen NCQE interaction cross section with a data set of 3.01 × 10 20 protons-ontarget (POT) and SK measured this process with its atmospheric neutrino data, which is a mixture of neutrino and antineutrino interactions. Both measurements suffer from large statistical and systematic uncertainties. This paper reports the updated result from T2K using neutrinos and the first measurement using antineutrinos. In this work the signal is termed "NCQE-like", to highlight the fact that the event selection may contain contributions from NC two-particle-two-hole (2p2h) interactions where two nucleons are involved in the interaction via meson-exchange currents. Previous studies [20,21] may have also included such events, though they were not addressed specifically. Further descriptions will be given in Section VII. In the analysis, data taken with exposures of 14.94 × 10 20 POT in neutrino mode and 16.35 × 10 20 POT in antineutrino mode are used. Both the statistical and systematic errors have been reduced with the present analysis.
The paper is structured as follows. First, Section II details the experimental setup of T2K. Section III explains the Monte Carlo (MC) simulation and is followed by descriptions of the event reconstruction and selection in Section IV. Estimates of uncertainties in the analysis are described in Section V before cross section results are given in Section VI. After discussion of the results in Section VII concluding remarks are given in Section VIII.

II. THE T2K EXPERIMENT
The T2K experiment [22] has been designed for precise measurement of neutrino oscillation parameters [23] and has a broad program of additional physics measurements. It consists of the J-PARC neutrino beamline, near detectors, and SK as its far detector. T2K has taken data in nine separate run periods, termed Runs 1−9, and its beam intensity has increased throughout. Protons are bundled into eight bunches (six in Run 1), referred to as a spill, and accelerated to 30 GeV/c by the J-PARC Main Ring synchrotron. Bunches are approximately 100 ns wide and separated by about 580 ns and spills are delivered to the neutrino production target with a repetition rate of 2.48 s. Hadrons produced in proton-target (graphite) interactions are efficiently focused and sign-selected by magnetic fields produced by three electromagnetic horns [24,25], before entering a decay volume. The polarity of the magnetic horns can be changed, allowing selection and focusing of either positively or negatively charged hadrons to produce beams composed of predominantly neutrinos or antineutrinos following the decay of the hadrons. The former is referred to as forward horn current (FHC) mode while the latter is referred to as reverse horn current (RHC) mode. Located 280 m away from the graphite target the two near detectors, INGRID [26] and ND280 [27,28], are placed on-axis and 2.5 • off-axis with respect to the proton beam direction, respectively. ND280 is used to measure the (anti)neutrino spectrum before the onset of neutrino oscillations and INGRID monitors the (anti)neutrino beam direction and intensity to ensure beam quality during data taking. In addition to the INGRID measurements a muon monitor placed just after the decay volume measures the beam direction and intensity on a bunch-bybunch basis by detecting muons from pion and kaon decays [29][30][31].
Super-Kamiokande is located 295 km away from the target and 2.5 • off-axis. Beam timing information is shared between J-PARC and SK via a GPS system. It is a cylindrical water Cherenkov detector located 1,000 m under Mt. Ikeno in Kamioka, Japan. The detector is divided into two parts, an inner detector (ID) and an outer detector (OD). The ID measures 33.8 m in diameter and 36.2 m in height and is instrumented with 11,129 20inch inward-facing photomultiplier tubes (PMTs) on its wall, while the entire detector volume, which includes the ∼2 m thick OD region, extends 2.75 m radially and 2.6 m above and below the ID. Serving primarily as a veto, the OD is equipped with 1,885 8-inch outward-facing PMTs attached on the back side of the ID wall. The entire volume is filled with 50 kton of ultra-pure water. In the present work, data from the fourth stage of the detector, known as SK-IV, are used. Further descriptions of SK can be found in Ref. [5].

III. EVENT SIMULATION
Simulation of the signal and background processes are essential to the optimization of the event selection and determination of systematic uncertainties in this analysis. Monte Carlo (MC) events generated according to models of neutrino beam, neutrino interactions, and the detector response including the γ-ray emission are considered.

A. Neutrino flux
The neutrino flux is estimated by simulation based on FLUKA2011 [32] and GEANT3 [33] for modeling hadronic interactions and particle transport and decays in the beamline. Pion and kaon production cross sections are renormalized using data from the NA61/SHINE experiment taken using both thin and T2K replica targets [34][35][36][37][38]. Oscillations are taken into account for neutrinos that produce charged-current (CC) interactions at SK, using parameters from the recent T2K measurements [23]. Figure 1 shows the predicted T2K fluxes in the FHC and RHC modes without neutrino oscillations.

B. Neutrino interaction
NEUT (version 5.3.3) [39] is used to simulate neutrinonucleon interactions and subsequent final state interactions inside the target nucleus. For NCQE interactions the nominal nucleon momentum distribution is based on the Benhar spectral function [19,40], while for CC quasielastic (CCQE) interactions the relativistic Fermi gas model [41] is used. The axial-vector mass is M QE A = 1.21 GeV/c 2 and the Fermi momentum for oxygen is 225 MeV/c. CC 2p2h interactions are modeled with the calculation in Ref. [42], but their neutral counterpart is not implemented in NEUT since no model is available in the literature. The simulation uses BBBA05 vector form factors [43] and a dipole axial-vector form factor. Single pion production is based on the model of Rein and Sehgal [44]. The axial-vector mass in the resonance interac- tion is M RES A = 0.95 GeV/c 2 . Deep inelastic scattering is simulated using the GRV98 parton distribution [45] with corrections by Bodek and Yang [46]. The final state interactions of hadrons inside the nucleus are simulated with a cascade model as described in Refs. [39,47]. Further simulation details are given in Ref. [47].

C. γ-ray emission and detector response
The emission of γ-rays from nuclear de-excitation is a key part of this analysis and is simulated separately for those produced by the neutrino-nucleus interactions (primary-γ) and those from nucleon-nucleus interactions (secondary-γ). These processes are schematically illustrated in Figure 2.
After the initial neutrino interaction an excited state of the remaining nucleus is selected based on the probabilities calculated in Ref. [19]. There are four possible states, (p 1/2 ) −1 , (p 3/2 ) −1 , (s 1/2 ) −1 , and others. Here (state) −1 represents the state of the nucleus after a nucleon that initially occupied states = p 1/2 , p 3/2 , s 1/2 is removed from the nucleus. The probability for each of four states to be produced is 0.158, 0.3515, 0.1055, and 0.385, respectively [19]. The (p 1/2 ) −1 state is the ground state of 15 O or 15 N and therefore leads to no γ-ray emission. Conversely, (p 3/2 ) −1 almost always emits one γray, with 6.18 MeV from 15 O and 6.32 MeV from 15 N being the most likely. Since (s 1/2 ) −1 is a higher excited state, the branching fraction to decays including nucleons or alpha particles may be large. After such decays, the resulting nuclei may decay with γ-ray emission if it is still in an excited state thereafter. The others state includes all other possibilities and mainly includes contributions from short-range correlations among nucleons. At present there is no data nor theoretical predictions of γ-ray emission for the states covered by others so in the nominal simulation they are integrated into (s 1/2 ) −1 . A systematic uncertainty stemming from this choice is described in Section V. Further detailed descriptions on the treatment of these states are given in Ref. [20]. The interactions of secondary particles inside SK and the response of its PMTs are simulated with a GEANT3based package [33]. Hadronic interactions are of particular importance to the present analysis, especially models of neutron-nucleus reactions and the resulting γ-ray emission. These are handled by GCALOR [48,49], which implements the MICAP model for neutrons below 20 MeV and NMTC above 20 MeV. The MICAP model uses experimental cross sections from the ENDF/B-V library [50], while NMTC is based on an intra-nuclear cascade model.

IV. RECONSTRUCTION AND SELECTION
Each event in SK is reconstructed with tools used for solar neutrino analysis [51][52][53]. The visible energy (E rec ) is reconstructed using the number of hit PMTs. At these energies PMTs usually have registered only one photoelectron and there are typically between 10 and 200 hit PMTs in the current analysis window. Note that the defi-nition of energy in the present work differs from the previous T2K work [20], where the electron mass (0.511 MeV) was added to the visible energy. The current definition is consistent with recent low energy analyses in SK [21,53]. The interaction vertex and direction are inferred from the PMT hit pattern and timing. A Cherenkov angle (θ C ) for each event is calculated as the most frequently occuring value in the distribution of opening angles to all three-hit combinations of PMTs. Various calibrations are used to evaluate the performance of the reconstruction as detailed in Refs. [54,55].
This analysis considers five event categories, neutrino NCQE interactions ("ν-NCQE"), antineutrino NCQE interactions ("ν-NCQE"), all other NC interactions ("NC-other"), CC interactions, and accidental (beamunrelated) backgrounds. Both the NC-other and CC categories include contributions from neutrinos and antineutrinos. Note that these event categories reflect the neutrino interaction prior to additional particle interactions within the nucleus. This means that, for example, the NC-other sample contains pion production events where a pion was produced but was later absorbed in the nucleus. The first four interactions are simulated using NEUT and beam-unrelated backgrounds are estimated using data outside of the T2K spill timing window. Event selection criteria are tuned to effectively select signal events, ν-NCQE andν-NCQE interactions, while removing other events as follows.
(1) Events are required to be in the energy range 3.49 to 29.49 MeV, above which CC interactions become dominant. Only data judged to be of good quality, based on the beam and detector conditions during each spill, are used [47]. To select beam-induced events with high purity, the reconstructed event timing is required to be within ±100 ns of the expected timing of each bunch ("on-timing"). A sample of beam-unrelated events is selected by applying the same energy and quality cuts in a time window [−500, −5] µs before the beam spill ("offtiming"). Events with hit clusters in a window spanning 20 to 0.2 µs before the event trigger which are consistent with activity from electrons produced in the muon or pion decay chain (decay-e's) are removed. The effect on the signal efficiency by this cut is negligible.
(2) Several additional event selection cuts are applied to remove backgrounds from radioactive impurities from the detector walls. First, a fiducial volume (FV) cut is applied to all events, which requires the distance between the reconstructed vertex position and the ID wall (dwall) to be more than 200 cm. Below 6 MeV radioactive backgrounds increase considerably, requiring tighter dwall and reconstructed event quality cuts. Cuts in this energy region are tuned (discussed below) using three variables, dwall, effwall, and ovaQ. Here effwall is the distance from the event vertex to the ID wall as measured backward along the reconstructed track direction. The ovaQ parameter is a measure of the reconstruction quality and is defined as the difference of two parameters, ovaQ = g 2 vtx −g 2 dir , where g vtx and g dir are the vertex and direction fit quality parameters, respectively [56]. Cuts on these parameters are optimized for five regions between 3.49 and 5.99 MeV with each 0.5 MeV bin width. The optimization is performed separately for each T2K run period because the detector condition and the beam power differ from run to run. A figure-of-merit (FOM) designed to maximize sensitivity to the NCQE signal is defined as: where N sig is the number of signal events predicted by the MC (ν-NCQE for FHC andν-NCQE for RHC) and N bkg is the total number of background events. The latter is composed of two components, N MC bkg and N beam-unrelated bkg , which represent non-signal neutrino events such as NCother and CC interactions, and beam-unrelated events from the off-timing data sample, respectively. Cuts on the three parameters above are chosen to maximize the FOM in each energy region. As an illustration the optimized values of dwall, effwall, and ovaQ for one of the FHC mode runs (T2K Run 8) are shown in Figure 3. A linear function is fit to each distribution to obtain the final cut criteria and is denoted by the red line in the figure.
For the dwall and effwall distributions, if the optimized value is 200 cm (the FV cut criterion) in two successive energy bins, the second and later bins are removed and the fit is repeated. In the end, each of these three parameters is required to be larger than the obtained line. That is, events with values in the upper right portion of the plots in the figure are kept. Note that at higher energies the optimum dwall and effwall values fall below 200 cm, but such events are already removed by the initial FV cut. Figure 4 shows the ovaQ distributions after the cuts described in (1), the FV cut, the optimized dwall cut, and the optmized effwall cut. There is clear separation between signal and background. Further description of the variables used in this selection are given in Refs. [20,56].
(3) The final phase of the event selection is focused on the removal of CC interaction events. A single charged particle whose momentum is large compared to its mass is likely to have a Cherenkov angle of ∼42 • in water. On the other hand if the particle momentum is lower, the reconstructed Cherenkov angle decreases. In this analysis low energy muons from CC interactions and still above Cherenkov threshold distribute around θ C = 20 • −35 • , whereas decay-e's have θ C ∼ 42 • . The contribution of each can be seen in Figure 5. To reduce these CC events, a linear cut in the reconstructed energy and Cherenkov angle plane is chosen by maximizing the FOM defined in Eq. (3). In the figure the resulting cut is shown with a red line. This is performed separately for the FHC and RHC samples. Using the optimized cut the signal efficiency is 99% (99%) while 63% (58%) of CC events are removed in FHC (RHC) mode. Some CC-other events still remain after this cut, which could be due to, for example, multiple-γ emission via neutron production (as explained later), but this fraction is small with respect to the total number of selected events. Similar population is seen also in the NC-other distribution. After all cuts, the event selection is more than 80% efficient for signal events, while reducing background events by more than two orders of magnitude. Figure 6 shows a comparison of the number of MC beam neutrino events against beam-unrelated events both before and after these cuts. The event selection summary for the beam data and MC is shown in Table I. After the event selection, 204 events are observed in the FHC data and 97 events are observed in the RHC data. These are compared with the number of predicted events in Table I. While the FHC sample has a high signal purity, the neutrino component forms nearly 20% of the RHC sample because of the difference between the neutrino and antineutrino cross sections. Figures 7 and 8 show distributions of the reconstructed energy, Cherenkov angle, and vertex position for the FHC and RHC samples, respectively. The observed E rec distributions agree well with the predictions in both FHC and RHC modes. A clear contribution from ∼6 MeV γ-rays is observed in both operation modes. In the FHC θ C distribution, the data at high angles is below the MC expectation, while no such MC excess is seen in the RHC data. This excess was also observed in the previous T2K measurement [20] although the statistical error was larger. At high angles this distribution is dominated by events with multiple γ-rays. Such events are caused mainly by fast neutron interactions with nuclei in the water. The excess in FHC may then be attributed to inaccurate modeling of secondary neutron reactions and their subsequent γ-ray emissions. The fact that the disagreement between observation and prediction is visible in FHC and not in RHC, may be understood by the difference in the out-going nucleon kinematics between neutrino and antineutrino interactions. Helicity conservation in antineutrino interactions produces more forward-going leptons in the final state and consequently lower momentum nucleons. The latter therefore goes on to produce fewer secondary γ-rays than that from its neutrino interaction counterpart. Comparing the ratio of the single-γ peak (∼42 • ) to the multiple-γ peak (∼90 • ) of the MC in each figure, there are relatively fewer events in the high-angle region of the RHC sample. The vertex positions of selected events in the data are found to be uniform and no bias relative to the beam direction is observed.

V. UNCERTAINTY ESTIMATES
Based on the observed number of events in Table I, the associated statistical error is 7.0% for the FHC sample and 10.2% for the RHC sample.
Systematic errors from six main sources are considered in the analysis, namely the neutrino flux prediction, the neutrino interaction model, the primary-γ and secondary-γ emission models, neutrino oscillation parameters, and the detector response. In this analysis, CC measurement results from the T2K near detectors are not used so as to ensure that flux and interaction systematics are treated independently. Only statistical uncertainties are considered for beam-unrelated events, 3.0% in the FHC sample and 3.9% in the RHC sample, since they are also part of the observed data and respond to detector uncertainties in the same way. The effect of possible rate fluctuations between the on-and off-timing windows is negligible. Table II summarizes the impact of each of these error categories on the different interaction modes populating the samples. Among them, systematic errors from the secondary-γ production model are the leading uncertainties. The error sources are described in detail below.

A. Neutrino flux and interaction model uncertainties
The impact of neutrino flux and interaction systematic uncertainties in this analysis is estimated by the change in the number of selected events relative to the nominal model under a 1σ shift in each error source. The procedure follows previous T2K analyses [20,57,58].
Flux uncertainties are evaluated for each neutrino flavor, horn polarity, and neutrino energy bin. Uncertainties in the hadronic interaction cross section are the dominant contribution to the assigned 6−8% flux uncertainties. This represents a large improvement over previous T2K analyses, due to improved hadron production and interaction constraints from NA61/SHINE measure- ments using a replica of the T2K target [38]. The value of the axial-vector mass used to generate quasielastic interactions with its 1σ error is M QE A = 1.21 ± 0.18 GeV/c 2 . Similarly the Fermi momentum in oxygen is taken to be 225 ± 31 MeV/c. Parameters describing contributions from 2p2h interactions, resonant pion production, and deep inelastic scattering follow the assignments in previous analyses [20,57,58]. These result in uncertainties of 8.2% (10.8%) for the NC-other and 16.5% (38.2%) for CC interaction backgrounds in the FHC (RHC) measurement. The larger uncertainty in the RHC CC component, as seen in Table II, is attributed to the different effect of M QE A . Since γ-rays are emit-ted isotropically and SK has 4π acceptance, the signal efficiencies are unaffected by neutrino interaction model uncertainties.
It should be noted that while NC inelastic scattering without nucleon emission, ν(ν) + 16 O → ν(ν) + 16 O * , should be present in the selected sample, it is not simulated in this analysis. According to Ref. [59], the sum of cross sections leading to 15 O * and 15 N * after the 16 O(ν, ν ′ ) interaction increases from 6.7 × 10 −42 cm 2 at E ν = 50 MeV to 481 × 10 −42 cm 2 at E ν = 500 MeV, while it is almost constant above ∼200 MeV. By comparing this to the expected NCQE cross section in Ref. [19], it is found that the NCQE process dominates over the NC inelastic process without nucleon knock-out above E ν ∼ 200 MeV. In addition, the former cross section is ∼40 times larger at 500 MeV and is expected to be even larger at higher energies. In the present measurement the signal is predominantly from neutrinos above E ν ∼ 500 MeV. Assuming that the detection efficiency of γ-rays produced from the de-excitation of nuclei recoiling from the NC inelastic interaction without nucleon emission is comparable to that of NCQE scattering, a 3% error on the signal channel is assigned conservatively in consideration of the expected interaction cross section differences. Another possible contribution is from NC interactions on hydrogen, ν(ν)+ 1 H → ν(ν)+ 1 H, where the final state protons may produce γ-rays via reactions with water. However, the contribution from such interactions is expected to be less than 1% of that from NCQE interactions on oxygen and therefore does not significantly affect the results of the present measurement.

B. Primary-and secondary-γ production uncertainties
Errors on the primary γ-ray emission come from the uncertainties on the spectroscopic factors. Calculation of the spectroscopic strength for the p 3/2 state has been found to be consistent with electron scattering data within 5.4% [19], which leads to an error on the observed event rate at T2K of less than 3%. The uncertainty due to the others state (all other states than (p 1/2 ) −1 , (p 3/2 ) −1 , and (s 1/2 ) −1 ) being included into the (s 1/2 ) −1 state in the nominal model is estimated by comparison with an extreme case. Since no significant deviation in the predicted p 3/2 strength has been observed in (e, e ′ p) and (p, 2p) experiments [60,61], others cannot behave like the (p 3/2 ) −1 state. In contrast, the possibility that the others state behaves like the ground state, (p 1/2 ) −1 , emitting no γ-rays, is considered, since this would not contradict any existing data. To model this, the others state is included into (p 1/2 ) −1 instead, and the change in the event rate relative to the nominal model is taken as the systematic error. This results in uncertainties in the 6−12% range for the signal and background modes. This extreme case covers the uncertainties of the p 1/2 and s 1/2 spectroscopic strengths. The total error on primaryγ production is taken to be the sum in quadrature of above two sources.
The secondary-γ emission rate is model-dependent and at present there is insufficient data on γ-ray emission from neutron-oxygen reactions at energies above 20 MeV [62], which are most relevant for the present work, making model selection difficult. Since different models predict different amounts of γ-ray emission, to reduce the impact of such model dependence, instead the total number of emitted Cherenkov photons from secondary emission processes is considered. First, the probability (P selected ) of an event being reconstructed in the 3.49−29.49 MeV energy region of this analysis is estimated as a function of the number of emitted Cherenkov photons using MC. The resulting probabilities for FHC and RHC are shown in Figure 9. The number of emitted Cherenkov photons (N C ) can be broken down into three parts, Here N primary-ν C denotes the contribution from the primary γ-ray emission and N secondary-n C (N secondary-p C ) is from secondary γ-rays produced by neutron (proton) interactions in water. The systematic uncertainty used in the analysis is estimated by varying the contributions from these secondary interactions and calculating the change in the selected sample using Figure 9. The source of uncertainty can be broken down into the initial nucleon-oxygen interaction and the subsequent nuclear de-excitation. In Ref. [63], proton-carbon data were fit  to obtain a constraint on the nucleon-nucleus scattering cross section. Their result showed a 30% difference between the measured and predicted (GCALOR) cross sections. In the present work, the target nucleus is different but the effect is found to be no larger than 5% in neutrino interaction measurements [58], so a conservative error of 40% is adopted. In order to estimate the impact of γ-ray emission from fast neutron reactions on oxygen, the results of a muon-induced spallation study at SK [64] are used. Since the selected sample contains contributions from such neutron interactions, and the measured energy distribution does not differ by more than 50% from the MC, this number is taken as the error estimate. For the uncertainty propagation the quadratic sum of these two contributions is used and a ±65% variation is applied to both N secondary-n C and N secondary-p C . The variation producing the largest change in the final sample is used to compute the final error and results in a ∼13% uncertainty for signal and roughly 20% for the NC-other and CC components. In addition, the impact of uncertainties from the final state interaction model has been evaluated to be as large as 3%. The total uncertainty for each is obtained by summing these two contributions in quadrature.

C. Oscillation parameter and detector response uncertainties
Errors on the oscillation parameters, θ 13 , θ 23 , and ∆m 2 32 , are taken from Ref. [23]. Varying each of these, the change in the selected number of CC events results in 3−4% errors for the FHC and RHC samples.
Errors on each reconstructed parameter used in the event selection, E rec , dwall, effwall, ovaQ, and θ C are considered as detector response uncertainties. These have been studied using detector calibrations [54,55], and their effect on the final sample is 1%. Similarly, the gain of the SK PMTs was found to vary over the observation period and its impact is considered as systematic error in this analysis. This gain shift changes the number of PMT hits used to reconstruct energy and produces a 3% error on the final sample. In total, 3−5% errors are assigned for each interaction mode.

VI. CROSS SECTION RESULTS
The number of observed events in the FHC and RHC data (D FHC obs and D RHC obs , respectively) are expressed as follows:  Tables I and  II. Correlations among the flux and cross section parameters are not considered in this analysis. The systematic uncertainty on primary-γ production is considered to be fully correlated among the different interaction types and operation modes, and the secondary-γ production error is treated in the same way, since the change of the γray emission rate should be common for the neutrino interaction types and T2K operation modes. Note that the primary-and secondary-γ production uncertainties are uncorrelated. Distributions of the calculated scale factors for one million pseudo experiments are shown in Figures 10 and 11. Here the dominant error is the secondary γ-ray model uncertainty as shown in Table II. The factors f ν-NCQE and fν -NCQE have a weak negative correlation for variations of the statistical uncertainty but a strong positive correlation under the influence of systematic uncertainties. In the end, the scale factors are measured as: The predictions of flux-averaged cross sections by NEUT for neutrino and antineutrino NCQE interactions on oxygen, ⟨σ NEUT ν-NCQE ⟩ and ⟨σ NEUT ν-NCQE ⟩, are calculated as: The nominal flux, ϕ ν = ϕ FHC ν is used for neutrinos and ϕν = ϕ RHC ν is used for antineutrinos in calculations of the flux-averaged NCQE cross sections. Note that summation is done over muon and electron (anti)neutrinos in Figure 1, though the actual flux at SK contains tau (anti)neutrinos due to neutrino oscillations. This treatment is justified because the NC cross section is flavorindependent. Here the integrations are conducted up to 10 GeV as higher energies have a negligible impact on the result. The measured flux-averaged NCQE-like cross sections on oxygen nuclei are obtained by multiplying the scale factors to each of Eqs. (8) and (9), These measurements are shown together with the predictions from NEUT in Figure 12. The neutrino measurement improves over the previous T2K result with FHC data, ⟨σ ν-NCQE ⟩ = 1.55 +0.71 −0.35 (stat. ⊕ syst.) × 10 −38 cm 2 /oxygen [20]. Covariance matrices of the neutrino and antineutrino flux-averaged NCQE-like cross sections are shown for both variations of the statistical and systematic uncertainties in Table III. Currently, there are no models available in the literature for the NC 2p2h interaction, so this channel is not simulated in the present analysis. Since NC 2p2h interactions involve multi-nucleon knock-out, not only multiple γ-rays are expected but additional secondary γ-rays from the recoil nucleons are expected as well. It should be noted that if this process exists then the selection in this analysis likely includes such events. However, if the ratio of the NC 2p2h and QE cross sections is similar to the corresponding CC ratio, roughly 5−10% [42], the present measurement will not be sensitive to these events.

B. Comparison with model predictions
The measured NCQE-like cross sections are tied to NEUT as the underlying model for signal and backgrounds. It is interesting to compare the current measurements with various theoretical models. Six models from Ref. [65] are used in the comparison: the Spectral  [40,[66][67][68][69]. The fluxaveraged NCQE cross sections for each model are compared in Figure 13. While the measured result for neutrinos is consistent with all of the models within the 1σ error, the SF, RMF, and SuSA models lie outside the 1σ region for antineutrinos. However, it is important to note that each model has its uncertainties and none of these models contains the NC 2p2h process.

C. Impact on supernova relic neutrino (SRN) searches
The present work can be used to estimate NCQE backgrounds from atmospheric neutrinos to SRN searches. Similarly, since γ-rays from NC 2p2h interactions are also a background to such searches the inclusive nature of the current measurement may provide useful constraints. Although the cross section results can be used directly, they suffer from large uncertainties from primary-and secondary-γ emission models as detailed above. If instead one uses the number of events in the expected SRN signal region, most uncertainties in Table II can be avoided and only errors arising from the difference between the T2K beam and atmospheric neutrino fluxes (<10%) and detector response error need to be considered. In the following, the present analysis sample is projected onto the E rec −θ C phase space used in the SK SRN search and divided into four regions: 1) E rec ∈ [3.49, 7.49] MeV and θ C ∈ [38,50] Figure 14 gives the E rec −θ C distributions from the FHC and RHC data and MC before the CC interaction cut and after all of the preceding cuts described in Section IV.  Figure 15. Similarly, Figure 16 shows the θ C distributions for E rec ∈ [3.49, 7.49] MeV and E rec ∈ [7.49, 29.49] MeV. Here also the FHC distributions for observation and prediction show discrepancies, which may be attributed to modeling of the secondary-γ emission. These distributions can be used to estimate the NCQE background to the SRN search by suitable weighting of the MC to data. Though beyond the scope of the present work, this is expected to significantly improve the current 100% error on this background used in the SK SRN analysis [8,9].  Events/MeV/2.7-degree FIG. 14. Two-dimensional Erec−θC distributions for FHC (top) and RHC (bottom) respectively before the CC interaction cut and after all of the preceding cuts described in Section IV. Magenta dots correspond to the observed data.

D. Future prospects
At present T2K has collected less than half of its expected POT and extensions of the experiment are being considered [70]. The larger statistics of future data sets motivate several possible improvements to the present work. Systematic errors from the secondary-γ production model can be reduced by incorporating recent measurements of γ-ray emission from neutron-oxygen interactions into MC. Measurements using 30, 80, and 250 MeV neutrons have been performed, but only results at 80 MeV are available at present [62]. Furthermore, neutron tagging at SK, particularly the high-efficiency tagging realized in the coming Gd-doped phase of Super-Kamiokande (SK-Gd), can be used to study the relationship of neutrons, their transport in water, and the production of secondary γ-rays. Information on the neutron capture vertex would further constrain the neutron kinetic energy  in NCQE interactions by measurement of the neutron flight distance from the primary interaction vertex. Neutron information would also allow for differential cross section measurements using the reconstructed Q 2 as well as studies of ∆s if proton and neutron final states can be distinguished. Finally, using the ∼8 MeV γ cascade following neutron capture on Gd, it may be possible to identify the NCQE interactions resulting in the ground state nucleus by requiring no activity by the primary-γ.

VIII. CONCLUSION
In this paper, neutrino-and antineutrino-oxygen neutral-current quasielastic-like interactions have been measured using nuclear de-excitation γ-rays at the T2K far detector, with data corresponding to 14.94×10 20 POT in FHC and 16.35 × 10 20 POT in RHC polarities. Compared to the previous T2K study, the present analysis has improved the event simulation and selection criteria, and reduced both systematic and statistical uncertainties. In addition, this work presents the first measurement of antineutrino interactions in this channel to date. The measured flux-averaged NCQE-like cross sections on oxygen nuclei are ⟨σ ν-NCQE ⟩ = 1.70 ± 0.17(stat.) +0. 51 −0.38 (syst.) × 10 −38 cm 2 /oxygen for neutrinos at a flux-averaged energy of 0.82 GeV and ⟨σν -NCQE ⟩ = 0.98 ± 0.16(stat.) +0. 26 −0.19 (syst.) × 10 −38 cm 2 /oxygen for antineutrinos at a flux-averaged energy of 0.68 GeV. Simultaneously treating both FHC and RHC data has resulted in similar sized errors for both the neutrino and antineutrino measurements. These results were found to be consistent with currently available models within the measurement precisions. In addition, MC and data comparisons in the kinematic regions of interest for SRN searches were performed. These measurements are expected to improve estimates of backgrounds to those searches not only in the present Super-Kamiokande experiment, but also in future water Cherenkov detectors such as SK-Gd and Hyper-Kamiokande. The data related to the results presented in this paper can be found in [71].
the WestGrid and SciNet consortia in Compute Canada, and GridPP in the United Kingdom. In addition, participation of individual researchers and institutions has been further supported by funds from ERC (FP7), "la Caixa" Foundation (ID 100010434, fellowship code