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

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.


INTRODUCTION
During the past few years, great advances in our understanding of neutrino interactions in the 100 MeV∼10 GeV energy range have been achieved. Experiments like T2K [1], MiniBooNE [2], and MINERνA [3] 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νA [4] and the future DUNE [5] and Hyper-Kamiokande [6].
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 [7], with a detection efficiency of approximately 20 %. 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 [8]. 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 [9,10] and geant4 [11]. 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 ν andν 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 algo-rithms 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 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 m radius spherical acrylic vessel (AV) volume containing heavy water nested into an 8.4 m radius spherical structure instrumented with 9456 photomultipliers (PMTs) [12]. The total mass of the detector enclosed in the PMT structure, adding the heavy and light water regions, was about 2.7 kt. 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 D 2 O 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 % by weight, which considerably boosted the neutron capture cross section and signal energy. Finally, in Phase III, 3 He 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 ± 0.02 days for Phase I and 499.45 ± 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: 1. 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. 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 keV (fast neutrons), they need to reach thermal energies (approximately 0.025 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 MeV and 1 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, 2 H, 35 Cl, and 16 O, with a subsequent emission of gamma rays of energies 2.2 MeV, 6.25 MeV, a cascade of 8.6 MeV, and a cascade of 4.1 MeV, respectively. Since the 2.2 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 snoman [12] detector simulation, which models the detector response. In the following section, we detail each stage of the simulation.

A. 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. [9], and the most relevant parameters for our analysis are summarized in Table IV. We input the unoscillated Bartol04 neutrino flux calculated for the SNOLAB location [14] 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.

B. 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) [11], 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 [15,16]. 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 % by using an implementation of the original model developed for the SNO experiment [12]. Neutron elastic scattering, crucial for the simulation of the thermalization process, is modeled using the NeutronHP package for energies below 20 MeV and the chips model for the higher energy range [11]. 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 2 H and 3 H, where a single energy state is present, but it is incorrect for 17 O and 36 Cl. 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 17 O deexcitations and for 36 Cl. The used energy levels and branching ratios for 36 Cl are extracted from Ref. [17].  A breakdown of the origin of the neutrons produced along with their energy distributions is shown in Table I 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 GeV, approximately 90% of them being below 50 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 II. 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 ν-ν 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 [18], 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 ± 0.3% for Phase I and 74.4 ± 0.4% for Phase II.

C. Detector simulation
The SNO detector is simulated with the package developed for the original SNO analyses, snoman [12]. 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 252 Cf to study the neutron detection response, 16 N to calibrate the energy scale, and a diffused laser source to measure the optical properties of the detector [12]. 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 megaelectron-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.

A. Reconstruction of atmospheric neutrino interactions
The atmospheric neutrino reconstruction algorithm called Ring Fitter [19] 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 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 [20] and Super-Kamiokande [21]. In the following, we give an overview of the algorithm.

Preliminary ring identification
We use the Hough transform technique [22] 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+[23] 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 GeV and at different positions in the detector are generated using snoman. The result is a map of position and total charge vs energy.

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 r, event time within the event window t e and event direction 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 x = ( r, d, t e ). The probability of observing n p.e. in a single PMT when λ p.e. are expected is assumed to follow a Poisson distribution: 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 r. The probability of observing a hit i with charge q i and time residual t i ( r) for a given n i and x hypothesis will be the product of the charge and time probabilities, P Q and P T : which are defined below. Then, the probability that a PMT i with λ i expected p.e. records a hit with a given q i and t i is obtained by summing over n: If the jth PMT is not hit, then n = 0 and the probability will simply be The likelihood function is obtained by multiplying the previous probabilities for all hit and unhit PMTs: 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 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 λ 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 λ, 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 x and using the minuit routine implemented in root [24]. We use the migrad algorithm to find the fit position and direction x f , for each of the two particle hypotheses.

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 ∆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 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 ( 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 electronequivalent 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 A 5).

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 ith PMT is defined by the average charge for the estimated number of p.e. λ i given by 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.

Estimation of neutrino energy
The neutrino energy is reconstructed according to the CCQE hypothesis, where m p , m n , and m l are the masses of the proton, neutron, and charged lepton, E b = 27 MeV is the effective binding energy of a nucleon in oxygen for leptonic interactions [25], E l is the energy of the charged lepton, and cos θ l is the angle between the outgoing lepton and the incoming neutrino. Since the atmospheric neutrino direction is unknown, we estimate cos θ l from the genie prediction as the mode of the cos θ 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σ) of the events in each energy bin.
Angle of produced lepton cos θ l vs lepton energy in a CCQE neutrino interaction for muons (top) and electrons (bottom). The red dots show the mode of the cos θ distribution at each energy bin with the 1σ uncertainty.

B. 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 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 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.0 ± 1.2) MeV, while the other type of interactions exhibit a significant deviation, as expected since they do not obey the CCQE hypothesis.

C. 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. [26]. These yield an approximately 15 % energy resolution and an approximately 20 cm position resolution for event energies of 6 MeV, estimated using an 16 N source [27].

V. SELECTION OF ATMOSPHERIC NEUTRINO EVENTS
Atmospheric neutrinos energies above 40 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.

A. 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 % 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 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 MeV is also applied. This criteria result in 204 selected neutrino interaction candidates in Phase I and 308 in Phase II. The (R/R AV ) 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.

B. Event classification
We divided the entire dataset into CCQE or non-CCQE and separately into ν µ or ν 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 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 ν e and ν µ interactions. The total number of selected events and the fraction of each component are shown in Table III for each selection.

VI. SELECTION OF NEUTRON CAPTURES
To identify neutron capture candidates, we require an event with energy larger than 4 MeV within a certain fiducial volume and in time coincidence with the neutrino interaction candidate, described in previous section. The main backgrounds are 8 B solar neutrinos, the highenergy 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 %, as measured using randomly generated detector triggers. Production of unstable isotopes with lifetime and energy of the order of the neutron captures (like 12 B) are expected to be more than 2 orders of magnitude smaller than that of neutrons [28].
We select all events within 0.25 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 m radius are rejected. Random coincidences are largely mitigated by the 4 MeV energy cut. We confirmed through an independent analysis that the detector trigger efficiency is well modeled above 4 MeV. Finally, events with a ∆t < 10 µ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.   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 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 35 Cl as compared to 2 H. The obtained efficiency values are compatible with the original neutron detection studies in [29]. 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 D 2 O volume, and therefore be closer to the AV. The modelling of the neutron detection efficiency is studied using dedicated 252 Cf calibration data (see Sec. VII A 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.
A. Detector systematic uncertainties

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 MeV. The latter provide calibration of the GeV energy scale since cosmic muons have a characteristic energy loss of approximately 2.35 MeV cm −1 in heavy water, so determination of the muon range provides a valuable calibration source for energies around approximately 1 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 µs. Instrumental backgrounds are reduced by requiring that 55 % of the triggered PMTs are within a 5 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 % 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 A) and the visible energy distribution is fitted with the expected analytical form [30] 3 where E is the energy which is constrained to E < E M , E M = 52.8 MeV is the maximum permitted energy, E 0 is an energy shift correction and the last term represents a Gaussian smearing of width σ 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 ± 4.1) MeV for data and (2.6 ± 0.7) MeV for MC with an energy resolution of (18.9 ± 4.7) MeV for data and (10.00 ± 0.65) 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 θ > −0.5, where θ 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.

Eν 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σ uncertainty is computed for every lepton energy bin, as shown in Fig. 5 and propagated into the final analysis.

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 A 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 R dt − R mc = −28 mm, where R dt and R mc are the radial position averages for data and MC. The quadrature difference between the width of the radial distribution for data σ dt R and MC σ mc R is 160.0 mm.

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 elec- tron events that misreconstruct as muon events is 11±1% for MC and 7 ± 3% for data. For the stopping muons, 26 ± 4% are tagged as electrons for MC, in good agreement with 28±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 ± 2% for MC and 19 ± 7% for data. For the Michel electron sample, the number of events reconstructed as multiring corresponds to 1±0.2% for MC and 15±7% for data. These discrepancies are propagated into the analysis as a systematic uncertainty.

Neutron capture energy and position systematic uncertainties
Reconstruction of the low-energy signal from neutron captures was extensively studied for the original SNO analyses [29]. 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 [31]. The impact of these uncertainties is negligible compared to the rest of the systematic uncertainties.

Neutron detection efficiency
The neutron capture efficiency for low-energy neutrons is characterized by the calibrations performed with a 252 Cf 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 252 Cf source. The estimated neutron detection efficiencies for both models agree within 1% for energies below 10 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 252 Cf 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.

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 A. 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 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 % for data and 2.1 % for MC, being compatible within statistical uncertainties.
The quadrature difference between these two values is propagated as a systematic uncertainty.

B. 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 ±1σ, returning a factor for every single event, which is applied as an individual event weight. In this way, we obtain the ±1σ boundaries for the number of predicted neutrons. The genie parameters of which the uncertainties have been propagated are shown in Table IV, classified in cross-section, hadronization or hadron transport model parameters. Their nominal values and 1σ 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 [32]; and the one associated to the hadron formation zone. The uncertainty in the cross section model is the dominant of the three categories.

C. 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 cosmicray fluxes and the solar modulation. These are provided by the Bartol Collaboration [14]. Uncertainties relating to neutrino oscillation parameters are included using the uncertainties provided by the PDG18 [30]. In addition, the oscillations depend on the production point of the neutrino, the uncertainties of which are estimated in Ref. [33] and included in the calculation of the oscillations.
The aforementioned parameters are shifted within 1σ, generating a set of toy MC used to calculate the 1σ 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.

D. 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: 1. Shift: the parameter is shifted by ±1σ.
2. Smearing: the observable is smeared using a Gaussian of width equal to 1σ.
3. Reweight: the event is given a weighted value, which corresponds to a ±1σ deviation from the nominal parameter.
The considered systematic uncertainties are shown in Table V, where the size of the 1σ uncertainty and its impact in the analysis are included, along with the propagation method. The fractional effect in Table V corresponds to the 1σ 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 V. The χ 2 /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 pro-duction 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 χ 2 /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 B, 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.

Visible Energy [MeV]
We compared the total number of produced neutrons obtained by this work with the SK results [8]. 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 ± 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.

A. 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 χ 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 ± 0.50 for primary neutrons and 0.95 ± 0.25 for secondary neutrons, with a best fit χ 2 /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.

B. Potential for ν −ν 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 χ 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 ± 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 ν µ vs ν 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σ in the entire range and for every subsample. Comparison with published SK results [8] 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σ. We compared data to predictions of primary and secondary neutrons with a χ 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 ν andν components using the number of detected neutrons shows that we can constrain thē ν component at the 46% level and increase the purity of ν 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 [34] will be very valuable to precisely study different neutrino-nucleus interactions and neutron production models as a function of interaction kinematics.

Energy [MeV] ν
Reco.  Neutron production measurement in this work compared to SK published results [8]. Black dots correspond to the present work, with gray boxes representing systematic uncertainties and solid lines being the total uncertainties. The estimation of SNO with pure light water (see the text for details) is shown with diamonds. The nominal SK measurement with light water is marked with circles, and it only displays statistical uncertainties.