Measurement of neutron production in atmospheric neutrino interactions at the Sudbury Neutrino Observatory

TL;DR

This paper measures neutron production in atmospheric neutrino interactions at the Sudbury Neutrino Observatory, comparing experimental data with Monte Carlo simulations to improve understanding of neutrino interactions and neutron yields.

Contribution

It provides the first detailed measurement of neutron multiplicities in atmospheric neutrino interactions at GeV energies, using data from SNO and comparing with simulations.

Findings

Good agreement between data and Monte Carlo simulations.

Neutron detection efficiencies of 15.3% and 44.3% for phases I and II.

Minor discrepancies at intermediate energies with less than 2σ significance.

Abstract

Neutron production in GeV-scale neutrino interactions is a poorly studied process. We have measured the neutron multiplicities in atmospheric neutrino interactions in the Sudbury Neutrino Observatory experiment and compared them to the prediction of a Monte Carlo simulation using GENIE and a minimally modified version of GEANT4. We analyzed 837 days of exposure corresponding to Phase I, using pure heavy water, and Phase II, using a mixture of Cl in heavy water. Neutrons produced in atmospheric neutrino interactions were identified with an efficiency of $15.3\%$ and $44.3\%$, for Phase I and II respectively. The neutron production is measured as a function of the visible energy of the neutrino interaction and, for charged current quasi-elastic interaction candidates, also as a function of the neutrino energy. This study is also performed classifying the complete sample into two pairs of event categories: charged current quasi-elastic and non charged current quasi-elastic, and $\nu_{\mu}$ and $\nu_e$. Results show good overall agreement between data and Monte Carlo for both phases, with some small tension with a statistical significance below $2\sigma$ for some intermediate energies.

Figures

Tables

Table 1. Table 1: Origin of neutrons produced by atmospheric neutrinos in SNO as predicted by the MC simulation. The processes below the single horizontal line correspond to the sources of secondary neutrons. The uncertainties in parentheses correspond to the MC statistical uncertainties.

Neutron origin	Fraction
Neutrino interaction	$33.0\left(0.2\right)\%$
Neutron inelastic	$34.9\left(0.2\right)\%$
$\pi$/$K$ inelastic	$15.0\left(0.1\right)\%$
Proton inelastic	$7.3\left(0.1\right)\%$
Hadron capture at rest	$6.4\left(0.1\right)\%$
$\mu$ capture at rest	$2.20\left(0.04\right)\%$
Photonuclear	$0.90\left(0.02\right)\%$
Other	$0.29\left(0.01\right)\%$

Table 2. Table 2: Percentage of events producing at least one neutron. The calculated uncertainties in parentheses corresponds to the MC statistical uncertainty.

Process	Fraction with at least
	one neutron produced
$\nu$ CCQE	38.4(2.2)%
$\overline{\nu }$ CCQE	99.9(0.1)%
$\nu$ CCOther	88.8(2.0)%
$\overline{\nu }$ CCOther	94.7(2.1)%
$\nu$ NC	84.8(1.8)%
$\overline{\nu }$ NC	82.4(2.3)%
$\nu$ total	61.5(1.1)%
$\overline{\nu }$ total	95.6(0.6)%
Total	69.5(0.8)%

Table 3. Table 3: Number of events selected in data by the different criteria for both phases together (top row) and fraction of interaction types and neutrino flavor in each selection together as calculated using MC. The quality cuts criteria select an inclusive sample of neutrino interactions; the CCQE criteria enhance CCQE events; the non-CCQE criteria enhance CCOther and NC events; and the electron- and muonlike criteria enhance the corresponding lepton type. Given the different FV cuts for CCQE and non-CCQE selections, some events do not fall in either of those two categories. The uncertainties in parentheses correspond to the MC statistical uncertainties.

Mode	Quality cuts	CCQE selection	Non-CCQE selection	Electronlike	Muonlike
No. events (data)	512	123	208	283	229
CCQE	51.1(0.5)%	64.5(1.2)%	28.7(0.6)%	47.4(0.7)%	55.6(0.8)%
CCRES	22.1(0.3)%	18.0(0.5)%	29.1(0.5)%	20.6(0.4)%	23.9(0.5)%
CCDIS	13.3(0.2)%	9.3(0.4)%	19.9(0.4)%	14.0(0.3)%	12.5(0.3)%
CC Other	0.18(0.02)%	0.15(0.04)%	0.34(0.05)%	0.15(0.03)%	0.21(0.04)%
NCES	0.23(0.03)%	0.20(0.05)%	0.23(0.04)%	0.20(0.03)%	0.26(0.04)%
NC Other	13.1(0.2)%	7.8(0.04)%	21.7(0.4)%	17.7(0.4)%	7.5(0.2)%
${\nu }_e$	48.9(0.5)%	50.2(1.0)%	49.4(0.8)%	74.9(0.9)%	17.5(0.4)%
${\nu }_{\mu }$	47.6(0.5)%	47.7(1.0)%	44.9(0.7)%	20.5(0.4)%	80.5(1.0)%
${\nu }_{\tau }$	3.5(0.1)%	2.1(0.2)%	5.7(0.2)%	4.6(0.2)%	2.1(0.1)%

Table 4. Table 4: Parameters adjusted in genie to estimate neutrino interaction systematic uncertainties. The parameters above the single horizontal line control the neutrino interaction cross section, while the ones below control the hadron transport models within the nucleus. See Ref. genie for more details.

genie label	Physical parameter	Nominal value	1$\sigma$ uncertainty
Cross sections
MaCCQE	CCQE axial mass	$0.990\mathrm{GeV}$	$-15\%+25\%$
MaCCRES	CC and NC resonance axial mass	$1.120\mathrm{GeV}$	$\pm 20\%$
MaCOHpi	CC and NC coherent pion production axial mass	$1.000\mathrm{GeV}$	$\pm 50\%$
MvCCRES	CC and NC resonance vector mass	$0.840\mathrm{GeV}$	$\pm 10\%$
R0COHpi	Nuclear size controlling pion absorption in Rein-Sehgal model	$1.000\mathrm{fm}$	$\pm 10\%$
CCQEPauliSupViaKF	CCQE Pauli suppression via changes in Fermi level	$0.225\mathrm{GeV}$	$\pm 35\%$
AhtBY, BhtBY	Higher-twist parameters in Bodek-Yang model scaling	$A=0.538,B=0.305$	$\pm 25\%$
CV1uBY	GRV98 PDF correction param in Bodek-Yang model	$0.291$	$\pm 30\%$
CV2uBY	GRV98 PDF correction param in Bodek-Yang model	$0.189$	$\pm 30\%$
Hadronization
AGKYxF1pi	Pion transverse momentum in AGKY model agky	See Appendix C of Ref. genie
AGKYpT1pi	Pion Feynman x for N$\pi$ states in AGKY model Ref. agky	See Appendix C of Ref. genie
FormZone	Hadron formation zone	See Appendix C of Ref. genie	$\pm 50\%$
Hadron transport
MFP_pi, MFP_N	Pion and nucleon mean free path	See Appendix C of Ref. genie	$\pm 20\%$
FrCEx_pi, FrCEx_N	Pion and nucleon charge exchange probability	See Appendix C of Ref. genie	$\pm 50\%$
FrAbs_pi, FrAbs_N	Pion and nucleon absorption probability	See Appendix C of Ref. genie	$\pm 20\%$

Table 5. Table 5: Summary of the different systematic errors propagated into the analysis. The first column details the source of systematic uncertainty. The second column is the 1 ​ σ 1 𝜎 1\sigma size of the propagated uncertainty or a reference to the relevant section if a single value cannot be given. The third column provides the 1 ​ σ 1 𝜎 1\sigma variation on the total number of produced neutrons per neutrino interaction. The fourth column is the method used to propagate the systematic uncertainty (see the text for details).

Systematic parameter	$\pm 1\sigma$ uncertainty	$1\sigma$ fractional effect	Type
High-energy scale	See Fig. 12	$0.7\%$	Shift
High-energy resolution			Smearing
Assumed $\cos \theta$ in $E_{\nu }$ reconstruction	See Fig. 5		Shift
Particle misidentification	$e=0\pm 5\%$, $\mu =4\pm 5\%$		Shift
Ring miscounting	$e=14\pm 14\%$, $\mu =11\pm 9\%$		Shift
High-energy radial bias	$28\mathrm{mm}$		Shift
High-energy radial resolution	$160\mathrm{mm}$		Smearing
Quality cuts efficiency	1.47%	$1.5\%$	Reweight
Neutron capture reconstruction	See Sec. VII.1.5		Shift, smearing & reweight
Neutron detection efficiency	See Sec. VII.1.6	15.9%	Reweight
Atmospheric neutrino flux	$\sim 15\%$	1.5%	Reweight
Neutrino interaction model	See Table 4	12.5%	Reweight
MC statistical error	–	1.9%	Reweight
Total	–	24.9%	–

Equations

$$P_N(n|\lambda )=\frac{e^{-\lambda }}{n!}{\lambda }^n.$$

$$P_Q(q_i|n_i)\times P_T(t_i(\vec{r})|n_i)$$

$$P_i^{hit}(q_i,t_i(\vec{r})|{\lambda }_i)=\sum _{n_i}{P}_N(n_i|{\lambda }_i)\times P_Q(q_i|n_i)$$

$$\times P_T(t_i(\vec{r})|n_i)$$

$$P_j^{unhit}=e^{-{\lambda }_j(\vec{x})}$$

$$\mathcal{L}(\vec{x})=\prod _i^{hit}{P}_i^{hit}(q_i,t_i(\vec{r})|{\lambda }_i(\vec{x}))\prod _j^{unhit}{e}^{-{\lambda }_j(\vec{x})}$$

$$\sum _{n_i>0}{q}_i\times P_N(n_i|{\lambda }_i)\times P_Q(q_i|n_i)$$

$$E_r^{\nu }=\frac{m_p^2-(m_n-E_b{)}^2-m_l^2+2(m_n-E_b)E_l}{2(m_n-E_b-E_l+p_l\cos {\theta }_l)}$$

$$\left[3{\left(\frac{E+E_0}{E_M}\right)}^2-2{\left(\frac{E+E_0}{E_M}\right)}^3\right]\otimes G(0,{\sigma }_E)$$

Full text

I Introduction

During the past few years, great advances in our understanding of neutrino interactions in the 100\text{\,}\mathrm{MeV}$$\sim$$10\text{\,}\mathrm{GeV} energy range have been achieved. Experiments like T2Kt2k , MiniBooNEminiboone , and MINER$\nu$Aminerva have shed light on the neutrino-nucleus interaction mechanisms. Nevertheless, the limited ability of the detectors used by these experiments to identify the neutrons produced in the neutrino interactions limits our understanding of the interaction processes. Development of neutron tagging techniques is useful for three main reasons. First, it would reduce atmospheric neutrino backgrounds in proton decay or supernova relic neutrino searches, boosting the sensitivity of current experiments. Second, it could help to distinguish neutrinos from antineutrinos in nonmagnetized detectors, since antineutrinos typically produce more neutrons. Third, it would provide crucial information on neutrino cross section models, which are the driving systematic uncertainty in neutrino oscillation experiments like T2K and NO$\nu$Anova and the future DUNEdune and Hyper-Kamiokandehyperk .

Water Cherenkov detectors have been proven to be of great value for solar and atmospheric neutrinos detection. Nevertheless, identification of neutrons generated in neutrino-nucleus interactions is challenging since it requires detection of the mega electron volt-scale deexcitation process that follows the neutron capture. Super-Kamiokande (SK) demonstrated that neutron detection is possible in water Cherenkov detectors sk_firstndet , with a detection efficiency of approximately $20\text{\,}\%$. In a later study, SK applied the new ability to measure the total number of generated neutrons in atmospheric neutrino interactions, as a function of the visible energy superk_neutrons . However, no comparison between interaction models and measurements is provided, and such a comparison does not currently exist in the literature. In addition, an inclusive analysis is performed, without distinction between different types of neutrino-nucleus interactions.

In this study, neutrons produced in atmospheric neutrino interactions are successfully identified with the Sudbury Neutrino Observatory (SNO), a heavy water Cherenkov detector. A measurement of the number of produced neutrons as a function of visible energy of the neutrino interaction for different neutrino interaction types is presented along with a comparison with a Monte Carlo (MC) model using genie genie ; genie_rw and geant4 geant4 . The number of produced neutrons as a function of reconstructed neutrino energy for charged current quasielastic events is also given. Finally, we study the potential for $\nu$ and $\bar{\nu}$ separation using neutron tagging.

This article is structured in the following way. A brief overview of the SNO detector is given in Sec. II, followed by a description of the MC model used in this analysis and a MC study in Sec. III. The reconstruction algorithms used to characterize the atmospheric neutrino interactions and neutron captures are explained in Sec. IV. The selection criteria for neutrino interactions and neutron captures are in Sec. V and Sec. VI, respectively. Sec. VII is dedicated to systematic uncertainties. The final measurements of neutron production in atmospheric neutrino interactions are presented in Sec. VIII, along with a comparison to results from SK. Sec. IX presents the final discussion and summary.

II SNO detector

SNO was a Cherenkov detector using heavy water located at a depth of $2092\text{\,}\mathrm{m}$ ($5890$ mwe) in INCO’s Creighton mine, near Sudbury, Ontario. The layout of the detector is shown in Fig. 1 and it consisted of a $6\text{\,}\mathrm{m}$ radius spherical acrylic vessel (AV) volume containing heavy water nested into an $8.4\text{\,}\mathrm{m}$ radius spherical structure instrumented with 9456 photomultipliers (PMTs) sno_det . The total mass of the detector enclosed in the PMT structure, adding the heavy and light water regions, was about $2.7\text{\,}\mathrm{k}\mathrm{t}$. The entire detector was suspended in a cavity and submerged in light water, which shielded against radioactivity from the rock. A cylindrical tube called the neck connected the inner part of the acrylic vessel with an external clean room, which served as the interface for filling and deploying calibration sources. The outer detector region featured 91 PMTs attached to the main structure but facing outward (referred as OWL), in order to provide a veto against external events. In addition, 8 PMTs (referred as NECK) were attached inside the neck, and 23 PMTs were suspended in a rectangular frame in the outer light water volume facing towards the neck region. The motivation was to veto possible light leaks occurring at the interface of the detector with the deck, and the flashes of light that were produced at the interface between the acrylic and the water surface.

The SNO experiment was designed for solar neutrino detection, and hence it was optimized for low-energy events. Neutron captures on heavy water provide a higher-energy signal than conventional water Cherenkov detectors. This increases their observable energy above the typical radioactive backgrounds, and allows a higher neutron detection efficiency. SNO was operated in three different phases. In Phase I (the D2O phase), the active volume was filled with pure heavy water. In Phase II (the salt phase), the heavy water volume was doped with chlorine in salt form (NaCl) at $0.2\text{\,}\%$ by weight, which considerably boosted the neutron capture cross section and signal energy. Finally, in Phase III, 3He proportional counters were deployed in the detector, which provided a completely independent means of neutron detection. However, this last phase is not used in the current analysis due to the added complexity to the geometry, which would require further study to determine the impact on our reconstruction of atmospheric neutrino interactions. The results reported in this analysis correspond to data collected during $337.25\pm 0.02$ days for Phase I and $499.45\pm 0.02$ days for Phase II.

III Neutron production and detection in atmospheric neutrino interactions

Production of neutrons in neutrino interactions is a complicated process that depends on neutrino-nucleon cross sections; on the interactions of the produced particles within the nuclear media, known as the final state interactions (FSIs); and on the hadronic interactions of the final state particles that propagate in the detector media. We differentiate two ways neutrons can be produced in atmospheric neutrino interactions:

In the final state of the neutrino-nucleus interaction (primary neutrons): this includes neutrons produced directly at the interaction vertex by antineutrinos, as well as those created due to FSIs. 2. 2.

As the byproduct of interactions of final state particles in the detector media (secondary neutrons): this includes neutron production due to hadronic inelastic scattering, photonuclear interactions of leptons and mesons, and muon captures.

The free neutrons propagate in the detector media undergoing nuclear collisions before they are captured. Since the energy of the produced neutrons is much higher than $1\text{\,}\mathrm{keV}$ (fast neutrons), they need to reach thermal energies (approximately $0.025\text{\,}\mathrm{eV}$) prior to being captured. The number of scatters they undergo strongly depends on the neutron energy. In heavy water, the energy loss is on average $44\%$ per collision, so the number of scatters for neutrons between $1\text{\,}\mathrm{MeV}$ and $1\text{\,}\mathrm{GeV}$ can range between 10 and 30, with higher-energy neutrons being more likely to exit the detector. Following the thermalization process, the neutron diffuses in the medium until it is captured. This diffusion is orders of magnitude slower than thermalization, so the neutron capture time is mostly determined by diffusion, which is specific to the capture material and independent of the energy at which the neutron was produced. Finally, the neutron is captured by a nucleus, which is left in an excited state and will deexcite, emitting particles on a very short timescale. The processes that could lead to a significant neutron detection in SNO are neutron captures on H, 2H, 35Cl, and 16O, with a subsequent emission of gamma rays of energies $2.2\text{\,}\mathrm{MeV}$, $6.25\text{\,}\mathrm{MeV}$, a cascade of $8.6\text{\,}\mathrm{MeV}$, and a cascade of $4.1\text{\,}\mathrm{MeV}$, respectively. Since the $2.2\text{\,}\mathrm{MeV}$ gamma-ray from H capture is below our analysis energy threshold, this detection channel is not relevant.

The entire process from neutron production to capture is simulated by our MC model. genie is used as a neutrino interaction generator, producing the final state particles, including primary neutrons. These particles are further propagated in the SNO geometry using geant4, which handles generation of secondary neutrons, neutron transport, capture, and gamma-ray emission. Finally, the detection process of gamma rays is handled by the snomansno_det detector simulation, which models the detector response. In the following section, we detail each stage of the simulation.

III.1 Generating neutrino interactions with GENIE

Atmospheric neutrinos interact in the different volumes of the SNO detector through charged current (CC) and neutral current (NC) interactions. Since the neutrino energies span several orders of magnitude, neutrinos will undergo several types of interactions: elastic scattering (ES), quasielastic (QE), resonant pion production (RES), deep inelastic scattering (DIS), or coherent scattering (COH). Pions and other hadrons will undergo a variety of FSI processes, such as: pion absorption, charge exchange, pion production, and elastic scattering, that modify the kinematics and nature of the original particles.

The neutrino interaction generator genie (version 2.10.2) is used to generate atmospheric neutrino interactions, the complex interaction models of which are described in Ref. genie , and the most relevant parameters for our analysis are summarized in Table 4. We input the unoscillated Bartol04 neutrino flux calculated for the SNOLAB location bartol04 and the SNO geometry and material composition for each phase. Neutrino oscillations are treated subsequently by reweighing the events. The total simulated data set contains 2 orders of magnitude more events than expected for the exposure of the analyzed data.

III.2 Secondary neutron generation and neutron propagation in GEANT4

The final state particles produced by genie are used as input into the geant4 tool-kit (version 10.0) geant4 , using the shielding physics list version 2.1. The same detector geometry used for genie is used in this step. The generation of neutrons is handled by a number of different models that simulate the processes: gamma photonuclear interactions; muon and electron nuclear processes; and inelastic scattering of mesons, protons and neutrons. Some of these processes have been compared against model predictions n_geant4_fluka ; n_geant4_proton . A limitation of GEANT4 is that it does not properly simulate deuteron photonuclear breakup. The impact of this process was estimated to be below $0.4\text{\,}\%$ by using an implementation of the original model developed for the SNO experiment sno_det . Neutron elastic scattering, crucial for the simulation of the thermalization process, is modeled using the NeutronHP package for energies below $20\text{\,}\mathrm{MeV}$ and the chips model for the higher energy range geant4 . This is a data-driven model that uses the Evaluated Nuclear Data File database. The relevant processes for neutron capture are also implemented in NeutronHP. A known problem with this model is that it randomizes the energy of the emitted gamma rays. As a result the sum of the total energy does not correspond to the actual total energy available for the deexcitation, violating energy conservation. This is not an issue for 2H and 3H, where a single energy state is present, but it is incorrect for 17O and 36Cl. A custom model based on the SNO implementation had to be introduced. We created a new neutron capture final state in our local geant4 installation that includes the actual branching ratios for 17O deexcitations and for 36Cl. The used energy levels and branching ratios for 36Cl are extracted from Ref. ncap_cl36 .

A breakdown of the origin of the neutrons produced along with their energy distributions is shown in Table 1 and Fig. 2, where we observe that roughly one-third of the neutrons is primary neutrons; one-third is produced as a result of neutron scattering, and one-third is due to other processes involving mainly protons, mesons and leptons. The energy of the produced neutrons ranges from a few mega-electron-volts to $1\text{\,}\mathrm{GeV}$, approximately $90\%$ of them being below $50\text{\,}\mathrm{MeV}$. The total number of produced neutrons in CCQE interactions, other CC interactions (CCOther) and NC interactions for neutrinos and antineutrinos is shown in Fig. 3. We observe that $69.5(0.8)\%$ of the neutrino interactions produces at least one neutron, as summarized in Table 2. On average, antineutrinos produce approximately one more primary neutron than do neutrinos in CC interactions, as can be seen at the bottom of Fig. 3. This difference is washed out by the production of secondary neutrons in CCOther interactions, but it still holds for CCQE interactions, highlighting the potential for $\nu$-$\bar{\nu}$ separation. The production of secondary neutrons is similar to the production of primary neutrons in CCQE interactions, but this is much larger in CCOther and NC interactions. The neutron production as a function of neutrino energy is shown in Fig. 4. Although the charged hadron production increases with the invariant hadronic mass, and hence neutrino energy genie_hadronic , the production of primary neutrons is practically constant over the entire energy range, and it is only the production of secondary neutrons that leads to an increase of the overall neutron multiplicity. According to our MC model, the fraction of neutrons that are produced within the AV and also captured inside the AV is $31.1\pm 0.3\%$ for Phase I and $74.4\pm 0.4\%$ for Phase II.

III.3 Detector simulation

The SNO detector is simulated with the package developed for the original SNO analyses, snoman sno_det . This package handles production and propagation of Cherenkov light in realistic detector conditions. The status of the electronics was recorded and simulated on a run by run basis, including the number of working PMTs and trigger conditions. Then, run-dependent efficiencies or reconstruction biases were modeled by snoman, which was extensively calibrated and validated using different deployed sources including AmBe and 252Cf to study the neutron detection response, 16N to calibrate the energy scale, and a diffused laser source to measure the optical properties of the detector sno_det . We also use snoman to simulate Cherenkov production from the final state particles produced by genie and geant4.

IV Event reconstruction

Two different classes of events need to be characterized: atmospheric neutrino interactions, which produce high-energy (approximately giga-electron-volt) leptons and hadrons in the final state with well-defined ringlike Cherenkov images in the detector, and neutron captures, which produce lower-energy (approximately mega-electron-volt) gamma rays that give a less well-defined Cherenkov signal. In order to properly deal with these different energy ranges, two event reconstruction algorithms are used and described below.

IV.1 Reconstruction of atmospheric neutrino interactions

The atmospheric neutrino reconstruction algorithm called Ring Fitter richie_thesis is designed to provide the position, direction, energy, particle identification (PID), and particle multiplicity from an atmospheric neutrino interaction occurring in the detector. The final state charged particles from a neutrino interaction are typically above approximately $50\text{\,}\mathrm{MeV}$, so the directional nature of the Cherenkov light creates well-defined ringlike structures. Characterizing these rings gives us critical information on the nature of the particle and consequently the neutrino interaction. The algorithm is based on the routines used by Cherenkov detectors such as MiniBooNE miniboone_fitter and Super-Kamiokande sk_fitter . In the following, we give an overview of the algorithm.

IV.1.1 Preliminary ring identification

We use the Hough transform technique hough to identify the center of the main ring in the spherical surface defined by the PMT structure. This will serve to give a preliminary estimate of the particle direction.

In order to obtain a first estimate of the event position, the fitter developed for the SNO+snop_nd water phase is used. Since it is optimized for low-energy events by design, its performance is poor at giga-electron-volt energies and it does not provide information on the particle type or multiplicity. The obtained position is used as a seed for the subsequent more complex algorithm.

The particle energy is also estimated at this stage by using the preliminary event position and the total amount of light collected in the event. This is done by building a lookup table using a complete MC simulation. Electrons and muons of energies up to $2\text{\,}\mathrm{GeV}$ and at different positions in the detector are generated using snoman. The result is a map of position and total charge vs energy.

IV.1.2 Determination of event position and direction

A likelihood fit is performed under the single-ring hypothesis to find the following observables related to the highest-energy particle, referred to as the main ring: event position $\vec{r}$, event time within the event window $t_{e}$ and event direction $\vec{d}$. The fit is run twice, once assuming an electron and again assuming a muon. The value of the likelihood in each case helps in identifying particle type, as described in the next section. The likelihood fit is based on the prediction of the number of photoelectrons (p.e.) that would be produced in each PMT for a specific position, direction, energy and particle hypothesis, represented by $\vec{x}=(\vec{r},\vec{d},t_{e})$. The probability of observing $n$ p.e. in a single PMT when $\lambda$ p.e. are expected is assumed to follow a Poisson distribution:

[TABLE]

For a given $n$, each PMT hit would present a different time and charge distribution, depending on its position with respect to the Cherenkov cone. The PMT time residual is defined as the PMT hit time corrected by the light’s time of flight assuming a position for the emission of the photon, which corresponds to $\vec{r}$. The probability of observing a hit $i$ with charge $q_{i}$ and time residual $t_{i}(\vec{r})$ for a given $n_{i}$ and $\vec{x}$ hypothesis will be the product of the charge and time probabilities, $P_{Q}$ and $P_{T}$:

[TABLE]

which are defined below. Then, the probability that a PMT $i$ with $\lambda_{i}$ expected p.e. records a hit with a given $q_{i}$ and $t_{i}$ is obtained by summing over $n$:

[TABLE]

If the $j$th PMT is not hit, then $n=0$ and the probability will simply be

[TABLE]

The likelihood function is obtained by multiplying the previous probabilities for all hit and unhit PMTs:

[TABLE]

For the PMT charge distribution we use the SNO single p.e. model and the averaged PMT gain measured at the detector. $P_{Q}(q_{i}|n_{i})$ is generated from MC using the measured signal p.e. charge distribution. On the other hand, the time distribution for single p.e. is parametrized as a prompt and prepulse peak, plus a uniform noise contribution and a flat scattering contribution for $t>$0\text{,}\mathrm{ns}$$. This distribution will be skewed towards earlier times for multi-p.e. hits, since the time registered by a PMT corresponds to the earliest photon. To model this effect, we created a two-dimensional probability distribution function (PDF) of $P_{T}$ as a function of $n$. This is done by extracting $n$ times from the single p.e. time distribution and populating the new PDF taking the time of the earliest p.e.

Estimation of $\lambda$ is done differently for muons and electrons. Muons created by atmospheric neutrino interactions are typically minimum ionizing particles during most of their range and suffer very little scattering. These two features are important since as a result the energy loss, path, and Cherenkov production per unit length are very reproducible for every muon; they typically travel on fairly straight lines, yielding a well-defined Cherenkov cone with a thickness proportional to their energy. Then, the Cherenkov yield and topology are determined very well by the position where the muon is created, along with its direction and energy. To estimate $\lambda$, we use a MC-generated PDF as a function of the PMT angle and distance from the muon track. For electrons, since their paths are much shorter, we approximate them as points. The angular dependence of the number of produced p.e. is calculated using the MC simulation for different electron direction and energy hypotheses.

Finally, we find the best fit value by floating $\vec{x}$ and using the minuit routine implemented in rootroot . We use the migrad algorithm to find the fit position and direction $\vec{x}_{f}$, for each of the two particle hypotheses.

IV.1.3 Particle identification and energy reconstruction

We identify whether the particle is electronlike or muonlike by exploiting the fact that the angular distribution of the emitted photons is much broader for electrons than for muons, due to the more pronounced electron scattering and secondary gamma-ray emission. We run the likelihood fit described above under the electron and muon hypotheses and calculate the likelihood difference $\Delta\mathcal{L}$ to determine particle type. The hypothesis with the best fit value is taken as the particle type. In cases where the fit for the position $\vec{r}$ is poor, the difference between the two hypotheses becomes small and the particle identification degrades. To overcome this problem, the likelihood is recomputed without the time residual term $P_{T}(t_{i}(\vec{r})|n_{i})$, and again, the hypothesis with the best fit value is chosen.

After the position, direction, and particle type have been precisely determined, we recalculate the particle energy by using MC-generated lookup tables for electrons and muons, binned in total PMT charge and radial position. The visible energy is defined as the electron-equivalent energy, i.e., the energy needed by an electron to produce the number of detected p.e. at the reconstructed radial position. The muon-equivalent energy is calculated in a similar fashion, and it is used to reconstruct the neutrino energy of muonlike events (see Sec. IV.1.5).

IV.1.4 Determination of ring multiplicity

Once the first ring has been identified and characterized, we predict the number of p.e. for each PMT and subtract them from the event. Then, the Hough transform is computed again in order to look for secondary rings. The predicted total charge for the $i$th PMT is defined by the average charge for the estimated number of p.e. $\lambda_{i}$ given by

[TABLE]

In order to reject false secondary rings, a Kolmogorov-Smirnov (KS) test against a flat background is performed. The used distribution is that of the PMT charge as a function of the angle between the PMT positions and the reconstructed center of the ring. An event is tagged as multiring if the total absolute charge and charge densities are above a certain threshold computed from MC and if the KS value is not significant. Should any of these conditions fail, the event is considered to be single ring.

IV.1.5 Estimation of neutrino energy

The neutrino energy is reconstructed according to the CCQE hypothesis,

[TABLE]

where $m_{p}$, $m_{n}$, and $m_{l}$ are the masses of the proton, neutron, and charged lepton, $E_{b}=$27\text{,}\mathrm{MeV}$$ is the effective binding energy of a nucleon in oxygen for leptonic interactions E_b , $E_{l}$ is the energy of the charged lepton, and $\cos\theta_{l}$ is the angle between the outgoing lepton and the incoming neutrino. Since the atmospheric neutrino direction is unknown, we estimate $\cos\theta_{l}$ from the genie prediction as the mode of the $\cos\theta_{l}$ distribution in a charged lepton’s energy bin (see Fig. 5). In this way, only the energy of the charged lepton is needed to estimate the neutrino energy. The uncertainties in these curves are computed by defining a symmetric region around the mode that encloses $68\%$ ($1\sigma$) of the events in each energy bin.

IV.2 Performance of reconstruction of atmospheric neutrino interactions

The Ring Fitter algorithm has been validated against MC simulation of single particles and neutrino interactions. Single muons and electrons are generated across the detector volume at energies between 0 and $2\text{\,}\mathrm{GeV}$. The energy resolution, position resolution, particle misidentification, and ring miscounting have been validated as a function of the energy and radius with electron and muon simulations. In the energy region of interest, the radial position resolution is $28\text{\,}\mathrm{cm}$ on average, the charged lepton energy resolution is below $7\%$, the particle misidentification rate is below $17\%$, and the rate of identification of single-ring events as multiring events is below $10\%$.

The reconstruction of atmospheric neutrino interactions was validated using simulated events by comparing the reconstructed radial position and neutrino energy with the true values. The bias in the radial position is very small and below the position resolution, as shown in Fig. 6(a). The bias in the reconstructed neutrino energy using the CCQE hypothesis is shown in Fig. 6(b). The CCQE events have a negligible bias of $7\pm 1.2\text{\,}\mathrm{MeV}$, while the other type of interactions exhibit a significant deviation, as expected since they do not obey the CCQE hypothesis.

IV.3 Reconstruction of neutron captures

To extract information on neutron captures, the official SNO reconstruction algorithms are used, which have been extensively validated with calibration sources. The position is reconstructed using the so-called path fitter, and the energy is measured by the ftk algorithm, described in Ref. ftk . These yield an approximately $15\text{\,}\%$ energy resolution and an approximately $20\text{\,}\mathrm{cm}$ position resolution for event energies of $6\text{\,}\mathrm{MeV}$, estimated using an 16N source sno_n16 .

V Selection of atmospheric neutrino events

Atmospheric neutrinos energies above $40\text{\,}\mathrm{MeV}$ are selected, so their interaction in the SNO detector produces charged particles well above the radioactive backgrounds. Atmospheric neutrino candidates are identified by criteria that start with the selection of events with more than 200 triggered PMTs (NHits). Additional cuts are designed to minimize instrumental backgrounds and external events (quality cuts). Finally, events are classified into different samples.

V.1 Quality cuts

We have designed a criteria to identify fully contained events, i.e. events of which the charged particles deposited their entire energy in the active volume of the detector. Our main backgrounds are external cosmic muons and instrumental events, both generating high NHits events. The former is eliminated by requiring fewer than three triggered OWLs. Events due to external light leaking into the detector were identified and eliminated by requiring that none of the NECK PMTs is triggered. Events due to random flashes of light created by the PMTs, electronic pickup or sparks produced by PMTs are largely reduced to less than $1\text{\,}\%$ of the final selection using dedicated low-level cuts relying on event topology and PMT charge and timing information. A spherical fiducial volume of less than $7.5\text{\,}\mathrm{m}$ radius is chosen, and events reconstructing at a larger radius are removed in order to eliminate events that reconstruct poorly, partially contained events, and the external cosmic muon contamination. A low-energy threshold of $50\text{\,}\mathrm{MeV}$ is also applied. This criteria result in 204 selected neutrino interaction candidates in Phase I and 308 in Phase II. The $\left(R/R_{AV}\right)^{3}$ distribution is shown in Fig. 7, where $R$ is the reconstructed radial position and $R_{AV}$ is the radius of the acrylic vessel. The visible energy distribution is shown in Fig. 7. The MC is normalized to match the number of selected atmospheric neutrino events in data in order to directly compare the shapes. The absolute MC normalization is irrelevant for this analysis.

V.2 Event classification

We divided the entire dataset into CCQE or non-CCQE and separately into $\nu_{\mu}$ or $\nu_{e}$. CCQE interactions are typically characterized by having a single charged particle in the final state. This would lead to single-ring events, so we rely on determination of ring multiplicity in order to enhance CCQE interactions (CCQE selection) or enhance CCOther and NC candidates (non-CCQE selection). For the former, we require a single-ring event within a reduced fiducial volume of $6.5\text{\,}\mathrm{m}$, while for the latter, we require just a multiring event. Hence, there are some events selected by the quality cuts that do not fall in any category. The PID capabilities of the reconstruction algorithm that separates showerlike events and tracklike events is sufficient to identify $\nu_{e}$ and $\nu_{\mu}$ interactions. The total number of selected events and the fraction of each component are shown in Table 3 for each selection.

VI Selection of neutron captures

To identify neutron capture candidates, we require an event with energy larger than $4\text{\,}\mathrm{MeV}$ within a certain fiducial volume and in time coincidence with the neutrino interaction candidate, described in previous section. The main backgrounds are 8B solar neutrinos, the high-energy tail of radioactive backgrounds, and events due to instrumental noise. The former two categories are eliminated by the coincidence criteria, and the latter is greatly reduced by the low-level cuts originally designed for the SNO analyses, which leave an accidental coincidence rate lower than $0.025\text{\,}\%$, as measured using randomly generated detector triggers. Production of unstable isotopes with lifetime and energy of the order of the neutron captures (like 12B) are expected to be more than 2 orders of magnitude smaller than that of neutrons sk_cosmo .

We select all events within $0.25\text{\,}\mathrm{s}$ after an atmospheric neutrino candidate. Given that the neutron capture lifetime is of the order of a few milliseconds, the impact of this cut on the neutron detection efficiency is negligible. Events outside a fiducial volume defined by a sphere with $6\text{\,}\mathrm{m}$ radius are rejected. Random coincidences are largely mitigated by the $4\text{\,}\mathrm{MeV}$ energy cut. We confirmed through an independent analysis that the detector trigger efficiency is well modeled above $4\text{\,}\mathrm{MeV}$. Finally, events with a $\Delta t<$10\text{,}\mathrm{\SIUnitSymbolMicro s}$$ are rejected in order to eliminate possible Michel electrons and low NHit instrumental backgrounds. We select 88 neutron capture candidates in Phase I and 388 in Phase II. The energy distribution and the distribution of the time difference with respect to the neutrino interaction are shown in Fig. 8 for both phases. The larger number of detected neutrons in Phase II is due the longer exposure and higher neutron detection efficiency with respect to Phase I. The MC is normalized to match the number of selected atmospheric neutrino events in data.

The total neutron detection efficiency was estimated from MC to be $15.3\text{\,}\%$ for Phase I and $44.3\text{\,}\%$ for Phase II. As shown in Fig. 9, it features a strong dependency on the radial position of the neutrino interaction. This is due to the fact that neutrons created close to the light water (large radius) are more likely to leave the AV and capture in H, yielding a $2.2\text{\,}\mathrm{MeV}$ gamma ray, which is below detection threshold. The neutron detection efficiency increases significantly for Phase II, as expected. The plateau region near the center of the detector is due to the larger neutron absorption cross section of 35Cl as compared to 2H. The obtained efficiency values are compatible with the original neutron detection studies in sno_ndet . The small differences are related to the fact that the energy of the neutrons produced by atmospheric neutrino interactions is typically higher than those produced by solar neutrinos, resulting in a higher chance of escaping the AV. The neutron detection efficiency decreases with energy since high energy neutrino interactions typically produce higher energy neutrons, which are more likely to exit the AV volume. In addition, the range of the particles produced in the neutrino interaction is larger at higher energies, so the production point of secondary neutrons could potentially be further from the neutrino interaction point inside the D2O volume, and therefore be closer to the AV. The modelling of the neutron detection efficiency is studied using dedicated 252Cf calibration data (see Sec. VII.1.6).

VII Estimation of systematic uncertainties

A number of possible sources of systematic uncertainties are considered and estimated using various calibration sources and control samples. We separated them in the following categories: detector-related systematic uncertainties, cross section model uncertainties and uncertainties on the atmospheric neutrino fluxes and oscillation parameters. They are described in detail in the following sections.

VII.1 Detector systematic uncertainties

VII.1.1 High-energy scale calibration

In order to characterize the detector response at higher energies and calibrate the Ring Fitter energy reconstruction algorithm, data from two different sources were used: Michel electrons and stopping cosmic muons. The former provide an understanding of the intermediate energy scale since they provide a well-known energy distribution with a sharp cutoff at $52.8\text{\,}\mathrm{MeV}$. The latter provide calibration of the $\mathrm{GeV}$ energy scale since cosmic muons have a characteristic energy loss of approximately $2.35\text{\,}\mathrm{MeV}\text{\,}{\mathrm{cm}}^{-1}$ in heavy water, so determination of the muon range provides a valuable calibration source for energies around approximately $1\text{\,}\mathrm{GeV}$.

Michel electrons are easily identified by looking for events with more than 100 triggered PMTs preceded by an event in a time window between 0.7 and $10\text{\,}\mathrm{\SIUnitSymbolMicro s}$. Instrumental backgrounds are reduced by requiring that $55\text{\,}\%$ of the triggered PMTs are within a $5\text{\,}\mathrm{ns}$ window. PMT afterpulsing also occurs on timescales of a few microseconds and therefore could introduce an energy bias. The after-pulsing probability was determined to be $1\text{\,}\%$ per p.e. To mitigate after-pulsing contamination, only Michel electrons that are preceded by stopping muons with less than 2500 NHits are selected.

The Michel electron candidates are reconstructed using the Ring Fitter (Sec. IV.1) and the visible energy distribution is fitted with the expected analytical form pdg18

[TABLE]

where $E$ is the energy which is constrained to $E<E_{M}$, $E_{M}=$52.8\text{,}\mathrm{MeV}$$ is the maximum permitted energy, $E_{0}$ is an energy shift correction and the last term represents a Gaussian smearing of width $\sigma_{E}$, which is interpreted as the energy resolution. The fit is done for data and for simulated Michel electrons generated using cosmic muons in snoman MC. The best fits are shown in Fig. 10 and correspond to an energy offset of $4.1\pm 4.1\text{\,}\mathrm{MeV}$ for data and $2.6\pm 0.7\text{\,}\mathrm{MeV}$ for MC with an energy resolution of $18.9\pm 4.7\text{\,}\mathrm{MeV}$ for data and $10.0\pm 0.65\text{\,}\mathrm{MeV}$ for MC. The energy bias is compatible between data and MC and the energy resolution for data is larger than predicted. The difference is attributed to the effect of unmodeled PMT after-pulsing, and to be conservative, it is propagated as a systematic uncertainty.

External stopping muons produce a Michel electron signal near the end point of the track, allowing estimation of the muon range within the detector active region using the Michel electron’s reconstructed position and the muon’s reconstructed direction. Stopping cosmic muons are selected by requiring only one Michel electron candidate following an external event with more than three triggered OWLs. Since we are interested in single muon events, we reduce the dimuon and shower component by requiring a maximum of 25 triggered external veto PMTs. The neutrino-induced muon component is reduced by selecting downward-going events with $\cos\theta>-0.5$, where $\theta$ is the reconstructed zenith angle of the muon. We measure $dE/dX$ as the Ring Fitter-reconstructed muon-equivalent energy divided by the estimated muon range. The $dE/dX$ distributions are shown in Fig. 11. We divided the dataset between low-energy (less than $1.35$ GeV) and high-energy (greater than $1.35$ GeV), in order to investigate any energy-dependent bias or resolution. Gaussian fits are performed for data and MC to estimate energy bias and resolution. The energy resolution is compatible between data and MC. We observe a small shift, which is attributed to a small difference in the averaged PMT gain used to reconstruct the energy along with possible misreconstruction of the muon track length. In the same fashion as was done with the Michel electron calibration, we err on the conservative side by propagating this difference as a systematic uncertainty.

The summary of the final energy biases and energy resolutions is shown in Fig. 12. The calculated energy bias is applied as a correction to data and MC. The differences between data and MC are propagated as a systematic uncertainty. To be conservative, the observed shift between data and MC is added to the uncertainty in the energy bias. The quadrature difference between the data and MC energy resolution is applied as a smearing to the MC, and the difference with the nominal MC is used to evaluate the systematic uncertainty.

VII.1.2 $E_{\nu}$ reconstruction

The uncertainty in the angle between the incoming neutrino and outgoing lepton induces a systematic uncertainty in the reconstructed neutrino energy calculated from Eqn. 7. The $1\sigma$ uncertainty is computed for every lepton energy bin, as shown in Fig. 5 and propagated into the final analysis.

VII.1.3 Atmospheric position bias and resolution

External cosmic muons enter the detector through the spherical structure that holds the PMTs and, hence, at a specific known radius. This is used as a control sample to study performance of the radial position reconstruction for data and MC. Cosmic muons are selected as described in Sec. VII.1.1. An extra cut to remove events with more than 4000 triggered PMTs is applied, in order to have clearer rings and to ensure that no other effect could inflate the estimation of the systematic uncertainty. The Ring Fitter algorithm is applied to these events in order to reconstruct the entrance radial position $R$. The agreement of the reconstructed radial position between data and MC is good, with $\left\langle R^{dt}\right\rangle-\left\langle R^{mc}\right\rangle=$-28\text{,}\mathrm{mm}$$, where $\left\langle R^{dt}\right\rangle$ and $\left\langle R^{mc}\right\rangle$ are the radial position averages for data and MC. The quadrature difference between the width of the radial distribution for data $\sigma_{R}^{dt}$ and MC $\sigma_{R}^{mc}$ is $160.0\text{\,}\mathrm{mm}$.

VII.1.4 Particle identification and ring multiplicity performance

We use the Michel electron and stopping muon candidate samples to test the performance on PID and ring multiplicity determination. The fraction of Michel electron events that misreconstruct as muon events is $11\pm 1\%$ for MC and $7\pm 3\%$ for data. For the stopping muons, $26\pm 4\%$ are tagged as electrons for MC, in good agreement with $28\pm 11\%$ for data. The difference is propagated in the analysis as a systematic uncertainty in the electron PID.

The rate of single particle events reconstructed as multiring events in the stopping-muon sample is $8\pm 2\%$ for MC and $19\pm 7\%$ for data. For the Michel electron sample, the number of events reconstructed as multiring corresponds to $1\pm 0.2\%$ for MC and $15\pm 7\%$ for data. These discrepancies are propagated into the analysis as a systematic uncertainty.

VII.1.5 Neutron capture energy and position systematic uncertainties

Reconstruction of the low-energy signal from neutron captures was extensively studied for the original SNO analyses sno_ndet . The systematic uncertainties associated with the capture position, the position resolution, the energy scale and the energy resolution were computed using dedicated calibration campaigns where the different sources mentioned above were deployed. Comparison between MC and data yields the systematic uncertainties propagated in this analysis sno_leta . The impact of these uncertainties is negligible compared to the rest of the systematic uncertainties.

VII.1.6 Neutron detection efficiency

The neutron capture efficiency for low-energy neutrons is characterized by the calibrations performed with a 252Cf source for both phases. The source was deployed at different radial positions, and the detection efficiency was measured and compared to the original MC simulation. It was found to agree within $1.9\%$ for Phase I and $1.4\%$ for Phase II, demonstrating that the neutron modeling built into snoman is well understood. We compared our simulation in geant4 to the one in snoman by comparing both models for single neutrons produced at different energies and reproducing the capture efficiency calculated for the 252Cf source. The estimated neutron detection efficiencies for both models agree within $1\%$ for energies below $10\text{\,}\mathrm{MeV}$ and within $3\%$ ($5\%$) for Phase I (II) at higher energies. To be conservative, we propagated the differences as systematic uncertainties by adding them in quadrature to the numbers extracted from the 252Cf calibration. The systematic uncertainty due to the detection efficiency for neutrons at energies relevant to this analysis is dominated by the width of the distribution at each energy and radius bin. The overall resulting systematic uncertainty is $15.9\%$ and is the dominant systematic.

VII.1.7 Quality cuts selection efficiency

External cosmic muons are used as a control sample in order to estimate the efficiency loss of the cuts described in Sec. V.1. Dark noise in the OWLs leads to valid events being rejected due to the OWL cut. This is estimated by measuring the OWL noise rate by randomly forcing the detector to trigger at a rate of $5\text{\,}\mathrm{Hz}$. Only $0.27\%$ of the forced triggered events have more than one OWL hit, and the random coincidence of 3 OWLs is below $0.05\%$. We conclude that the loss in efficiency due to this effect is negligible. A similar study is applied to the NECK PMTs concluding that none of these effects has an appreciable impact. The inefficiency of the quality cuts for the cosmic muon sample is $1.5\text{\,}\%$ for data and $2.1\text{\,}\%$ for MC, being compatible within statistical uncertainties. The quadrature difference between these two values is propagated as a systematic uncertainty.

VII.2 Neutrino interaction model uncertainties

The number of predicted primary neutrons depends on the interaction models. genie implements a system to vary the different parameters that impact neutrino cross sections and FSI. We change each relevant parameter by $\pm 1\sigma$, returning a factor for every single event, which is applied as an individual event weight. In this way, we obtain the $\pm 1\sigma$ boundaries for the number of predicted neutrons. The genie parameters of which the uncertainties have been propagated are shown in Table 4, classified in cross-section, hadronization or hadron transport model parameters. Their nominal values and $1\sigma$ uncertainty are also shown. For this work, we varied the axial and vector masses for the CCQE, CCRES, and NC interactions; the parameters in the Bodek-Yang model for DIS; the mean free path, absorption probability, and charge exchange probability for hadrons traveling through the nucleus; the parameters associated to the AGKY hadronization model agky ; and the one associated to the hadron formation zone. The uncertainty in the cross section model is the dominant of the three categories.

VII.3 Neutrino flux uncertainties

Uncertainties on the neutrino production model and the neutrino oscillation parameters affect the theoretical prediction of the neutrino flux at SNOLAB. The model uncertainties are mostly driven by the uncertainty in the composition and energy spectrum of the primary cosmic-ray fluxes and the solar modulation. These are provided by the Bartol Collaboration bartol04 . Uncertainties relating to neutrino oscillation parameters are included using the uncertainties provided by the PDG18 pdg18 . In addition, the oscillations depend on the production point of the neutrino, the uncertainties of which are estimated in Ref. vprod_height and included in the calculation of the oscillations.

The aforementioned parameters are shifted within $1\sigma$, generating a set of toy MC used to calculate the $1\sigma$ error bands of the neutrino energy spectra. Those boundaries are used to propagate the flux systematic uncertainties into the analysis by reweighting the different components and taking the difference with respect to nominal as the estimated effect of these uncertainties.

VII.4 Systematic uncertainties propagation and summary

The overall strategy of propagating systematic uncertainties consists of defining parameters that control the different uncertainties and redoing the analysis for different values of these parameters. The difference with the nominal value is interpreted as the size of the effect of the specific systematic uncertainty. There are three types of parameters depending on the nature of the propagation:

Shift: the parameter is shifted by $\pm 1\sigma$. 2. 2.

Smearing: the observable is smeared using a Gaussian of width equal to $1\sigma$. 3. 3.

Reweight: the event is given a weighted value, which corresponds to a $\pm 1\sigma$ deviation from the nominal parameter.

The considered systematic uncertainties are shown in Table 5, where the size of the $1\sigma$ uncertainty and its impact in the analysis are included, along with the propagation method. The fractional effect in Table 5 corresponds to the $1\sigma$ variation on the total number of produced neutrons per neutrino interaction. Bin by bin uncertainties are considered in the final measurement.

VIII Results

The number of neutron capture candidates after an atmospheric neutrino interaction is shown in Fig. 13 for both phases. The agreement between data and MC is good, although we identified four events with abnormally large neutron multiplicity in Phase II, compared to MC. Their energies and radial positions for the neutrino and neutron events are within the bulk of the population and the MC expectation.

After correcting for the calculated neutron detection efficiency shown in Fig. 9, we estimate the average number of produced neutrons as a function of the visible energy in each phase, as shown in Fig. 14. The error bars on the data correspond to the statistical uncertainties while the size of the MC boxes represent the systematic uncertainties listed in Table 5. The $\chi^{2}/\mbox{ndof (number of degrees of freedom)}$ values are $8.17/6$ for Phase I and $10.8/6$ for Phase II, which include bin-to-bin correlations and correspond to p-values of $0.23$ and $0.09$, respectively. We performed a consistency check by comparing the efficiency-corrected neutron production in MC (red band) with the true neutron production (green line). This shows an excellent agreement, demonstrating that the efficiency correction is properly applied. The figure separates out the number of primary neutrons (blue line) to show how the production is dominated by secondary neutrons at higher energies, as discussed in Sec. III. The measured neutron production shows good agreement between both phases, despite the different neutron detection efficiencies.

Based on the compatibility between phases, we performed an analysis on the combined dataset. The $\chi^{2}/\mbox{ndof}$ value on the average number of produced neutrons vs visible energy is $6.66/6$, which corresponds to a p-value of $0.35$. After classifying the full dataset as defined in Sec. V.2, the average number of produced neutrons is calculated and shown in Fig. 15 for each selection, allowing the study of neutron production for different interaction scenarios. The CCQE selection has a purity of $64.5\%$. For the non-CCQE selection, a purity of $71.3\%$ is achieved. Finally, the predicted neutron production for electronlike and muonlike events is overall in good agreement with the prediction. The neutrino energy is reconstructed for the CCQE-enhanced selection, and the neutron multiplicities are calculated with respect to this observable, as shown in Fig. 16.

We compared the total number of produced neutrons obtained by this work with the SK results superk_neutrons . Since our measurement of neutron production is a combination of light and heavy water, we estimated the neutron production in a SNO detector filled with light water, in order to compare to the SK results. We calculate the expected neutron production difference between light water and heavy water by generating neutrino interactions in two SNO configurations: one with the AV filled with heavy water (nominal) and another with the AV filled with light water. genie vertices are produced in each geometry, and the final state particles are propagated in geant4 as described in Sec. III. According to our MC model, the total neutron production rate inside the analysis FV is $9.8\pm 2.8\%$ larger for SNO with heavy water than for SNO with light water, driven by the larger production from neutron inelastic scattering. We estimated the neutron production in SNO with light water by scaling our measurement by $0.9$. In Fig. 17, we show the comparison of the SNO measurement with the SNO with light water estimation and the nominal SK results. Our results are reasonably in agreement with SK data.

VIII.1 Fit to primary and secondary neutrons

The production of primary and secondary neutrons as a function of energy is very different—secondary neutrons production is larger at higher energy, while primary neutron production is rather flat (see Fig. 4). We estimate the contribution of each component by defining two normalization parameters (one for primary and another one for secondary neutrons) and constraining them with a $\chi^{2}$ fit. The difficulty of this analysis resides in the large correlations between these two parameters, given the uncertainties on the neutron production. We can break the degeneracy by fitting the CCQE and non-CCQE samples together, since the ratio between primary and secondary neutrons is quite distinct for CCQE and non-CCQE interactions (see Fig. 3). Before the fit, the nominal distributions show a p-value of $0.19$. The best fit for the normalization factors is $0.41\pm 0.50$ for primary neutrons and $0.95\pm 0.25$ for secondary neutrons, with a best fit $\chi^{2}/\mbox{ndof}=14.4/12$. The fit was performed using stand-alone CCQE, non-CCQE, electronlike and muonlike selections. The case presented here is the one that yields the lowest relative uncertainties. The uncertainty on the primary neutron production parameter is driven by a combination of the small production of primary neutrons and large uncertainties on the low-energy bins caused mainly by the neutron detection efficiency. Fig. 18 shows the corresponding distributions before and after the fit. The difference with respect to the nominal prediction is small and features a p-value of $0.43$. The secondary production is compatible with the MC model prediction, while the fit prefers lower primary neutron production, being in slight tension with the nominal prediction. Similar fits to the different phases and selections yield compatible results. The systematic uncertainties described in Sec. VII and the bin-to-bin correlations are taken into account in the fit.

VIII.2 Potential for $\nu-\bar{\nu}$ separation

In Sec. III, we showed how the simulation predicts that antineutrinos typically produce more neutrons than do neutrinos. This effect is enhanced in the CCQE case, since secondary neutron production is minimal, and antineutrinos produce on average one more primary neutron than do neutrinos (see Fig. 3). This feature is exploited to explore identification of neutrinos and antineutrino events by studying the distribution of the number of detected neutrons. Two normalization parameters are defined for the neutrino and antineutrino components and a $\chi^{2}$ fit is applied to the CCQE selection. The distributions before and after the fit are shown in Fig. 19. It is important to notice the difference in shape between the two contributions, which breaks the degeneracy of the two components. We found a best fit value of $0.81\pm 0.37$ for the normalization of the antineutrino component, in good agreement with the unity. This shows that we can constrain the antineutrino component at the $46\%$ level.

On the other hand, by selecting events with one or more detected neutrons, we enhance the number of antineutrino events from $23.6\%$ to $34.4\%$, according to the MC simulation.

IX Summary and discussion

We have measured the number of produced neutrons in atmospheric neutrino interactions as a function of the visible energy using the SNO detector. The neutrino interactions have been classified as $\nu_{\mu}$ vs $\nu_{e}$, and a subset has been classified as CCQE-like vs non-CCQE, in order to study the neutron production in each sample. The predictions from a MC model built using genie and geant4 are in reasonable agreement with our measurements, although there are small tensions in certain energy regions. Data and MC are compatible within $2\sigma$ in the entire range and for every subsample. Comparison with published SK resultssuperk_neutrons shows a good agreement. We provided the neutron production as a function of the neutrino energy for CCQE events, showing that data and MC agree within $1\sigma$. We compared data to predictions of primary and secondary neutrons with a $\chi^{2}$ fit to the number of produced neutrons as a function of visible energy for the CCQE and non-CCQE selections. Our study of the separation of $\nu$ and $\bar{\nu}$ components using the number of detected neutrons shows that we can constrain the $\bar{\nu}$ component at the $46\%$ level and increase the purity of $\bar{\nu}$ events by $10.8\%$ by selecting neutrino events in coincidence with neutrons captures.

The projected future phase of SK with Gd-loaded water will be very interesting to better understand neutron production models. Furthermore, an experiment with larger statistics and higher neutron detection efficiency like ANNIE annie will be very valuable to precisely study different neutrino-nucleus interactions and neutron production models as a function of interaction kinematics.

Acknowledgements.

This research was supported by Natural Sciences and Engineering Research Council of Canada; Industry Canada; National Research Council Canada; Northern Ontario Heritage Fund; Atomic Energy of Canada, Ltd.; Ontario Power Generation; High Performance Computing Virtual Laboratory; Canada Foundation for Innovation; Canada Research Chairs program (Canada); U.S. Department of Energy Office of Nuclear Physics; National Energy Research Scientific Computing Center; Alfred P. Sloan Foundation; National Science Foundation; the Queen’s Breakthrough Fund; Department of Energy National Nuclear Security Administration through the Nuclear Science and Security Consortium (United States); Science and Technology Facilities Council (formerly Particle Physics and Astronomy Research Council) (United Kingdom); Fundação para a Ciência e a Tecnologia (Portugal). We thank the SNO technical staff for their strong contributions. We thank INCO (now Vale, Ltd.) for hosting this project in their Creighton mine. We also thank the Super-Kamiokande Collaboration for allowing us to use their data in our comparison.