Cosmogenic Neutron Production at the Sudbury Neutrino Observatory

Neutrons produced in nuclear interactions initiated by cosmic-ray muons present an irreducible background to many rare-event searches, even in detectors located deep underground. Models for the production of these neutrons have been tested against previous experimental data, but the extrapolation to deeper sites is not well understood. Here we report results from an analysis of cosmogenically produced neutrons at the Sudbury Neutrino Observatory. A speciﬁc set of observables are presented, which can be used to benchmark the validity of GEANT4 physics models. In addition, the cosmogenic neutron yield, in units of 10 − 4 cm 2 / (g · µ ), is measured to be 7 . 28 ± 0 . 09 (stat.) +1 . 59 − 1 . 12 (syst.) in pure heavy water and 7 . 30 ± 0 . 07 (stat.) +1 . 40 − 1 . 02 (syst.) in NaCl-loaded heavy water. These results provide unique insights into this potential background source for experiments at SNOLAB.


I. INTRODUCTION
High energy muons created in cosmic-ray interactions in the Earth's atmosphere penetrate deep underground, where they induce electromagnetic and hadronic showers. These produce, among other particles of interest, free neutrons with an energy spectrum spanning several GeV. These cosmogenic neutrons form a direct background to searches for rare processes, such as neutrinoless double beta decay, nucleon decay, and dark matter interactions. The development and realization of next-generation detectors targeting these physics topics require unprecedented levels of background reduction. The prerequisite deep-underground location of such experiments reduces the rate of spallation backgrounds, but even the small number of remaining events can prove limiting to the potential physics reach of the experiments. It thus becomes critical to advance the understanding of the production and properties of cosmogenic neutrons. The average energy of the surviving cosmic muons increases with depth, and the extrapolation of cosmogenic neutron production rates from measurements made at shallow sites to greater depths is not well understood. Measurements at deep locations are critical to the success of future experiments.
The Sudbury Neutrino Observatory (SNO) experiment offers a unique data set to study cosmogenic neutron production deep underground. The SNO detector was a kiloton-scale heavy water detector, located at a depth of 5890 ± 94 m.w.e. Using the parameterization found in [15], the average muon energy at this depth is (363.0 ± 1.2) GeV , higher than those in many other published studies [1,2,[4][5][6][7][8][9][10][11][12][13][14], and comparable to that at LSD [3]. The SNO data can thus provide information in the high-energy regime, and further the understanding of how models for neutron production scale with muon energy.
Here we present results derived from the observation of cosmogenic neutrons in the SNO detector, namely a comparison of observables to model predictions and a measurement of the neutron production rate. Section II describes the SNO detector; Section III describes the Monte Carlo simulation used; Section IV describes the analysis methods, including the selection criteria for muons and neutrons, and backgrounds to this measurement; Section V presents comparisons of characteristic observables seen in the data to those predicted by simulations; and Section VI presents the results of the cosmogenic neutron yield measurement.

II. THE SNO DETECTOR
The SNO detector was a water Cherenkov detector located in INCO's (now Vale's) Creighton mine, near Sudbury, Ontario, at a depth of (2.092 ± 0.033) km. It consisted of a spherical acrylic vessel (AV) 12 m in diameter, filled with 1000 metric tons of 99.92% isotopically pure heavy water ( 2 H 2 O, or D 2 O). Surrounding the AV were 9456 Hamamatsu R1408 photomultiplier tubes (PMTs), each 20 cm in diameter, arranged onto a support structure (PSUP) of diameter 17.8 m. Each PMT was outfitted with a light-concentrator which increased the total photocathode coverage to approximately 55%. The AV was surrounded by 7.4 kt of ultra-pure H 2 O. The detector arrangement is shown in Figure 1.
Data taking proceeded in three phases. During Phase I, the inner volume was filled with pure D 2 O, with the neutron detection signal being the emission of a 6.25-MeV gamma following radiative capture on the deuteron. In Phase II, neutron detection was enhanced with the addition of 2 t of NaCl; 35 Cl has a larger neutron capture cross section, and a cascade of photons totaling 8.6 MeV in energy is emitted upon neutron capture. In Phase III, an array of 3 He proportional counters was deployed for neutron detection. The present analysis considers only data taken during the first two phases, with livetimes of 337.25 ± 0.02 and 499.45 ± 0.02 days, respectively.

III. MONTE CARLO SIMULATION
The existing SNOMAN Monte Carlo and analysis code [16] incorporates a detailed, high-precision model of the SNO detector, including geometry, material and optical properties, and the response of the PMTs and electronic readout system. This model was based on measurements of microphysical parameters, and tuned and verified using calibration data from deployed radioactive and optical sources in the context of previous neutrino analyses [17][18][19][20][21][22][23][24][25][26]. However, the code relevant to the production and propagation of muons and neutrons evolved to become a compilation of algorithms from various sources. In particular, neutron propagation was based principally on the MCNP package [27], which in SNOMAN is applicable only for neutron energies below 20 MeV. For the purposes of both improved accuracy in the high-energy regime, and ease of interpretation by the scientific community, in the present analysis SNOMAN is used only for the purposes of modeling detector response and event reconstruction in the context of measuring the neutron yield; the propagation of muons and neutrons is performed in GEANT4 [28] (version 10.00.p02), using the standard "Shielding" physics list with two modifications described below.
In the course of this analysis, two issues concerning the treatment of deuterons by the standard physics processes included in the Shielding list were discovered. One of the most prominent neutron-producing reactions relevant to this analysis is the photonuclear reaction γd → pn, which can occur in electromagnetic showers initiated by a cosmic muon. GEANT4 tabulates photonuclear cross sections as a function of the mass number of the nucleus, but, when calculating the cross section for a given isotope, uses a mass number corresponding to the average mass of the naturally occuring isotopes of the given element. For heavy isotopes of hydrogen, this incorrectly returns the cross section on a free proton, which for energies below the pion threshold is 0, as no nuclear break-up can occur for a single nucleon. This issue has been reported to the GEANT4 development team and has been corrected in release version 10.5. In this work, a patch was implemented to disable this behavior for deuterons, for which a cross section tabulation already exists.
It was further discovered that the default model for photonuclear final state generation, the Bertini Intranuclear Cascade, fails to properly model photodisintegration of the deuteron below the pion threshold. Indeed, while γd → γγd and similar reactions occur, γd → pn reactions do not. For the present analysis, we reimplemented the deuteron photodisintegration model developed for SNOMAN [29] as a GEANT4 physics process, which is applied only to γd reactions below the pion threshold. In short, this model treats deuteron break-up as a two-body problem subject to conservation of energymomentum. A summary of the contributions of various cosmogenic neutron-producing processes in GEANT4 is shown in Table I. TABLE I. Breakdown of cosmogenic neutron producing processes at SNO, as modeled by GEANT4. All processes labeled "inelastic" refer to inelastic scattering, and "µ-nuclear" refers to direct muon-nucleus interactions via virtual photon exchange.
The first step in the Monte Carlo is to generate muons on a spherical shell approximately 4 m outside the PSUP. Given the spherical geometry of the SNO detector, the track can be specified using three coordinates: the impact parameter, which is the distance from the center of the detector to the midpoint of the line connecting the entrance and exit points; the zenith angle, which is the angle of the track measured from vertical; and the corresponding azimuthal angle. The impact parameters and entrance angles are sampled from the muons reconstructed in data, convolved with the resolution of the muon track reconstruction algorithm used in previous cosmic analyses [30]. The initial muon energy is sampled from an analytic form taken from [15], namely where b = 0.4/km.w.e., ε = 693 GeV, γ = 3.77, are constants which parameterize the shape of the spectrum, h = 5.89 km.w.e./ cos θ is the slant depth parameterized by the incident zenith angle θ, and A is the normalization. This distribution is the result of propagating muons from surface [31], neglecting their angular dependence, through a depth h, in the approximation of continuous energy loss. While the angular dependence of the energy spectrum at surface is neglected, the angular dependence due to the flat rock overburden is a larger effect, and is included. The propagation of muons and all daughter particles is handled by GEANT4, subject to the two corrections to photonuclear reactions described above. To mitigate poor performance due to the great number of low-energy photons created by high-energy muons, optical photon tracking is disabled and no detector response is simulated. All observables extracted from the Monte Carlo are thus taken as truth information, as output solely of the physics models.

IV. ANALYSIS
There are two goals of this study. The first is to provide a detailed comparison of the data to model predictions across a number of observables, including the capture time and the reconstructed position of the captured neutrons, offering validation of the models implemented in GEANT4. The second goal is a measurement of the neutron yield, defined as the number of neutrons produced per unit muon track length per unit target material, in the D 2 O target.
Use of a heavy water target in SNO offered a higher energy signature for neutron capture than the morecommonly used light water and liquid scintillator: neutron capture on the deuteron results in a 6.25-MeV gamma, in comparison to the 2.2-MeV gamma from capture on hydrogen. As a result the efficiency for detecting neutron capture events is greater than 95% in the data set under consideration (Sec. IV F). The signal energy is also well above internal radioactive backgrounds, leading to effective neutron identification. Conversely, due to the relatively low muon flux at this depth, the data set is limited in statistics when comparing to studies performed of shallower sites.

A. Muon reconstruction
The reconstruction of a muon candidate event is performed under the through-going hypothesis and outputs several parameters that specify the muon track, including the impact parameter (b) and zenith angle (θ).
The details of the reconstruction algorithm are described in [25]. The reconstruction is performed in two stages, where a preliminary fit from the first stage is used as the seed to a more sophisticated algorithm in the second stage. The first stage is a purely geometric construction: the entrance point is identified with the cluster of earliest hit PMTs, and the exit point with the charge-weighted position of all hit PMTs. The second stage, which takes this seed track as input, is a likelihood fit containing terms for the number of detected photoelectrons, and the PMT multiphotoelectron charge and hit-times. Using an external muon-tracking system to validate the fits, the muon reconstruction algorithm was found to perform with a resolution of less than 4 cm in impact parameter and 0.5 • in zenith angle [30].

B. Data selection
The data used in this analysis was taken during Phases I and II, with the AV filled with pure heavy water and salt-loaded heavy water, respectively. It is thus a subset of the data used in the SNO cosmic muon flux measurement [25], which also considered data taken during Phase III, and the 13-day period between Phases II and III when the detector contained pure heavy water. Phase I data was collected between November 2, 1999 and May 28, 2001, and Phase II data was collected between July 26, 2001 and August 28, 2003, for a combined livetime of 836.7 ± 0.03 days.
The selection criteria for muon events are designed to select through-going muons and reject instrumental backgrounds. Specifically, to qualify as a muon, events must have had at least 500 calibrated PMTs fired, with fewer than three of them in the neck of the AV, which is a characteristic of external light entering from the top of the detector. Events that occur within 5 µs of another event in which 250 PMTs fired, or within a 2-s window containing 4 or more such events, are identified as a class of instrumental events called "bursts," and are removed from the analysis. Furthermore, events with uncharacteristically low total PMT charge and/or broad timing distributions are inconsistent with the muon hypothesis, and are similarly identified as instrumental events. Further high-level cuts are made, among which are the requirements that the reconstructed impact parameter b < 830 cm to ensure the validity of the track fit, and the reconstructed energy loss −dE/dX ≥ 200 MeV/m to reject muons that stop inside the detector volume. Finally, cuts are imposed on the fraction of photoelectrons geometrically contained inside the predicted Cherenkov cone for the muon track, and on the timing of these in-cone photons.
These criteria are identical with previous cosmic muon analyses [25,30] with one exception. A Fisher discriminant was previously used to reject stopping muons, but was found to incorrectly exclude muons with high light production -potentially the most interesting from the standpoint of neutron production -from the analysis. For the present analysis, we omit this linear discriminant cut; stopping muons are unlikely to contaminate neutron selection due to their relatively prompt decays, as discussed below. A total cross-sectional area of 216.4 m 2 is considered in this analysis, for which Monte Carlo studies of cosmic muons in SNOMAN show the total event selection cut efficiency to be greater than 99% for throughgoing muons [25].
The average capture time for thermal neutrons is known to be on the order of tens of ms in pure D 2 O, and was decreased to a few ms with the addition of NaCl in Phase II. We thus search for cosmogenic neutrons in a time window of 20 µs < ∆t < ∆t max following any through-going muon. The lower bound of 20 µs was chosen both to exclude Michel electrons from the decay of daughter muons from pions produced in hadronic showers, and to veto a period of several µs following particularly energetic muons in which the PMTs experienced significant afterpulsing. Imposing this lower bound reduces the livetime for neutron selection by less than 0.5%. The upper bound ∆t max was chosen to accept > 99% of neutron captures in each phase, and is set to 300 ms in Phase I and 40 ms in Phase II. Low-level cuts to identify candidate events are identical to those used in previous analyses [21,22,24]. Neutron events are identified by reconstructing Compton scatters of the capture gammas under a single-scatter hypothesis, yielding a total effective electron energy E eff and reconstructed radial position r. Neutron events are selected by requiring 4.0 MeV < E eff < 20.0 MeV and r < 550.0 cm.
These high-level selection criteria differ from previous neutron selection in using a widened energy window consistent between the two phases, compared to the 6-10 MeV window used previously for Phase II data [18], intended to maximize neutron acceptance. Table II shows the number of muons accepted for the cosmogenic neutron search, and the percentages for which a follower was detected in both the data and Monte Carlo. The scarcity of neutron followers as shown in the table results in fewer than 3000 muons with detected neutron followers across both phases.

C. Tests of model predictions
In order to validate the models of cosmogenic neutron production and propagation in the GEANT4 Shielding physics list at SNO depth and muon energies, we compare the data with model predictions for a number of observable distributions, including the properties of muons after which neutrons were observed, detected neutron multiplicity, neutron capture position, capture distance from the muon track, clustering of capture positions, and capture time.
These quantities offer benchmarks of different aspects of the models implemented in GEANT4, and unique measurements of the physics involved in neutron production. For example, measurement of the per-muon neutron multiplicity yields insight into the validity of the cross sections of different neutron-producing reactions, while the capture time is sensitive to different neutron energies. Understanding these complementary observables in the simulations and the data will lead to improved physics modeling, imperative for more precise physics measurements.
Furthermore, a measurement of the neutron production rate, using Monte Carlo information as input, requires the reliable simulation of several effects: direct and secondary production of neutrons, typically through electromagnetic and hadronic channels; the energy spectrum of produced neutrons, which can range up to several GeV; the transport of neutrons both at high and thermal energies; and the detection of capture gammas.
As the neutrons are thermalized and then detected after radiative capture, this analysis is not directly sensitive to the energy of the neutrons, nor their production mechanisms. The observables listed above, however, allow a means to verify the reliability of the Monte Carlo implementations of neutron propagation and capture, in the context of measuring the neutron production rate.

D. Neutron yield
The "neutron yield" is defined as the production rate of neutrons per unit muon track length per unit material density. Here we measure yields in heavy water, both pure and with the NaCl loaded at 0.2% by weight. We define the track length of each muon as µ = 2 R 2 AV − b 2 through the target volume of density ρ, where R AV = 600 cm is radius of the AV, and N (µ) n to be the number of neutrons produced by the muon. The yield is then where N µ is the total number of muons and avg is the average muon track length.
In principle, the number of neutrons can be determined by simply counting neutron-like events following muons, with the following corrections: we express the probability for a neutron produced by a muon of impact parameter b to be captured in the fiducial volume, the "capture efficiency", as ε Cap (b); and the probability for a neutron capture at radius r to trigger the detector and survive the event selection cuts, the "observation efficiency", as ε Obs (r). With a background count of N (µ) bkg , the number of produced neutrons is then where N (µ) f is the number of follower events, and we account for the relevant efficiencies on a per-neutron and per-muon basis. The number of background counts is comprised of neutrons originating external to the inner volume, radioactive backgrounds coincident with the follower selection window, and radioisotopic backgrounds also produced in spallation reactions, respectively. Estimates for the number of background counts in both phases are given in Section IV G. The first expression in Eq. (2) is an idealized production rate, measured under the assumption that neutron production is a Poisson process, occurring constantly along the path of the muon. This is largely untrue, however, as the majority of production actually occurs during showering [32]. The Poisson rate is equal to the mean permuon yield were each muon to have equal track length. This is, in general, distinct from the mean of the true per-muon yield values calculated using the track length appropriate to each muon, which we denote byȲ n . Because SNO is able to reliably reconstruct individual muon tracks, we calculate a per-muon yield unique to each muon, and computeȲ n as the mean Y (µ) n .

E. Capture efficiency
The capture efficiency is defined as the fraction of neutrons produced by a muon that are captured in the fiducial volume, parameterized as a function of the impact parameter of the muon. A 252 Cf source was deployed in SNO to measure the capture efficiency of MeV-scale neutrons (see Figure 13), but the energy spectrum from cosmogenic production extends much higher, and the capture efficiency in this regime may be different. We thus evaluate this efficiency solely using GEANT4 simulations. An uncertainty on the capture efficiency due to the spectrum of starting neutron energies, shown in Figure 2, is calculated by computing the efficiency in ten bins in energy, ranging from 0 to 5 GeV, and computing the RMS difference of these binned efficiencies from the nominal value, weighted by each bin's integral of the energy spectrum. The capture efficiencies in both phases are shown in Figure 3. The cosmogenic capture efficiency curves differ from those measured with the 252 Cf source (see Figure 13) for two reasons: principally, the cosmogenic capture efficiency is parameterized by the muon impact parameter, not neutron starting position, and also differences in the neutron energy spectra.

F. Observation efficiency
The observation efficiency is defined as the probability for a neutron capture through a visible capture mode to trigger the detector and pass the event selection criteria outlined in Section IV B. We evaluate this efficiency by propagating and reconstucting capture gammas in SNO-MAN. This efficiency is shown in Figure 4. Because the energy threshold used in this analysis is lower than that used in past solar neutrino analyses, this efficiency is comparable in both phases, and relatively stable with respect to position in the detector.

G. Backgrounds
The yield measurement as defined in Equations (2) and (3) is subject to three general classes of background, namely cosmogenic neutrons from sources other than the detector volume, radioisotopes produced in conjunction with neutrons, and random coincident events, each of which is discussed below.

External captures
One background to measuring the rate of neutron production in heavy water is contamination from cosmogenic neutrons produced in other materials, which we define as "external captures." At SNO, the principal external sources are the AV and surrounding light water. We assess this contamination as a function of impact parameter, and find, using GEANT4, that the average number of external neutrons capturing in the fiducial volume per muon is at most (5.3 ± 0.2) × 10 −3 in Phase I and (1.5 ± 0.1) × 10 −2 in Phase II, where the larger capture efficiency in Phase II determines the difference.

Cosmogenic radioisotopes
The passage of a muon can result in the production of various unstable isotopes [33], as well as the neutrons that are the focus of this analysis. While the usual concern for cosmogenic production centers on long-lived isotopes, such as 16 N with a half-life of roughly 7 s, the timing cut used to select followers makes this analysis sensitive to the production of short-lived isotopes. From both calculations and measurements of isotope production at Super-Kamiokande [33,34], we determine the expected dominant isotope background to be 12 B, a beta-emitter with a half-life of 20 ms and Q-value of 13 MeV. Our approach to assessing the contribution of this background is data-driven: we search for contamination from 12 B decays using a maximum likelihood fit of both the timing and energy distributions of events following cosmic muons. Explicitly, where t and E are the time delay and energy of each event, we construct a likelihood function where τ 1 = 20 ms/ ln 2 is the 12 B lifetime, and P NC and P B are the reconstructed energy spectra for neutron captures and 12 B β-decays, respectively. The fit parameters are τ , the neutron capture time, and f B , the fractional 12 B contamination. The fit is performed separately on the samples of follower events in each phase; the results of the fit in energy space are shown in Figure 5. The best fit capture time constants are consistent with those fit under the boron-free hypothesis (Section V F). We compute an upper limit on the fractional 12 B contamination at the 90% confidence level by marginalizing over the free time constant. This results in limits on the radioisotopic contamination of 2.4% and 0.67% in Phases I and II, respectively, which are included as uncertainties on the measured neutron yield.

Random coincidences
All remaining backgrounds are uncorrelated with the passage of a muon, and are classified as random coincidences. We assess this class of backgrounds by imposing neutron selection criteria on events in a 3-s time window immediately preceding the trigger time of each muon. Doing so determines the average coincidence rates to be 7.89 × 10 −4 s −1 and 9.73 × 10 −4 s −1 in Phases I and II, respectively, which translate to average numbers of coincident events per muon of 2.4 × 10 −2 and 3.9 × 10 −3 , respectively.

V. STUDY OF EVENT DISTRIBUTIONS
To aid in the development and improvement of physical models, both strictly theoretical and those implemented in simulation packages, we present distributions of observables of cosmogenic neutrons and their relation to their leading muon in the data, and a comparison to model-based predictions. Specifically, we show distributions of the track parameters of muons for which neutron followers were observed, follower multiplicity, the capture positions measured both in the detector and in relation to the leading muon, and the time delay between the muon and follower event. In all cases, the MC has been scaled to the normalization of the data, for easy comparison of the shapes of the distributions.

A. Follower selection
The number of muons that have follower events passing the selection criteria described in Section IV B is shown in Table II. The enhanced proportion of muons after which followers were observed in Phase II reflects the higher capture cross section. Figure 6 shows the distributions of muon impact parameter, both for all muons and only those with followers. The pre-selection distributions agree because the input to the Monte Carlo is sampled from the population of muons observed in the data. The shapes of the post-selection distributions are roughly proportional to the muon track length in the detector. With regard to the zenith angle, the subset of muons with followers is representative of the larger population, and is shown in Figure 7.

B. Follower multiplicity
The distributions of the number of neutron-like events following a muon are shown in Figure 8. Muons with hundreds of followers were observed in each phase; indeed, events of such high multiplicity are reproduced using existing simulation tools. The potential disagreement in the number of high-multiplicity events in Phase II, however, may indicate that some reactions on chlorine are mismodeled. This could be attributed to incorrect cross sections for the dominant, low-multiplicity, neutron-producing processes, i.e. photonuclear and neutron inelastic scattering, or incorrect final-state generation after near-complete nuclear breakup at high energies.
Distinct identification of cosmic muons as showering either electromagnetically or hadronically has been demonstrated by studying the distribution of multiplicities of neutron followers in high energy (> 90 GeV) muon-induced showers in liquid scintillator detectors [35].   When imposing shower selection criteria, the multiplicity distribution analagous to those shown in Figure 8 exhibited two peaks, corresponding to electromagnetic and hadronic showering, with the hadronic case corresponding to larger multiplicities. Our data set includes neutrons of all origins, and the distributions shown in Figure 8 do not exhibit the bimodal topography charac-teristic of such shower separation.  Figure 9 shows the distributions of the radial position of neutron captures in the detector. Because the muon flux is uniform in area and, in aggregate, neutrons are produced uniformly along a track, they are, in aggregate, produced uniformly in the volume of the detector. This is reflected in Phase II, where there is a large capture cross section and the capture position is more strongly correlated with production position. In Phase I, where the effective capture cross section is reduced by 2 orders of magnitude, neutrons are more likely to diffuse out of the fiducial volume; this effect grows as the muon and, hence, neutrons are located closer to the edge of the AV, which has a relatively high hydrogen content, and results in a deficiency of captures in the outer fiducial volume compared to the center. The agreement of the comparison shown in Figure 9 constitutes a partial validation of the propagation of neutrons in the GEANT4 detector model, but is complicated by the finite size of the detector. More ideal tests would use large volumes where boundary effects are suppressed.

D. Capture clustering
The majority of neutron production occurs in electromagnetic showers. The initiation of a shower usually entails a very localized energy deposition by the muon, in contrast to the smaller, constant ionization losses. In the electromagnetic case, this energy deposition has a characteristic profile in the direction of the muon track, which at cosmic-muon energies in light water has a width typically on the order of several meters; see [32] for a discussion.
In an attempt to profile the energy deposition relevant to neutron production, we investigate the clustering of muon-induced neutrons. Specifically, we use the neutron capture positions as proxies for their production positions, which act as proxies for energy deposition. We de-fine a clustering metric, σ Long , as the standard deviation in the coordinate of the followers' capture positions measured longitudinally along the muon track. More specifically, we define r n as the reconstructed position of a neutron capture event, r µ entrance and r µ exit as the positions where the muon enters and exits the PSUP, respectively, and x n as the coordinate of the neutron capture measured along the track; that is, x = 1 and The distributions of this clustering metric in both phases are shown in Figure 10. The shapes of the distributions in the top panel are determined as the sum of χ-distributions; a well-known result states that the variance of n normally distribution samples follows a χ 2distribution for n − 1 degrees of freedom. Indeed, the bottom panel of Figure 10 shows the distributions of clustering metrics for muons broken down by multiplicity -those followed by 2 neutrons, and those followed by greater than 2 neutrons -and shows that the 2-neutron widths follow a falling distribution, unlike the bell-shaped curves shown for multi-neutron events.
The mean capture profile width is (1.28 ± 0.06) m in Phase I, and (1.08±0.04) m in Phase II. If interpreted as a length scale over which energy is deposited into hadronic channels, this is smaller than the expected scale for electromagnetic deposition, which in light water occurs over a range of several meters [32].

E. Lateral capture distance
The distributions of the lateral capture distance from the leading track are shown in Figure 11, which follow an anticipated exponential form. The offset in exponen-tial behavior from 0 is due both to neutrons being produced away from the track, and the distance traveled by the neutrons before thermalizing. The characteristic distances, both in data and simulation, in Phase II are reduced in comparison to Phase I, which is expected on the basis of the larger capture cross section for 35 Cl than that for 2 H. A single muon in Phase I preceeded a follower candidate observed more than 12 m away, an extreme not predicted by the Monte Carlo. The muon did not enter the AV, and traveled only through the surrounding light water.
The data from Phase II exhibit a rather gross difference in shape from the Monte Carlo prediction, a phenomenon not observed in Phase I. Indeed, we believe that this points to a problem with GEANT4's treatment of cosmogenic neutrons. While validations of low energy neutron transport have been performed, opportunities to benchmark models of high energy transport are scarce. It is also possible that the energy spectrum of primary neutrons determined in GEANT4 is incorrect, or that the cross sections for scattering from chlorine at high energy are not valid. No such discrepancy is observed in Phase I because low energy neutrons in deuterium experience appreciable random walks, typically several meters in length, before capturing. Any sub-meter difference in the path length traveled at high energy is masked by the effect of this relatively long random walk. Indeed, using a simple toy MC which samples high-energy transport lengths from the Phase II distributions in Figure 11 and low-energy transport lengths from the distribution of random walk lengths that a neutron may experience in pure D 2 O, the resulting distributions exhibit a similar level of agreement as in Phase I, in which no discrepancy is seen.

F. Time delay
Distributions of the delay between a muon's passage through the detector and its follower captures are shown in Figure 12. The data during each phase may be fit with a pure exponential, yielding maximum-likelihood estimators of the characteristic capture time of 48.5 ± 1.3 ms in Phase I, and 5.29 ± 0.07 ms in Phase II. While muoninduced neutrons may be produced with very high energies, this is in agreement with the previously measured capture time for 252 Cf neutrons in the salt phase of 5.29 ± 0.05 ms [18]. As the thermalization time is small in comparison to the overall capture time, this agreement suggests that the modeling of low-energy neutron transport and capture are valid in the presence of chlorine, further indicating that the source of the discrepancy in lateral capture distance is in the high energy regime.

VI. RESULTS FOR NEUTRON YIELD
The measured neutron yield values in pure heavy water and salt-loaded heavy water are found to be, in units of 10 −4 cm 2 / (g · µ), 7.28 ± 0.09 (stat.) +1.59 −1.12 (syst.) and 7.30 ± 0.07 (stat.) +1.40 −1.02 (syst.), respectively. These are to be compared with the respective values predicted by GEANT4 of 7.01 ± 0.014 (stat.) and 7.29 ± 0.014 (stat.), respectively, though it should be noted that systematic uncertainties on the simulated values may be quite large; see the extensive discussion in [11].
The systematic uncertainties for this measurement are shown in Table III, including uncertainties from the  Monte Carlo-based capture and observation efficiencies, as well as the number of neutron-like background counts coincident with a through-going muon.
The dominant uncertainty is due to the Monte Carlobased capture efficiency. A 252 Cf fission source was deployed in both phases to measure a per-neutron capture efficiency for low energy (< 15 MeV) neutrons as a function of position in the detector [18]. We assess an additional uncertainty on the muon-induced capture efficiency by computing a volume-weighted average of the relative error between the capture efficiency for 252 Cf neutrons as reported by GEANT4 and the results of the calibration campaign, which are shown in Figure 13. While the simulation is able to reproduce the gross features of the low-energy capture efficiency in both phases, the disagreement at high radii, where the efficiency decreases substantially, causes this to be the dominant uncertainty.

A. Evaluation of the Poisson hypothesis
The yield value presented above is the measurement of Y n (see Eq. (2)), which is standard in the literature, and is the value appropriate when describing neutron production as a Poisson process. This can be compared to the mean per-muon yield,Ȳ n (see Eq. (5)), which in units of 10 −4 cm 2 / (g · µ) is 7.62 ± 0.89 (stat.) and 9.32 ± 1.22 (stat.), in Phases I and II, respectively. The two rates are consistent in pure heavy water, but not in Phase II, where the discrepancy is 24.4%. The mean permuon yield is more sensitive to high-multiplicity muons than the idealized rate, and indeed the few muons in the tail of the Phase II distribution shown in Figure 8 are the source of this difference. Monte Carlo sampling indicates that a discrepancy this large is not unusual, and suggests that a Poisson rate, while useful for summarizing a gross production rate, should not be interpreted as a parameter fundamental to neutron production.

B. Comparison to other experiments
While no cosmogenic neutron yield measurements have been published for heavy water, several have been performed using liquid scintillator targets. The nuclear composition of heavy water, abundant with weakly bound deuterons, differs from that of the carbon chains typically found in organic liquid scintillators, and so the results should not be compared directly. Still, the average numbers of nucleons per unit volume are comparable, and so the yields should be of similar scale. Figure 14 shows several yield measurements performed with liquid scintillator targets as a function of average muon energy, and a fit to a scaling law of the form Y n = aE b µ recently performed by the Daya Bay Collaboration [14], with both the LSD [3] and this measurement overlaid. The average muon energy at SNO depth was determined using the parameterization in [15]. It is observed that while cosmogenic neutron production in heavy water occurs on a similar scale to the extrapolation from liquid scintillator measurements, it is enhanced, consistent with the greater average mass number. With the SNO+ experiment cur-rently running in the original SNO cavern with plans to record data with both light water and liquid scintillator targets, it will be possible to perform additional yield measurements at this same site using multiple different materials, to further elucidate the nature of neutron production at such high energies. FIG. 14. Power-law fit for the cosmogenic neutron yield in liquid scintillator, performed by the Daya Bay Collaboration [14], with the SNO Phase I and LSD measurements overlaid. The SNO and LSD measurements are not included in the fit, and the target material used in SNO is different.

VII. CONCLUSIONS
Although the production and propagation of cosmogenic neutrons are modeled in publicly available software, such as GEANT4 [28], these models have not been exhaustively tested, particularly at the depth of SNO, due to the scarcity of experimental data. Extrapolations from more shallow experimental sites are not well understood. SNO offers a unique opportunity to test models at this depth, and in this muon energy regime, as well as to understand this source of background events for other experiments at SNOLAB. Community-standard simulation tools are seen to reproduce many characteristic observables of muon-induced neutrons in the SNO detector. However, some discrepancies indicate that these tools may be improved, particularly in the high energy regime. Using these simulation tools, the cosmogenic neutron yield at a depth of 5890 km.w.e. in heavy water, and heavy water loaded with 0.02% NaCl by mass, is found to be, in units of 10 −4 cm 2 / (g · µ), 7.28 ± 0.09 (stat.) +1.59