A with a suppressed dependence on the axial mass M A , which has not yet been determined. For a comprehensive prediction of the neutron emission associated with neutrino interactions, we establish a simulation of particle emission via nuclear deexcitation of 12 C, a process not considered in existing neutrino Monte Carlo event generators. Energy spectrum fitting for each neutron multiplicity gives g s A ¼ -0.14 þ0.25 -0.26 , which is the most stringent limit obtained using NCQE interactions without M A constraints. The two-body current contribution considered in this analysis relies on a theoretically effective model and electron scattering experiments and requires future verification by direct measurements and future model improvement. DOI: 10.1103/PhysRevD.107.072006
Various experiments have measured neutrino-nucleon interactions, and our understanding of these interactions gradually deepens. Among many neutrino interaction channels, the neutral-current quasielastic (NCQE) interaction contains fundamental and interesting information about nucleons. The NCQE interaction, ν l þ N → ν l þ N, where N denotes either a proton or neutron, does not change the lepton charge between the initial and final states. In contrast, the charged-current quasielastic (CCQE) interaction, ν μ þ n → μ -þ p, does. The CCQE interaction only involves isovector weak currents, while the NCQE interaction is sensitive to isoscalar weak currents. Therefore, searching for strange quarks existing as sea quarks in nucleons through their isoscalar contribution to the NCQE interaction is possible. In experiments, one measures the strange axial coupling constant g s A , which is the strange axial form factor at four-momentum transfer squared Q 2 ¼ 0. Since the Q 2 dependence of the axial form factor is parametrized by an axial mass M A , the measured value of g s A generally depends on the value of M A .
The BNL E734 experiment performed the first measurement of g s A using the NCQE interaction [1,2]. They used accelerator neutrinos and measured the ν þ p → ν þ p and ν þ p → ν þ p differential cross sections as a function of Q 2 . They confirmed a strong positive correlation between g s A and M A . They obtained g s A ¼ -0.15 AE 0.07 with the strong constraint of M A ¼ 1.061 AE 0.026 GeV, the world average at the time. In the 1970s and 1980s, various measurements from deuteron-target bubble chambers appeared to be consistent with obtained results of M A ∼ 1.0 GeV [3]. However, recent experiments using carbon and oxygen targets have found results as large as M A ¼ 1.1-1.3 GeV, and the discrepancy has become an issue [4]. It is becoming clear that a two-body current contribution, called twoparticle two-hole (2p2h), must be considered to explain this discrepancy [5][6][7]. A direct measurement of the 2p2h interaction has not yet been realized, so there is a modeldependent uncertainty. Therefore, it is difficult to determine a reasonable constraint on M A .
The MiniBooNE Collaboration measured the fluxaveraged NCQE differential cross section [8]. Assuming M A ¼ 1. 35 GeV, obtained from their CCQE analysis [9], they found g s A ¼ 0.08 AE 0.26. In this analysis, they did not simultaneously fit M A and g s A . Using the results provided by the MiniBooNE Collaboration, Golan et al. performed an independent simultaneous-fit analysis using the NuWro Monte Carlo event generator [10,11]. This analysis also took into account the 2p2h contribution and obtained M A ¼ 1.10 þ0. 13 -0.15 GeV and g s A ¼ -0.4 þ0.5 -0.3 , confirming a positive correlation between these parameters.
The strange axial coupling constant g s A corresponds to the strange quark-antiquark contribution to the nucleon spin, commonly represented by Δs. Several experimental results have been obtained using polarized-lepton deep-inelastic scattering: Δs ¼ -0.18 AE 0.05 from EMC [12,13], Δs ¼ -0.085 AE 0.018 from HERMES [14], and Δs ¼ -0.08 AE 0.02 from COMPASS [15]. These results rely on SUð3Þ f flavor symmetry. The SUð3Þ f flavor symmetry is violated by a maximum of 20%, in which case these results are shifted by AE0.04 [15]. This uncertainty is approximately equal to or larger than the statistical and systematic errors of the experiments mentioned above. It is clearly of interest to measure g s A ðΔsÞ in a way that is independent of SUð3Þ f flavor symmetry, namely by measuring the NCQE interaction.
One challenge in measuring g s A using the NCQE interaction is the strong correlation with M A . In the BNL E734 and MiniBooNE experiments, a proton target was used primarily because of the difficulty of measuring NCQE on a neutron target. The value of g s A significantly changes the ratio of proton and neutron NCQE cross sections. Therefore, a measurement exclusively on a proton target (or neutron target) depends highly on M A and other normalization uncertainties. Conversely, when measuring the ratio, the normalization cancels out, and we can measure g s A with only a slight dependence on M A . In practice, nucleons measured by detectors are affected by final-state interactions (FSI), nuclear deexcitation, and secondary interactions (SI). These effects somewhat smear the information about the target nucleon. Nevertheless, information about the target nucleons and g s A can be extracted by measuring the neutron multiplicity with high efficiency.
This paper aims to measure the neutron multiplicity of atmospheric neutrino NCQE interactions at KamLAND and to obtain g s A with a slight dependence of M A . In addition to the 2p2h contribution, the nuclear deexcitation process, which can emit neutrons, is considered in our analysis. Our paper is organized as follows: Sec. II describes the formalism of the NCQE interaction; Sec. III introduces a simulation of particle emission via nuclear deexcitation; Sec. IV describes the KamLAND detector and data analysis; Sec. V contains details of the Monte Carlo simulation; Sec. VI, the analysis method and results; and our conclusions are presented in Sec. VII.
The Llewellyn-Smith formula [16] is commonly used to describe CC and NC QE interactions. The hadronic current is composed of vector and axial parts. Assuming a dipole form, the axial form factors of CC and NC are expressed as follows:
where g A denotes the axial coupling constant, and the sign þð-Þ is for proton (neutron). A value of g A ¼ 1.2723 AE 0.0023 is determined by nucleon β decay experiments [17].
The strange quark contribution g s A only appears in the form factors of NC.
The relationship between the vector and electromagnetic form factors can be written as follows:
where θ W is the Weinberg angle, and the indices p and n represent the proton and neutron, respectively. The vector form factors for the proton and neutron (F pðnÞ 1;2 ) can be written in terms of the electric G E and magnetic G M form factors:
where M is the average of the proton and neutron masses.
The electric and magnetic form factors, G E and G M , are formulated from electron scattering data. The dipole form was commonly adopted in the past, as in the axial form factors. However, as the deviation from the dipole form became apparent, a more sophisticated parametrization, BBBA05 [18], has recently been used. The strange vector form factor F s 1;2 ðQ 2 Þ in Eq. ( 4) can be expressed assuming a dipole form:
where the vector mass M V ¼ 0.84 GeV is determined by electron scattering experiments. A global analysis of the polarized electron elastic-scattering experiments shows that the values of strange vector form factors are consistent with zero [19]. Thus, we set F s 1 ¼ F s 2 ð0Þ ¼ 0 in this analysis. Generally, the vector form factors can be precisely determined from high-statistics electron-scattering data. In contrast, the axial form factors are uncertain because they can be measured only through neutrino interactions.
As can be seen from Eq. ( 2), the extraction of g s A , the purpose of this paper, depends on both g A and M A . Since g A is precisely determined, the uncertainty in M A is the larger issue.
The strange axial coupling constant g s A significantly changes the relative proton and neutron NCQE cross sections with little change in the total cross section. Figure 1 shows the NCQE cross section on carbon per nucleon in NuWro [10]. For lower values of g s A , the neutron contribution to the total cross section becomes smaller while the proton contribution increases. This trend is also evident by looking at the neutron cross section as a fraction of the total NCQE cross section, as shown in Fig. 2. The value of M A changes the shape of the NCQE differential cross section and the overall cross section normalization. These changes are almost equal for proton-and neutron-target cross section contributions. Therefore, measurements of only the proton-target (or neutron-target) NCQE interaction depend highly on uncertainties in normalization factors such as M A and the neutrino flux. In contrast, measuring the neutron-target cross section as a fraction of the total NCQE cross section makes it possible to measure g s A with less dependence on these normalization factors. The nucleons measured by detectors are affected by FSI and SI, so it is impossible to strictly identify the target nucleons on an event-by-event basis. However, by measuring nucleon multiplicity, it is possible to statistically separate the contribution of target nucleons using the distribution, within the uncertainty of these nuclear effects. This method requires high nucleon detection efficiency.
In this analysis, we measured neutron multiplicity using KamLAND, which has a neutron detection efficiency of over 80%.
The default values of g s A adopted in common neutrinointeraction Monte Carlo generators are different: g s A ¼ -0.08 in NEUT version 5.4.0.1 [20,21], g s A ¼ -0.12 in GENIE version 3.00.06 [22], and g s A ¼ 0 in NuWro version 21.09. These differences change the neutron fraction of the total NCQE cross section on carbon by about 10%.
Nuclear deexcitation often occurs associated with neutrino-nucleus interactions. The typical excitation energy is about 20 MeV in the case of the 12 C target [23]. The excitation energy is higher than the separation energies of various particles, including neutrons, protons, and α particles. Various types of particles can therefore be emitted via deexcitation processes. It is important to predict these nuclear processes, especially for experiments measuring neutron multiplicity. However, current sophisticated neutrino Monte Carlo event generators, such as NuWro, NEUT, and GENIE, do not take them into account. Here, we have established a systematic method to predict nuclear de-excitation [24]. This method can be used with the results of neutrino Monte Carlo event generators. Since this study is intended for use in liquid scintillator detectors, including KamLAND, we only discuss the 12 C target.
Neutrino Monte Carlo event generators are event-byevent simulations, so we need an event-by-event deexcitation model to use them. We use two simulation software packages in this prediction, TALYS version 1.95 [25] and a modification of Geant4 version 10.7.p03 [26].
TALYS is an open-source software package for the simulation of nuclear reactions. It provides a complete and accurate nuclear reaction simulation up to 200 MeV, including fission, scattering, and compound reactions. Given any nucleus and excitation energy, it provides the branching ratios of all nuclear deexcitation processes. Although TALYS provides branching ratios, it does not perform event-by-event simulations.
Geant4, a widely-used software package for simulating the passage of particles through matter, makes it possible to do the event-by-event simulation. Within Geant4, "G4RadioactiveDecay" simulates nuclear deexcitation and radioactive decay. An event-by-event simulation of deexcitation decay chains is performed by loading the branching ratios obtained from TALYS into G4RadioactiveDecay with several modifications. In addition to the branching ratios from TALYS, various parametrizations related to the shell model, including excitation energies and spectroscopic factors, are necessary for the simulation.
B. Shell model picture of 12
In the simple shell model picture of the 12 C ground state, two nucleons lie in the s 1=2 shell, four nucleons lie in the p 3=2 shell, and no nucleon lies in the p 1=2 shell. When a nucleon in the p 3=2 shell disappears, the excitation energy is Cross section per nucleon (10
1. NCQE cross section on carbon per nucleon as a function of neutrino energy. The black and orange lines represent the neutrino cross sections on protons and neutrons, respectively. The green and blue lines represent the antineutrino cross sections on protons and neutrons, respectively. The solid (dashed) lines are the cross sections with g s A ¼ 0 (-0.3). These results are obtained using NuWro with M A ¼ 1.2 GeV [10].
A . The solid (dashed) line represents the neutrino (antineutrino) cross-section fraction. This result is obtained using NuWro with M A ¼ 1.2 GeV at 0.5 GeV neutrino energy. The red, violet, and blue vertical lines represent the default values adopted in neutrino Monte Carlo generators, NEUT [20,21], GENIE [22], and NuWro [10], respectively. Lower values of g s A lead to a lower neutron contribution to the total cross section.
zero, leading to no deexcitation. Assuming the same probability for all nucleons, the spectroscopic factors for s 1=2 and p 3=2 are 1=3 and 2=3, respectively. However, it is known that the actual spectroscopic factor of s 1=2 is smaller than this value because it is more tightly bound than p 3=2 . We adopt 0.296 for s 1=2 and 0.704 for p 3=2 from electron scattering data [27].
In a more precise shell model picture, the p 1=2 shell is partially occupied by a nucleon pair due to nucleonnucleon correlation. From various shell model calculations, this partial occupation, called the pairing effect, is expected to occur with a probability of 40 AE 10% [23]. Therefore, 20 AE 5% of the time, the disappearance of a single nucleon from the p 3=2 or p 1=2 shell will leave the residual nucleus,
Nucleon disappearance from the s 1=2 shell is more complicated than disappearance from the p shell. Because of the high excitation energy, typically more than the separation energies, we need to consider various particle emissions, including multistep processes as well as single-step deexcitations. The branching ratios for γ, α, n, p, deuteron (d), triton (t), and 3 He emissions are extracted from TALYS, including the full decay chains of the daughter nuclei. Since the excitation energy of an s 1=2 -hole is large, the impact of the pairing effect is negligible.
Figure 3 shows the branching ratios of 11 B Ã decay as a function of excitation energy calculated with TALYS. The spin-parity is J π ¼ 1=2 þ for single nucleon disappearance from the s 1=2 shell. At the typical excitation energy of 23 MeV, neutron emission accounts for about 65% of deexcitations. This process strongly affects the neutron multiplicity associated with neutrino-nucleus interactions. In contrast, the neutron branching ratio for 11 C Ã decay at a 23 MeV excitation energy is about 6%. This branching ratio is similar to that of proton emission for 11 B Ã .
The excitation energy of s 1=2 -hole state has a finite width and is commonly parametrized with a Lorentzian distribution. We adopt E ¼ 23 AE 1 MeV as the mean and Γ ¼ 14 þ10 -2 MeV as the FWHM width from electron scattering data [27,28]. We briefly mention how the uncertainty of these values affects the branching ratios at the end of this section.
We simulate the deexcitation decay chain event by event with Geant4 using branching ratios extracted from TALYS and the excitation energy distribution. The original Geant4 code does not treat emissions of tritons, deuterons, or 3 He, so we modified the code to implement these decay modes. The kinematics of the deexcitation process, such as separation energies and recoil, is taken into account properly.
We compare our prediction with experimental data and other predictions. The orange histograms show the predicted results from Hu et al. using TALYS [30], and the greens represent our results. deexcitation modes: 11 9 Be, and
Li. The published result does not distinguish between d and α, so for comparison, we calculate the relative branching ratios of n and d=α. Another prediction result from Hu et al. uses TALYS version 1.95 [30], the same version used in our analysis. The excitation energy and spin-parity configurations may cause the difference between Hu's result and ours. The branching ratio to n is the most important parameter in this analysis. Our result agrees with the experimental data within a relative uncertainty of 20%.
Figure 5 compares the measured and predicted branching ratios for 11 B Ã in the same excitation energy range with the experimental result from Yosoi et al. [31]. The 3 He branching ratio is not shown because it is less than 1%. The n branching ratios are consistent within a 20% relative uncertainty. There is a large difference in the single-step decay of triton, where the experimental result has a much larger value than the predictions. It is seen from Fig. 3 that such a high branching ratio can not be explained by the model implemented in TALYS. The authors also discussed this issue, but the causes are still unclear. Further checks are needed, such as validation experiments and model evaluations. We also confirmed a large difference in the multistep α decay. Our result gives almost 0% while others show about 5%. The α emission process is dominant at low excitation energies around 10 MeV, which lead to low α kinetic energies and low excitation energies of the daughter nuclei. Since the neutron separation energy of 7 Li is as high as 7.3 MeV, multistep α deexcitations do not contribute significantly to neutron emission. All these differences between our prediction and experimental results and with other predictions are considered model-dependent uncertainties.
We also compare the branching ratios of 11 C Ã with another prediction by Kamyshkov et al. using SMOKER [23]. The SMOKER code does not consider the deexcitation modes of d, t, and 3 He, which account for about 15% of the total. We therefore only compare the n, p, and α branching ratios. In contrast to 11 B Ã , neutron emission is a minor deexcitation mode, while proton emission is a major one. The total branching ratio for single-step and multistep neutron decays is 5.7%, while SMOKER predicts 13.8%. This difference is also treated as a model-dependent uncertainty.
Finally, we check the impact of the mean and width of the excitation energy distribution on the branching ratios. The relative changes in the branching ratios are within 15% when each parameter is changed within its uncertainty. We assign this uncertainty from the excitation energy in addition to the model-dependent uncertainty derived from Figs. 4 and5.
KamLAND, a 1000-ton liquid-scintillator (LS) detector, is located 1000 m underground in the Kamioka mine, Japan. The cosmic muon flux is suppressed by a factor of 10 -5 relative to sea level. The detector consists of an 18 m diameter stainless-steel spherical tank that defines the boundary of inner and outer detectors (ID and OD, respectively). The inner surface of the tank is instrumented with 1325 17-inch photomultiplier tubes (PMTs) and 554 20-inch PMTs facing the center of the detector. A 13 m diameter EVOH/nylon balloon is suspended containing 1000 tons of LS. The elemental composition of the LS is approximately CH 2 [32]. The space between the balloon and the tank is filled with nonscintillating mineral oil, operating as a buffer (BO). The OD is a cylindrical vessel filled with pure water. This region is instrumented with 140 20-inch PMTs, acting as a cosmic-ray muon veto. Further details of the detector are in [33].
The KamLAND data used in this paper are based on a total live time of 10.74 years, acquired between January The branching ratios of n are multiplied by a factor of 1=2. The green histograms represent our result using TALYS, and the orange histograms represent the prediction by Hu et al. using TALYS [30]. The experimental data in black are from Yosoi et al., and the authors also provide the predicted result using the CASCADE code [31]. The hatched histograms represent the branching ratios for single-step decays, and the open histograms represent those from multistep decays. balloon (IB) at the center of the KamLAND [34]. Period III (3.66 years of live time) refers to data during KamLAND-Zen 400 experiment from October 2011 to August 2015. After this period, we extracted the IB and refurbished the OD system in 2016 [35]. Period IV (1.52 years of live time) started after the OD refurbishment in April 2016.
After the end of period IV, we started KamLAND-Zen 800 by installing a new 190 cm radius inner balloon. The data acquired during KamLAND-Zen 800 is not included in this analysis.
KamLAND detects neutrino interactions via scintillation light. There is no threshold for scintillation light, unlike Cherenkov light. As a result, a scintillator detector like KamLAND can detect not only charged leptons and pions but also protons and neutrons with low-energy thresholds. Protons directly produce scintillation light through ionization, while neutrons are detectable via proton recoils and later through capture on nuclei. The primary energy deposition by the proton recoils occurs very quickly on the order of ns, while the capture has a much longer lifetime of several hundred μs, making it possible to perform delayed coincidence measurements. Since the NCQE interaction mainly emits protons and neutrons, this feature of scintillator detectors makes it possible for us to measure the NCQE interaction.
A neutrino interaction in KamLAND produces a prompt event caused by the energy deposit of charged particles and neutron recoils. Neutrons are then captured by protons (or 12 C) with a lifetime of 207.5 AE 2.8 μs [32], emitting a 2.2 MeV (4.9 MeV) gamma ray which produces a delayed event. We can observe the neutron capture events with high accuracy by performing delayed coincidence measurements using time and spatial correlations of prompt and delayed events.
We give some notes on the energy and vertex used in this paper. We use visible energy to evaluate the atmospheric neutrino events here. For CC events, the visible energy includes the energy deposit of the final-state lepton (electron or muon). On the other hand, in the case of NC events, the visible energy does not include that of the final-state lepton (neutrino), leading to lower-prompt visible energy than CC events. The vertex used in this paper is almost equivalent to the centroid of the energy deposition. Since we cannot distinguish the energy deposit of different particles produced by a neutrino interaction, we treat all the energy deposition at the same point source. A new fitter for reconstructing neutrino interaction points and end points of the charged particle is currently under development.
We select prompt events with visible energies (E prompt ) in the range of 50-1000 MeV, where the charge linearity of the PMTs and electronics has been confirmed by dye-laser calibration. Furthermore, NCQE and CCQE interactions are dominant in this energy region. We apply two spherical fiducial volume selection criteria with different radii: A 450 cm radius for 50 < E prompt < 200 MeV (low-E selection), and a 500 cm radius for 200 < E prompt < 1000 MeV (high-E selection). Because fast neutron events are present as a background below 200 MeV, we apply a tighter radius cut for the low-E selection. Detailed information about the fast neutron background is described in Sec V. We also apply OD cuts using the number of hit OD PMTs within a 200 ns time window N 200OD to cut cosmic muon backgrounds: N 200OD < 5 for periods I-III and N 200OD < 9 for period IV. Since we refurbished the OD system before the beginning of period IV, we adjust the threshold, so that veto efficiencies are equal. The OD cuts reject atmospheric neutrino events where the final-state particles exit the ID. All the events selected in this analysis are fully contained in the ID. Overall, we find 425 events for the high-E selection and 114 events for the low-E selection. The event rate in each period is stable within statistical errors.
We select delayed events, i.e., neutron capture gamma rays, using the delayed coincidence method. We use the radius (R delayed ), the time difference from the prompt event (ΔT), and the number of hit 17 inch PMTs within a 125 ns time window (N sumMax). We set R delayed < 600 cm, which is well inside the LS region (R < 650 cm). Immediately after a high-charge event, PMT afterpulses cause many noise events. The high event rate leads to channel-level electronics deadtime effects, and many PMT waveforms are not recorded, making accurate energy reconstruction difficult. Thus, we set 10 < ΔT < 1000 μs and exclude events with a time delay less than 10 μs. We select delayed events using N sumMax instead of the visible energy as it is less affected by these issues. We set N sumMax > 275 hits, a sufficiently low threshold to detect 2.2 MeV gamma rays.
Figure 6 shows the time difference between atmospheric neutrino events (prompt) and neutron capture events (delayed). The detection inefficiency that occurs for ∼50 μs immediately after atmospheric neutrino interactions can clearly be seen. The ΔT distribution is fitted with a function,
between 200 < ΔT < 1000 μs, where τ n ¼ 207.5 μs.
The constant term N const corresponds to the background contamination in delayed events. It is consistent with zero within a large uncertainty: N const ¼ 0.56AE 0.74 events=50 μs. The background event rate is also estimated using a long off-time window (2 < ΔT < 3002 ms). The result is ð3.61 AE 0.08Þ × 10 -2 events=50 μs. This low event rate means we have negligible contamination in the delayed events, ð0.160 AE 0.003Þ%. The neutron tagging efficiency ϵ can be calculated from the actual number of observed neutrons (N obs ) and the integral of the fit result,
estimating the inefficiency caused by the channel-level electronics deadtime effects. The selection inefficiency caused by radius cut is taken into account in the detector simulation described in Sec. V. The simulation also shows that the inefficiency associated with the gamma ray escaping the LS is insignificant. We obtain ϵ ¼ 89.7 þ8.4 -7.3 % from Fig. 6. It is known that the neutron tagging efficiency in KamLAND has a prompt energy and a time dependence. Since the leading causes of this inefficiency are afterpulses and the overshoots that occur ∼100 ns after a high-charge event, these effects depend on the charge intensity of the prompt event. In addition, PMT aging has gradually decreased the efficiency. However, due to low statistics, the analysis performed here using neutrons associated with atmospheric neutrino events cannot evaluate the prompt energy and time dependence. We therefore need an alternative way to estimate the efficiency more precisely. A more precise analysis using cosmic muons is described in Sec. IV C.
Figure 7 shows the spatial difference between the prompt atmospheric neutrino interaction and delayed neutron capture events. The spatial difference ΔR is the distance between the reconstructed positions of the center of energy deposition for the prompt and delayed events. Since the neutrons emitted via the neutrino interaction have high energy, the ΔR distribution spreads widely. The KamLAND data are compared with Monte Carlo simulation without any spectral fitting, with M A ¼ 1.2 GeV and g s A ¼ 0. The Monte Carlo simulation and KamLAND data are in good agreement. This consistency indicates that the Geant4 neutron transport model, used in the detector simulation, reproduces the data very well. The simulation details are described in Sec. V.
After cuts, we find 356 delayed events in the high-E selection and 91 delayed events in the low-E selection, with negligible background contamination. Note that the presence or absence of delayed events is irrelevant to the selection of prompt events.
As mentioned in Sec. IV B, the neutron tagging efficiency in KamLAND has prompt energy and time dependence. For a more precise analysis, we parametrize the neutron tagging efficiency as a function of prompt energy for each period. We use cosmic muons with high statistics as prompt events and apply the same selection criteria for the delayed events as in Sec. IV B. The method of calculating the neutron tagging efficiency is the same. For each prompt energy bin, the ΔT distribution is fitted with the function of Eq. ( 9). The obtained ΔT distributions are similar to Fig. 6, but differ in shape in the region ΔT < 150 μs, where the inefficiency occurs. Using the fit results, we calculate the neutron tagging efficiency according to Eq. (10). We confirm that the efficiency monotonically decreases over the experimental livetime of KamLAND and the prompt energy, within statistical uncertainty, as expected. Figure 8 shows the efficiency obtained as a function of prompt energy for period IV. The uncertainty is smaller than that obtained in Sec. IV B due to higher statistics. The energy dependence is parametrized with a second-order polynomial for each period, FIG. 6. Time difference between an atmospheric neutrino event (prompt) and a neutron capture event (delayed). All KamLAND atmospheric neutrino data sets are shown: both high-E and low-E selections during periods I-IV. The red line represents the fit result by Eq. ( 9) in the region 200 < ΔT < 1000 μs. The blue dashed line represents the selection criteria corresponding to ΔT < 1000 μs. FIG. 7. Spatial separation between atmospheric neutrino events and neutron capture events. All KamLAND atmospheric neutrino data sets are shown, as in Fig. 6. The red line represents the result of the Monte Carlo simulation before spectral fitting; the simulation assumes M A ¼ 1.2 GeV and g s A ¼ 0. The rightmost bin includes overflow. The simulation reproduces KamLAND data well.
where E prompt has units of GeV. The efficiency averaged over period I-IV is about 80% at E prompt ¼ 1 GeV and 88% at E prompt ¼ 0.1 GeV. These values are consistent with the result obtained in Sec. IV B, ϵ ¼ 89.7 þ8.4 -7.3 %. To take into account the prompt energy and time dependence of the efficiency, we use the values of p 0 , p 1 , p 2 , and error matrices under this parametrization in the fits to energy spectra described in Sec. VI.
Atmospheric neutrino events at KamLAND are estimated using Monte Carlo simulations. We use the atmospheric neutrino flux calculations of HKKM 2014 above 100 MeV [36] and Battistoni et al. below 100 MeV [37]. While seasonal variations of the flux are negligible, less than 1%, the effects of the solar cycle are not. We discuss this effect in Sec. VA. We calculate the neutrino oscillation effect propagating through the Earth using the Prob3++ package developed by members of Super-Kamiokande Collaboration [38]. The atmosphere is modeled as a vacuum. The Earth is modeled as a sphere of radius 6371 km with a simplified version of the preliminary reference Earth model (PREM) [39]. This simplified version of PREM has four layers with a spherical density profile. We use the three-flavor oscillation parameters assuming the normal hierarchy from [17]. Note that the neutrino oscillation affects only CC interaction, namely the background of this analysis. The uncertainty on the neutrino oscillation parameters gives no perceptible change in the sensitivity of this analysis.
We use NuWro version 21.09 to simulate neutrino interactions. For CCQE and NCQE interactions, the Llewellyn-Smith formalism with BBBA05 vector form factors is adopted. Resonant pion-production (RES) processes are simulated with the Adler-Rarita-Schwinger formalism [40,41] with dipole form factors [42]. The cross section discussed in this model has a 10% normalization uncertainty. For the 2p2h interaction, we choose the transverse enhancement model (TEM) [7], which is the only NC 2p2h model available in NuWro. This model has a 20% normalization uncertainty on the cross section. The TEM model does not predict the fraction of np pair targets in 2p2h interaction. In electron scattering, experiments have confirmed np pairs are dominant, with measured fractions of as 0.96 þ0.04 -0.22 [43] and 0.92 þ0.08 -0.18 [44], but the case of pure weak interactions is uncertain. The theoretical calculation for the weak interaction predicts 67% for the fraction of np pairs [45]. Based on the mean and deviation of these values, we assign 0.85 þ0. 15 -0.20 in this study. The cross section data used to model nucleon FSI is a custom fit model, which improves the agreement with the experimental data [46]. The one used in pion FSI is based on [47]. We use the local Fermi gas model, which is more accurate than the relativistic Fermi gas model.
After the neutrino interaction and nuclear deexcitation simulations, the detector response is simulated using a Geant4-based Monte Carlo simulation called KLG4. KLG4, which uses Geant4 version 9.6.p04, is a full optical detector simulation, including detailed descriptions of the KamLAND geometry and optical parameters. We adopted a hadron physics package called "QGSP_BIC_HP", suitable for sub-GeV hadronic interaction and precise thermal neutron transportation. The optical parameters, such as the light yield, quenching effect, and attenuation length, are tuned to reproduce the KamLAND data. We primarily used radioactive source calibration data for tuning, including 60 Co, 68 Ge, and 137 Cs sources. The quenching effect is parametrized by Birk's formula [48]. The energy peaks of these sources agree within 3.5%, and the vertex bias is less than 3 cm. We also estimate the energy scale uncertainty using spallation products of cosmic muons, 12 B and 12 N. Using the energy spectra of these β decays, which have endpoints at around 15 MeV, the uncertainty is estimated to be almost equal to that of the source calibration data. Finally, we checked the charge scale uncertainty for the high-energy region using minimum-ionizing cosmic muons. The charge peak of minimum ionization agrees within 8%, and the value is used in the fit described in Sec. VI. It is known that Birk's formula does not properly describe the quenching effects for heavier charged particles such as protons. The proton quenching effect of KamLAND LS is precisely measured using a monochromatic neutron beam [49]. The quenching factor obtained by the experiment is parametrized by a formula proposed by Chou [50] that empirically extends Birk's formula. KLG4 8. Neutron tagging efficiency as a function of prompt visible energy during period IV. The red (magenta) lines represent the best fit (1σ uncertainty) of the parametrization with the second-order polynomial of Eq. (11). Period IV has the lowest efficiency of the four periods. implements Chou's formula to describe the quenching effect for protons.
Fast neutrons induced by cosmic muons in the surrounding rock and water are a dominant background below 200 MeV. The neutrons scatter on protons and carbon nuclei in the LS, mimicking prompt events. Then, they are thermalized and captured on protons and carbon, creating delayed events. Neutrons produced outside the detector are exponentially attenuated by the shielding layers of water, BO, and LS. However, a contribution remains within the fiducial volume of this analysis. We estimate the fast neutron background using the KLG4 and cosmic muon profile at the KamLAND site [32,51]. The uncertainty depends on the neutron production yield in rock, and the simulation takes considerable computation. We conservatively assign a 100% uncertainty to our estimate.
HKKM 2014 provides atmospheric neutrino flux data at the solar minimum and maximum. The minimum and maximum are defined using the count rate of a specific neutron monitor (NM), the Climax NM [52]. This parameter is widely used to characterize the degree of solar activity. There is a linear and inverse correlation between these parameters. It is assumed that while the correlation gradient will depend on the location of various NMs, a linear correlation applies. HKKM defines 4150 counts=hour=100 as the solar minimum, and 3500 counts=hour=100 as solar maximum. From the Climax NM data trend, we can adequately consider the solar cycle's effect on the atmospheric neutrino flux. However, because the Climax NM was shut down in 2006, we need to calculate an equivalent Climax NM count, termed the NM parameter, using other NM data.
We use five NM datasets in addition to the Climax NM, which cover the entire analysis period: the Moscow, Apatity, Thule, Newark, and Oulu neutron monitors [53,54]. Their count rates have a linear correlation with the Climax NM data. We fit the correlation between each dataset and the Climax NM with a first-order polynomial during the period for which both monitors were available. We then convert the count rate of each monitor to the NM parameter, which is directly comparable to the Climax NM count rate. Figure 9 shows the trend of the NM parameter. Our data set indicates that solar cycle 24 had low solar activity. This result is consistent with that obtained by the Super-Kamiokande Collaboration [55]. We calculate the livetime-averaged NM parameter for each period as shown in Table I. The relative normalization change due to the solar cycle is calculated to be about 3%. The uncertainty of these count rates is 110 counts=hour=100 from the standard deviation of the five converted count rates.
Before discussing the fits to event energy spectra, we briefly introduce how g s A affects the KamLAND data. The effect of g s A appears as a change in the distribution of neutron multiplicities, while no apparent change is seen in FIG. 9. Time variation of the NM parameter. The shaded regions denote the four analysis periods in this paper. The dashed lines represent the solar minimum and maximum as defined in the flux calculation of HKKM 2014 [36]. The error bars are calculated from the standard deviation of the five converted count rates. the visible energy distribution. Since the NCQE interaction is dominant below 200 MeV, corresponding to the low-E selection in this analysis, the neutron multiplicity in that region is sensitive to g s A . Figure 10 shows the neutron multiplicity distribution of atmospheric neutrino events in the low-E selection (50 MeV < E prompt < 200 MeV). Since negative g s A increases the NCQE cross section with protons, the total cross section with KamLAND LS with its CH 2 composition also increases. The NCQE interaction with free protons is not accompanied by neutron emission via FSI and nuclear deexcitation and typically leads to zero neutron multiplicity. Thus, negative g s A enhances the rate of NCQE events with zero neutron multiplicity. Based on these considerations, we emphasize the importance of considering neutron multiplicity in the analysis.
We simultaneously extract M A and g s A from a fit of visible energy spectra. We used a binned χ 2 method incorporating systematic uncertainties. The χ 2 is composed of a Poisson term χ 2
Poisson and a penalty term χ 2 penalty :
The Poisson term is defined using the number of observed events n ijk and the number of expected events ν ijk ,
where the indices i, j, and k represent the ith period, jth visible energy, and kth neutron multiplicity bins. We have four data collection period bins corresponding to periods I-IV. We also have thirteen visible energy bins, eight for the high-E selection and five for the low-E selection. We divide the data into four neutron multiplicity bins, neutron multiplicity 0, 1, 2, and 3 or more. The analysis can consider neutron multiplicity by including the neutron multiplicity bins in the Poisson term. The penalty term is defined as
where l represents a systematic uncertainty parameter other than the neutron tagging efficiency, E l is the expected value, O l is the observed value in the fit, and σ l is expected uncertainty of the parameter l. The indicides n and m denote parameters of neutron tagging efficiency, and M -1 nm represents the inverse of error matrix described in Sec. IV C. The systematic uncertainties considered in this analysis are summarized in Tables II, III, and IV.
Figure 11 shows the best-fit visible energy spectra. The KamLAND data, which measure the neutron multiplicity with almost 80% efficiency, are well described by the simulations over a wide energy range, 50-1000 MeV. The NCQE interaction dominates in the low-energy region, roughly below 200 MeV. The neutron multiplicity in this energy region determines the value of g s A . Figure 12 shows the two-dimensional allowed regions for M A and g s A . We obtain M A ¼ 0.86 þ0. 31 -0.20 GeV and g s A ¼ -0.14 þ0.25 -0.26 . Our result is consistent with the result by Golan et al. using MiniBooNE data [11]. The plot shows little dependence on M A , as expected. This feature, realized by measuring neutron multiplicity, is important in the present experimental situation where measured values of M A vary from experiment to experiment. FIG. 10. Neutron multiplicity of atmospheric neutrino events in the low-E selection (50 MeV < E prompt < 200 MeV). The orangeshaded region represents the expected fast neutron background, the gray-shaded region shows expected atmospheric neutrino events from interaction modes other than NCQE, and the blue solid (dashed) lines denote NCQE interactions with g s A ¼ 0 (-0.30). The rightmost bin includes overflow. The simulation data are shown prior to the spectral fit, assuming M A ¼ 1.2 GeV. polarized-lepton deep-inelastic scattering experiments is not included. The second point concerns the treatment of M A and the 2p2h contribution. As described in Sec. I, it is difficult to determine a reasonable constraint on M A in the current experimental situation. The MiniBooNE and KamLAND results were obtained without M A constraints and included consideration of the 2p2h interaction. In contrast, the BNL E734 result did not consider any 2p2h interaction, and M A was strongly constrained. The BNL E734 result could therefore be affected by the contribution of the 2p2h interaction and a larger M A uncertainty. Our result gives the most stringent limit on g s A among NCQE measurements without M A constraints. The experimental NCQE data prefer smaller values than the results of polarized-lepton deep-inelastic scattering experiments and those adopted in neutrino Monte Carlo generators. However, they still have large uncertainties and are not yet accurate enough to claim adequate theoretical inputs. Further improvements in both experimental accuracy and theoretical modeling will be necessary.
In our analysis, the systematic and statistical uncertainties on g s A are almost comparable: g s A ¼ -0.14 AE 0.17ðstatÞ þ0. 18 -0.20 ðsystÞ. The dominant systematic uncertainties come from FSI and the 2p2h interaction, followed by the KamLAND neutron tagging efficiency. In order to improve the sensitivity to g s A , we need to optimize the FSI model. Electron scattering data can validate it, but it is not easy to model the dynamics of strong interactions in nuclei. Recently, S. Dytman et al. reported detailed comparisons of the FSI models implemented in NuWro, NEUT, and GENIE in [57]. They show significant variations between the generators, and further discussion is necessary. We also need to measure the 2p2h interaction directly and check the model's validity. Direct measurements of 2p2h interactions using detectors with good track resolution are planned and ongoing [58]. A combined analysis with the data from those experiments will be able to constrain g s A further. Nieves et al. [5] and Martini et al. [6] have been developing microscopic models to describe the 2p2h interaction. However, since they mainly focus on the CC, only the TEM is currently available for NC 2p2h in the generators. This situation has forced us to rely on the TEM in this analysis. The TEM is a theoretically effective model and relies on electron scattering experiments. The model uncertainties calculated in deriving the TEM are accurately considered in this analysis. Since the 2p2h interaction is due to nuclear effects, it is natural to consider an analogy with electron scattering experiments. Nevertheless, verification by direct measurement is required. We also expect that NC 2p2h models other than the TEM will be developed and implemented into the generators to allow verification of various models.
We report a new measurement of the strange axial coupling constant g s A using neutron multiplicity associated with the NCQE interaction of atmospheric neutrino at KamLAND. A simulation method for nuclear deexcitation, which is required to predict the neutron multiplicity accurately, is established. We use KamLAND atmospheric neutrino data from January 2003 to May 2018, corresponding to 10.74 years of total live time. By fitting the visible energy spectrum for each neutron multiplicity, we obtain FIG. 12. Two-dimensional allowed regions for M A and g s A . The red contour and dot are the result of this work. The side panels show the one-dimensional Δχ 2 -profiles projected onto M A and g s A . The violet contour and dot display the 1σ CL and best-fit value from Golan et al. using MiniBooNE data [11]. The parameter M A is treated as a free parameter in both results.
A (Δs) measurements. In addition to the result of this work, results from EMC [12,13], HERMES [14], COMPASS [15], BNL E734 [2], and Golan at al. using MiniBooNE data [11] are shown. Results with orange symbols are polarized-lepton deep-inelastic scattering experiments, and those with green symbols are neutrino NCQE scattering experiments. The red, violet, and blue vertical lines represent the default values adopted in several neutrino Monte Carlo generators.
g s A ¼ -0.14 þ0.25 -0.26 , which is the most stringent limit obtained using NCQE interactions without M A constraints. The experimental data on NCQE interactions, including this result, favor slightly smaller values than those used in the neutrino Monte Carlo generators. However, further improvements in accuracy are necessary to claim an appropriate value.
The main future tasks for enhancing the accuracy are detailed investigations of the FSI models and a direct measurement and model validation of 2p2h interaction. We need careful investigations on the FSI models to understand the large differences among the generators and to optimize the models to give good consistency with experimental data. Since only the TEM is currently available for NC 2p2h in the generators, this analysis has been forced to rely on the TEM, which relies on electron scattering experiments. Validation of the TEM by directly measuring the 2p2h is a high-priority future task. Although direct measurement of the 2p2h interaction is quite challenging at KamLAND, a combined analysis with other experiments, which aim at the direct measurements of the 2p2h, will be effective. We also expect that NC 2p2h models other than the TEM will be implemented into the generators.
In recent neutrino physics, the importance of accurately determining g s A and comprehensively predicting neutron multiplicity, including nuclear deexcitation, has increased dramatically. Detectors capable of measuring neutron multiplicity have been rare, but recent and next-generation detectors, such as Super-Kamiokande Gadolinium [59], Hyper-Kamiokande [60], and JUNO [61], will make it possible. These experiments plan to use neutron tagging information to significantly reduce the main background, atmospheric neutrino events, in searches for supernova relic neutrinos and proton decay. The dominant systematic uncertainty in these analyses comes from neutrino-nucleus interactions, especially the nuclear effects related to neutron emission. Therefore, the prediction accuracy of neutron emission in neutrino interactions affects the accuracy of these observations. Since the NCQE interaction of atmospheric neutrino is the main background to the supernova relic neutrino search, the determination of g s A will be essential. Next-generation detectors will significantly improve the measurement statistics, so reducing these systematic uncertainties is essential.
This analysis is the first to measure neutron multiplicity with a detection efficiency of ∼80%. It is also the first to compare measured neutron multiplicity with simulations that consider nuclear deexcitation. This analysis will add significant knowledge to the many recent and nextgeneration experiments mentioned above. All of them must consider nuclear deexcitation processes when conducting these studies. We expect to integrate the nuclear deexcitation simulation developed here into neutrino event generators for use in other experiments.
Tables II, III, and IV show the systematic uncertainties in this analysis and their best-fit values.