Measurement of the neutrino-oxygen neutral-current quasielastic cross section using atmospheric neutrinos in the SK-Gd experiment

We report the first measurement of the atmospheric neutrino-oxygen neutral-current quasielastic (NCQE) cross section in the gadolinium-loaded Super-Kamiokande (SK) water Cherenkov detector. In June 2020, SK began a new experimental phase, named SK-Gd, by loading 0.011% by mass of gadolinium into the ultrapure water of the SK detector. The introduction of gadolinium to ultrapure water has the effect of improving the neutron-tagging efficiency. Using a 552.2 day data set from August 2020 to June 2022, we measure the NCQE cross section to be 0.74 $\pm$ 0.22(stat.) $^{+0.85}_{-0.15}$ (syst.) $\times$ 10$^{-38}$ cm$^{2}$/oxygen in the energy range from 160 MeV to 10 GeV, which is consistent with the atmospheric neutrino-flux-averaged theoretical NCQE cross section and the measurement in the SK pure-water phase within the uncertainties. Furthermore, we compare the models of the nucleon-nucleus interactions in water and find that the Binary Cascade model and the Liege Intranuclear Cascade model provide a somewhat better fit to the observed data than the Bertini Cascade model. Since the atmospheric neutrino-oxygen NCQE reactions are one of the main backgrounds in the search for diffuse supernova neutrino background (DSNB), these new results will contribute to future studies - and the potential discovery - of the DSNB in SK.

We report the first measurement of the atmospheric neutrino-oxygen neutral-current quasielastic (NCQE) cross section in the gadolinium-loaded Super-Kamiokande (SK) water Cherenkov detector.In June 2020, SK began a new experimental phase, named SK-Gd, by loading 0.011% by mass of gadolinium into the ultrapure water of the SK detector.The introduction of gadolinium to ultrapure water has the effect of improving the neutron-tagging efficiency.Using a 552.2 day data set from August 2020 to June 2022, we measure the NCQE cross section to be 0.74 ± 0.22(stat.)+0.85  −0.15 (syst.)× 10 −38 cm 2 /oxygen in the energy range from 160 MeV to 10 GeV, which is consistent with the atmospheric neutrino-flux-averaged theoretical NCQE cross section and the measurement in the SK pure-water phase within the uncertainties.Furthermore, we compare the models of the nucleon-nucleus interactions in water and find that the Binary Cascade model and the Liège Intranuclear Cascade model provide a somewhat better fit to the observed data than the Bertini Cascade model.Since the atmospheric neutrino-oxygen NCQE reactions are one of the main backgrounds in the search for diffuse supernova neutrino background (DSNB), these new results will contribute to future studies -and the potential discovery -of the DSNB in SK.

I. INTRODUCTION
Neutrinos emitted from all past core-collapse supernovae comprise an integrated flux called the diffuse supernova neutrino background (DSNB) [1].Detecting rays.The gamma-rays give their energy to electrons or positrons via Compton scattering or pair production, then Cherenkov photons are emitted.By detecting the positron signal (prompt signal) and the neutron signal (delayed signal), we can remove a large number of backgrounds that do not contain neutrons.However, backgrounds that contain neutrons cannot be completely removed.
One of the main backgrounds in the DSNB search is caused by the atmospheric neutrino-oxygen neutralcurrent quasielastic (NCQE) reactions.NCQE reactions can be expressed as ν(ν) + 16 O → ν(ν) where the atmospheric neutrino knocks out a nucleon of the oxygen nucleus, and the residual nucleus may emit one or more de-excitation gamma-rays with a few MeV.When a neutron is knocked out, the combination of de-excitation gamma-rays and neutron mimics the IBD event, making it difficult to distinguish between NCQE and IBD events.Therefore, the precise estimation of NCQE events is essential for the DSNB discovery in SK-Gd.
To estimate the NCQE events precisely, the behavior of neutrons in water must be understood.In IBD events, the outgoing neutron has at most a few MeV, while in NCQE events, the knocked-out neutron may have hundreds of MeV.In the latter case, it can knock out other nucleons of oxygen nuclei in water, and additional de-excitation gamma-rays and neutrons are generated.Therefore, it is crucial to understand the nucleon-nucleus interactions in water (secondary interactions).
The T2K experiment measured the accelerator neutrino-oxygen NCQE cross section with a large uncertainty mainly coming from the de-excitation gamma-rays by secondary interactions [3,4].In previous SK DSNB searches, the expected number of NCQE background events, which was scaled using the T2K NCQE cross section, had large systematic uncertainties (60-80%) [5,6].
In the past, the secondary interaction model based on the Bertini Cascade model (BERT) [7] was the only choice in SK.However, now other secondary interaction models like the Binary Cascade model (BIC) [8] and the Liège Intranuclear Cascade model (INCL++) [9] can be employed and compared with data.In this paper, we discuss the reproducibility of the observed data in each secondary interaction model using atmospheric neutrino events.Then we report the first measurement of the atmospheric neutrino-oxygen NCQE cross section in the Gd-loaded SK water Cherenkov detector.

II. THE SUPER-KAMIOKANDE EXPERIMENT
The SK experiment [10] is located in Kamioka, Gifu, Japan, with the large water Cherenkov detector placed 1,000 m underground, resulting in a 2,700 m water equivalent overburden.The rate of cosmic ray muons is reduced by a factor of 10 5 compared to that at ground level.The detector consists of a stainless steel cylindrical water tank with a diameter of 39.3 m, a height of 42.0 m, containing 50 kilotons of ultrapure water.The water tank is optically separated into an inner detector (ID) and an outer detector (OD).The ID has 11,129 20inch PMTs to reconstruct the energy, vertex position, direction, and kind of charged particles, while the OD has 1,885 8-inch PMTs to veto incoming cosmic ray muons.Radioactive backgrounds are concentrated near the detector wall.Thus events more than 2 m away from the ID wall are used in the analyses, resulting in a fiducial volume of 22.5 kilotons.
SK started its observation in April 1996, and so far, the observation is categorized into seven phases (from SK-I to SK-VII).Since the start of SK-IV in 2008 we have been able to search for neutron signals up to 535 µs after the trigger thanks to an electronics upgrade [11,12].However, through SK-V, the neutron signal was a 2.2 MeV gamma-ray from neutron capture on free protons (hydrogen in the water), and the neutron-tagging efficiency was low.To increase the efficiency, we loaded 0.011% by mass of Gd in SK, at which time the SK-VI (SK-Gd) phase started in July 2020 [13].The time constant of neutron capture on Gd at this concentration is about 115 µs [13].Now 0.03% of Gd has been loaded in SK, and we continue the observation as SK-VII since mid-2022.
The previous atmospheric neutrino-oxygen NCQE cross section measurement [14] was performed using 2,778 days of SK-IV pure-water data from October 2008 to October 2017.This study uses 552.2 days of SK-VI data with 0.011% Gd-loaded water from August 2020 to June 2022.This data set is the same as the one used for the DSNB search in SK-VI [6].

III. SIMULATION
The atmospheric neutrino flux at the SK detector is predicted using the HKKM11 model [15], which shows good agreement with the observation in SK [16].
Neutrino interactions are simulated using NEUT [17] (version 5.4.0.1).The NCQE cross section on oxygen is based on the model using the oxygen spectral function reported by Ankowski et al. [18,19] with the BBBA05 vector form factor [20] and the dipole axial form factor [20].The state of the residual nucleus after neutrinooxygen nucleus interaction (primary interaction) is selected based on the probabilities computed in Ref. [19].There are four states, (p 1/2 ) −1 , (p 3/2 ) −1 , (s 1/2 ) −1 , and others, where (state) −1 shows the state of the nucleus after a nucleon initially occupying the state (p 1/2 , p 3/2 , or s 1/2 ) is removed.The production probability of each state is 0.1580, 0.3515, 0.1055, and 0.3850, respectively.(p 1/2 ) −1 state is the ground state of 15 O or 15 N, thus no gamma-ray is emitted.Mainly 6.18 MeV or 6.32 MeV gamma-rays are emitted from (p 3/2 ) −1 state of 15 O or 15 N, respectively [21,22].In the case of (s 1/2 ) −1 state, nucleons and gamma-rays are emitted because the excitation energy is high.The de-excitation mode is selected based on the 16 O(p, 2p) experiment [23].The others state includes all other possibilities that are not in (p 1/2 ) −1 , (p 3/2 ) −1 , and (s 1/2 ) −1 states, and there are no data nor theoretical predictions covered by this state.In our simulation, the others state is set to be integrated into (s 1/2 ) −1 state by default.
In the past, a GEANT3-based [24] SK detector simulation where only BERT was implemented for neutron tracking in water was used.However, a Geant4based [25] (version 10.05.p01)SK detector simulation has been newly developed for the SK-Gd experiment.In this simulation, BERT (FTFP_BERT_HP physics list), BIC (QGSP_BIC_HP physics list), and INCL++ (QGSP_INCLXX_HP physics list) can be used as the secondary interaction model.Here, BERT is a traditional cascade model used in GEANT.BIC uses a large set of hadron data to choose interaction processes to improve the accuracy.INCL++ is an advanced binary cascade model including phase space and quantum mechanical processes.The features of each secondary interaction model are described in Sec.V.In this NCQE cross section measurement, BERT is used as the baseline model.

IV. EVENT SELECTION
In this study, we search for NCQE events that consist of prompt signals from de-excitation gamma-rays and delayed signals from neutrons.We select the events where the visible energy of the prompt signal (E vis ) is between 7.49 MeV and 29.49MeV because there are many NCQE events in this E vis region from simulation.In each candidate event, delayed signals within 535 µs from the prompt signal are searched for.The event reconstruction and neutron-tagging method follow the DSNB search in SK-Gd phase [6,26].The event reduction is described below.
First, several backgrounds in the energy range from a few MeV to tens of MeV are removed using the same reduction code as in the solar neutrino analysis [27].The removed background events are those caused by PMT noise hits, radioactive backgrounds, and decay electrons from cosmic ray muons.These are removed by the goodness of reconstruction, the distance from the ID wall to the reconstructed vertex, and the time difference from preceding cosmic ray muons, respectively.Second, the spallation events, which are prominent backgrounds in the energy range from a few MeV to tens of MeV, are removed.These are decays of radioactive isotopes produced by nuclear spallation of oxygen nuclei induced by energetic cosmic ray muons, and are removed using the time and track of muons close to candidate events.Moreover, atmospheric neutrino events other than NCQE events are removed.For example, events associated with muons and pions are removed using two hit peaks close in FIG. 1. NCQE cumulative signal efficiencies.Event reductions are performed in the order shown in the legend, and before the reduction for spallation events is taken as 100%.Details of these event reductions are described in Ref. [5,6].
time and the cleanliness of Cherenkov rings.Finally, we select NCQE events using the reconstructed Cherenkov angle of the prompt signal (θ C ) and the number of delayed signals per event (N delayed ).Here, θ C is determined from the distribution of angles calculated by a combination of 3-PMT hits.
The difference between this study and the DSNB search in SK-VI [6] is the cut criteria for θ C and N delayed .In IBD events, only one relativistic positron and one neutron are emitted.Therefore, θ C and N delayed tend to be about 42 degrees and one, respectively.In contrast, in an NCQE event, multiple gamma-rays and multiple neutrons are easily emitted.When multiple gamma-rays are emitted, θ C tends to be larger because of the uniform distribution of the hit PMTs.Therefore, in this study, we select events for which θ C is greater than 50 degrees and N delayed is greater or equal to one.The cut criteria of θ C , E vis and N delayed are the same as the study in SK pure-water phase [14].After applying all event selections to 552.2 days of SK-VI data, 38 events remain.We confirmed that these events are uniformly distributed.The NCQE cumulative signal efficiencies are shown in Fig. 1.
This final sample includes not only NCQE events but also other events such as atmospheric neutrino NC non-QE events, atmospheric neutrino charged-current (CC) events, remaining spallation events, reactor neutrino events, and accidental coincidence events.To determine the final number of NCQE events, it is necessary to estimate the number of all those events.The expected numbers of atmospheric neutrino events in each secondary interaction model estimated by the simulation are summarized in Table I.After applying all event selections, NCQE and NC non-QE events account for about 60% and about 30% of total events, respectively.The expected number of spallation, reactor neutrino, and ac- cidental coincidence events, which are calculated by the same method as Ref. [6], are 0.9, 0.1, and 1.6, respectively.Note that the number of DSNB events predicted by the Horiuchi + 09 model [28], which is not included in the expected number of events, is 0.1.

V. SECONDARY INTERACTION MODELS
In Ref. [3,4], it was found that agreements of the secondary interaction model based on BERT remain poor and result in significant systematic uncertainty as described in Sec.I. Therefore, we compare the observed data with the other secondary interaction models using the newly developed Geant4-based SK detector simulation.Here, we use three secondary interaction models: BERT, BIC, and INCL++.Fig. 2 shows the distributions of θ C , E vis , and N delayed in each secondary interaction model.
The distributions of θ C , E vis , and N delayed strongly depend on the number of de-excitation gamma-rays and neutrons.For example, the direction of Cherenkov photons becomes more isotropic as the number of deexcitation gamma-rays is increased.Moreover, the total energy of de-excitation gamma-rays is correlated to the number of de-excitation gamma-rays.Therefore, θ C and E vis become larger as the number of de-excitation gamma-rays gets larger.Furthermore, N delayed is correlated to the number of neutrons.
In the distributions of E vis and N delayed , the number of events is larger for BERT than other two models.Furthermore, in the θ C distribution, the differences between BERT and other two models are large in high-angle regions.These differences come from the number of deexcitation gamma-rays and neutrons by secondary interactions.The number of de-excitation gamma-rays and neutrons is similar between BIC and INCL++.While, in BERT, the number of de-excitation gamma-rays and neutrons is larger than in the other two models.
We calculated the chi-square χ 2 for θ C , E vis , and N delayed distributions by using the Poisson- likelihood [29].Here, χ 2 is defined as where bin is the number of bins, N obs,i is the observed number of events of i-th bin and N exp,i is the expected number of events of i-th bin.The values are summarized in Fig. 2. Due to the small statistics, the chi-square cannot give conclusive results; however, the values are smaller for BIC and INCL++ than for BERT in all distributions.
As described in Sec.I, an accelerator neutrino-oxygen NCQE cross section measurement was conducted as part of the T2K experiment [4].The observed and expected number of events in θ C ∈ [78, 90] degrees and E vis ∈ [7.49, 29.49] MeV obtained in the T2K data analysis are shown in Table II.We have also performed an analysis using the same criteria of θ C and E vis .The results are also summarized in Table II.With these selection criteria, the expected number of events in BERT is larger than the observed number of events.The similar discrepancy was observed in T2K because the secondary interaction model is based on BERT [4,30].In both cases, the differences of the expected and observed number of events in BERT are larger than that in BIC and INCL++, which shows similar trend as above.

VI. NCQE CROSS SECTION A. Measured NCQE cross section
The flux-averaged theoretical neutrino-oxygen NCQE cross section is where φ i (E) is the atmospheric neutrino flux [15] at neutrino energy E and σ i (E) theory NCQE is the theoretical NCQE cross section [19].The integral is performed between 160 MeV and 10 GeV because the NCQE cross section is small below 160 MeV and the atmospheric neutrino flux is small above 10 GeV.The systematic uncertainty related to the energy cutoff is described in Sec.VI B. The measured neutrino-oxygen NCQE cross section is where N obs (= 38) is the observed number of events, N exp NCQE is the expected number of NCQE events, and N exp Non-NCQE is the expected number of non-NCQE events, including NC non-QE, CC, spallation, reactor neutrino, and accidental coincidence.

B. Systematic uncertainties
Systematic uncertainties on the expected NCQE, NC non-QE, and CC events are summarized in Table III.We follow the estimation methods of measurements in SK pure-water phase and T2K [3,4,14].The estimation of each systematic uncertainty is described below.
Primary interaction uncertainty arises from the spectroscopic strengths of the oxygen nucleus.Computation of the p 3/2 spectroscopic strength is consistent with 16 O(e, e ′ p) experiment within 5.4% [19,21].For the others state, there is no reliable predictions as written in Sec.III, thus the uncertainty is conservatively estimated by comparing with an extreme case, that is the difference between the default state ((s 1/2 ) −1 ) and the ground state ((p 1/2 ) −1 ).
Secondary interaction uncertainty arises from the secondary interaction model used in the detector simulation.As described in Sec.V, the chi-square differences were inconclusive.Therefore, the uncertainty is taken to be the difference in the expected number of events from BERT to BIC or INCL++.
In the calculation of Eq. ( 3), the integral is performed between 160 MeV and 10 GeV, while the expected number of atmospheric neutrino events is estimated using full energy range.Energy cutoff uncertainty is estimated by applying the energy cutoff to the expected number of atmospheric neutrino events.Since the expected number of events decreases by the energy cutoff, only negative systematic uncertainty is considered.
The uncertainty of the measured atmospheric neutrino flux in SK differs in each energy bin, as shown in Fig. 6 and Table IV of Ref. [16].In this measurement, we chose the conservative value and 18.0% in [160 MeV, 10 GeV] is applied to atmospheric neutrino flux uncertainty.
Systematic uncertainty of the measured NCQE cross section is estimated by performing toy MC considering the systematic uncertainties.
As a result, the 1σ confidence level region becomes [0.59, 1.59] × 10 −38 cm 2 /oxygen, and the measured NCQE cross section is determined as σ measured NCQE = 0.74 ± 0.22(stat.)+0.85 −0.15 (syst.)× 10 −38 cm 2 /oxygen. (5) The measured NCQE cross section, the theoretical NCQE cross section [19], and the atmospheric neutrino flux predicted using the HKKM11 model [15] are shown in Fig. 3.The measured NCQE cross section is consistent with the flux-averaged theoretical NCQE cross section within the uncertainties.Furthermore, the measured NCQE cross section is consistent with the measurement in the SK pure-water phase within the uncertainties (1.01 ± 0.17(stat.)+0.78 −0.30 (syst.)× 10 −38 cm 2 /oxygen) [14].The systematic uncertainty on the measured NCQE cross section in this study is larger than that in the measurement of the SK pure-water phase.The reason is that we take the difference of secondary interaction models into consideration, conservatively estimated by the comparison among these models.The uncertainty will be reduced with better understanding of secondary interaction models in future.

VII. CONCLUSION AND FUTURE PROSPECTS
We reported the first measurement of the atmospheric neutrino-oxygen NCQE cross section in the Gd-loaded SK water Cherenkov detector.Using a 552.2 day data set, the NCQE cross section was measured to be 0.74 ± 0.22(stat.)+0.85 −0.15 (syst.)× 10 −38 cm 2 /oxygen in the energy range from 160 MeV to 10 GeV, which was consistent with the atmospheric neutrino-flux-averaged theoretical NCQE cross section (1.02 × 10 −38 cm 2 /oxygen) and the measured NCQE cross section in the SK pure-water phase (1.01 ± 0.17(stat.)+0.78 −0.30 (syst.)× 10 −38 cm 2 /oxygen).Moreover, from the comparison of three different secondary interaction models, we found that BIC and INCL++ provide a somewhat better fit to the observed data than BERT.
As described in Sec.II, we continue the observation with a 0.03% Gd-loaded SK detector, the phase known as SK-VII.Since the neutron-tagging efficiency in SK-VII is higher than that in SK-VI (35.6%) [6,26], more delayed signals can be detected, and the observed number of events can be accumulated faster in SK-VII than in SK-VI.Assuming that the neutron-tagging efficiency in SK-VII is about 60%, the statistics increases by about 1.4 times with the same live time as SK-VI.After one more year of observation in SK-VII the statistical uncertainty will reach the NCQE cross section measurement in 3. The measured neutrino-oxygen NCQE cross section, the theoretical neutrino-oxygen NCQE cross section [19], and the atmospheric neutrino flux predicted using the HKKM11 model [15].Vertical bars show the statistical uncertainty (short bar) and the total uncertainty (long bar).Horizontal bars show the 1σ from the mean (0.60 GeV) of the theoretical NCQE cross section multiplied by the atmospheric neutrino flux.
the SK pure-water phase, and the secondary interaction models will be able to be verified more precisely.Additional measurement using T2K's accelerator neutrino beam interactions in SK-Gd will help to further refine the physics models for the secondary interactions.

5 FIG. 2 .
FIG.2.The distributions of θC (left), Evis (center), and N delayed (right).Each color-filled histogram shows the expected events in BERT.Non-NC includes CC, spallation, and reactor neutrino events.Solid and dashed line show the total expected events in BIC and INCL++, respectively.In all distributions, θC is greater than 50 degrees, Evis is between 7.49 MeV and 29.49MeV, and N delayed is greater or equal to one.The values of the chi-square are also shown for each distribution.
FIG.3.The measured neutrino-oxygen NCQE cross section, the theoretical neutrino-oxygen NCQE cross section[19], and the atmospheric neutrino flux predicted using the HKKM11 model[15].Vertical bars show the statistical uncertainty (short bar) and the total uncertainty (long bar).Horizontal bars show the 1σ from the mean (0.60 GeV) of the theoretical NCQE cross section multiplied by the atmospheric neutrino flux.

TABLE I .
The expected number of atmosphric neutrino events in each secondary interaction model.The fractions are summarized in parentheses.The expected number of spallation, reactor neutrino, and accidental coincidence events are common to each model and are described in the text.

TABLE III .
[3]tematic uncertainties of the expected NCQE, NC non-QE, and CC events.NC non-QE and CC cross section are from Ref.[3].As for the secondary interaction and energy cutoff, only negative values are considered.Atmospheric neutrino flux, atmospheric neutrino/antineutrino ratio, data reduction, and neutron-tagging uncertainties are common to NC, NC non-QE, and CC and are described in the text.