Cosmogenic neutron production at Daya Bay

Neutrons produced by cosmic ray muons are an important background for underground experiments studying neutrino oscillations, neutrinoless double beta decay, dark matter, and other rare-event signals. A measurement of the neutron yield in the three different experimental halls of the Daya Bay Reactor Neutrino Experiment at varying depth is reported. The neutron yield in Daya Bay's liquid scintillator is measured to be $Y_n=(10.26\pm 0.86)\times 10^{-5}$, $(10.22\pm 0.87)\times 10^{-5}$, and $(17.03\pm 1.22)\times 10^{-5}~\mu^{-1}~$g$^{-1}~$cm$^2$ at depths of 250, 265, and 860 meters-water-equivalent. These results are compared to other measurements and the simulated neutron yield in Fluka and Geant4. A global fit including the Daya Bay measurements yields a power law coefficient of $0.77 \pm 0.03$ for the dependence of the neutron yield on muon energy.

Neutrons produced by cosmic ray muons are an important background for underground experiments studying neutrino oscillations, neutrinoless double beta decay, dark matter, and other rareevent signals. A measurement of the neutron yield in the three different experimental halls of the Daya Bay Reactor Neutrino Experiment at varying depth is reported. The neutron yield in Daya Bay's liquid scintillator is measured to be Yn = (10.26 ± 0.86) × 10 −5 , (10.22 ± 0.87) × 10 −5 , and (17.03 ± 1.22) × 10 −5 µ −1 g −1 cm 2 at depths of 250, 265, and 860 meters-water-equivalent. These results are compared to other measurements and the simulated neutron yield in Fluka and Geant4. A global fit including the Daya Bay measurements yields a power law coefficient of 0.77 ± 0.03 for the dependence of the neutron yield on muon energy.

I. INTRODUCTION
Neutrons and other hadrons produced by cosmic ray muons are an important source of background for underground low-background experiments studying neutrino oscillations, double beta decay, dark matter and other rare events. There have been several studies of cosmogenic neutron production. Muon-induced neutron and isotope production has been studied with the CERN Super Proton Synchrotron (SPS) muon beam in 2000 [1]. Studies on neutron and isotope yields in various materials in underground detectors have been performed by the INFN large-volume detector (LVD) [2], Borexino [3], KamLAND [4], and many others [5][6][7][8][9][10][11][12]. This paper reports a measurement of the neutron production rate in liquid scintillator at three different values of average muon energy by the Daya Bay Reactor Neutrino Experiment, an underground low-background neutrino oscillation experiment.
Daya Bay, located near the city of Shenzhen in the Guangdong province in China, is designed to study neutrino oscillations by measuring the survival probability of electron antineutrinos from nuclear reactors [13]. Daya Bay has made increasingly precise measurements of sin 2 2θ 13 [14][15][16][17][18] and the effective neutrino mass-squared difference |∆m 2 ee | [16][17][18]. Figure 1 shows a diagram of the Daya Bay experimental site. The Daya Bay Nuclear Power Plant complex consists of six reactors, producing 17.4 GW of total thermal power. The experiment has three experimental halls (EHs), two halls near the reactors cores (EH1, EH2) and one hall far from the cores (EH3). Relative measurements in multiple detector sites are used to predominantly cancel reactor flux and spectral shape uncertainties. In its full configuration, the experiment employs eight functionally-identical antineutrino detectors (ADs) to decrease detector-related errors, with two placed in each near hall and four in the far hall. The ADs are enclosed in water to shield against backgrounds and located underground to reduce the cosmic ray muon flux. Each site has redundant muon detectors to identify the residual muons. The ADs are designed to identify neutron captures, providing the possibility to identify muon-induced neutrons. The three EHs are at vertical depths of 250, 265, and 860 meters-waterequivalent (m.w.e.) allowing for a measurement of the neutron yield at three different values of average muon energy within the same experiment.

A. Antineutrino Detectors
The Daya Bay ADs [19,20] detect antineutrino interactions via the inverse beta decay (IBD) reaction, ν e + p → e + + n. The ADs, shown in Fig. 2, consist of three concentric cylindrical regions separated by transparent acrylic vessels. The central target region of each AD is 3 m in height and 3 m in diameter and is filled with 20 tons of liquid scintillator loaded with 0.1% gadolinium by weight (GdLS). The second layer, called the gamma catcher, is filled with liquid scintillator (LS) to detect neutron capture gamma rays that escape from the target A map of the layout of the Daya Bay Reactor Neutrino Experiment, including six reactor cores (Daya Bay, Ling Ao I, and Ling Ao II cores) and three experimental halls (two near and one far). The antineutrino detectors (ADs) are located in the underground experimental halls, with two ADs at the Daya Bay Near Hall (EH1), two at the Ling Ao Near Hall (EH2), and four at the Far Hall (EH3).
region. The outer layer is a mineral oil buffer for additional shielding against radioactivity. The stainless steel vessel surrounding the outermost layer is 5 m in height and 5 m in diameter. Each detector contains 192 8-inch photomultiplier tubes (PMTs) distributed uniformly on the inside wall of the containment vessel. Reflectors on the top and bottom improve the light collection and uniformity.
The IBD reaction is characterized by two timecorrelated triggers, the prompt signal coming from the energy loss of the positron in the scintillator and its annihilation, and the delayed signal from the capture of the neutron. The liquid scintillator is loaded with gadolinium (Gd) to increase the capture rate of thermal neutrons which suppresses backgrounds from accidental coincidences. The neutrons are preferentially captured on Gd, producing an 8-MeV gamma ray cascade, which is much higher than the energy range of most background radioactivity.

B. Muon Detectors
The ADs are surrounded by a water shield, which also serves as a water Cherenkov counter providing 4π veto coverage for muons traversing any AD. The water shield is covered by a cosmic-ray detector array made of resistive plate chambers (RPCs). See Ref. [21] for a detailed description of the muon system.
Each AD is shielded from natural radioactive background and cosmogenic neutrons by at least 2.5 m of water in every direction. The water shield is instrumented with 288 (384) PMTs in the near (far) halls to detect muons via Cherenkov radiation. It is optically separated into two individual water Cherenkov detectors, the inner water shield (IWS) and outer water shield (OWS), using a thin layer of diffusely reflecting Tyvek. The OWS is 1 m thick on the sides and bottom of the water pool. A water circulation and purification system is used to maintain water quality and detector performance [22]. Figure 2 shows a diagram of the near site muon detectors.
A system of RPC modules is installed above the water shield. Each module has four RPC layers and has dimensions 2.17 m × 2.20 m × 0.08 m. There are 54 modules in each near hall and 81 modules in the far hall. In addition, two telescope RPC modules positioned 2 m above the main RPC systems in each hall are used for precise muon track reconstruction for a smaller portion of the solid angle to benchmark the muon simulation. The position resolution of the RPCs is approximately 10 cm. The RPCs and telescope RPCs are shown in Fig. 3.

III. SIMULATION
The neutron yield is defined as the number of neutrons produced per muon, per path length of the muon through the material, per density of the material, which in this case is GdLS. Neutrons captured on Gd following an identified AD muon are selected in the data, but corrections to this number are necessary to determine the yield. For example, some neutrons produced by a muon in the GdLS will escape without being detected. In this analysis, the corrections between the produced and detected number of neutrons are derived from a Monte Carlo (MC) simulation. Muon track reconstruction is challenging in Daya Bay due to the detector size and the complicated geometry of the water shield. Instead of determining the average path length of muons through the GdLS based on reconstructed tracks, which would introduce large uncertainties, the average muon path length as determined by simulation is used to calculate the yield. Aspects of the simulation that are important for this analysis are described below.

A. Muon Flux Simulation
The sea level muon flux is well-described by Gaisser's formula [23,24]. For Daya Bay, Gaisser's formula is modified to better describe the muon spectrum at low energies and large zenith angles [25]. A digitized mountain profile is generated based on topographic maps of the Daya Bay Nuclear Power Plant region. The Music code [26,27] is used to propagate muons from the top of the mountain to the underground halls and uses the digitized mountain profile to calculate the path length in rock. Music's standard rock properties (Z = 11, A = 22, and ρ = 2.65 g/cm 3 ) are used in the simulation. Table I shows the underground muon flux for each hall from the Music simulation. Figures 4 and 5 show muon angle and energy distributions at each hall, respectively.  The zenith angle is defined as a muon's angle from the vertical, and azimuth is a muon's horizontal compass angle from true North. Differences in angular distributions at each hall are due to the mountain profile. The error in the total simulated muon flux is about 10%, which includes the uncertainties in the mountain profile mapping, rock composition, density profiling, and Music simulation.

B. Detector Simulation
The muons generated with Music are used as the incident muon sample in the detector simulation to study neutron production by muons. Approximately 4 × 10 9 , 2 × 10 8 , and 2 × 10 9 muons are simulated in EH1, EH2, and EH3, respectively.
The Daya Bay detector MC simulation is based on Geant4 [28,29]. For the purpose of the neutron yield study, neutrons are also simulated in the Daya Bay detectors using Fluka [30,31] as a cross-check of the default Geant4 simulation. The features of both Geant4 and Fluka relevant to the neutron yield analysis are described in this section. Unless otherwise stated, the nominal Geant4-based simulation is used in what follows.

Geant4
Geant4 is a widely-used toolkit for simulating the passage of particles through matter [28,29]. Geant4 version 9.2p01 is used for this analysis. The physics list used in the simulation is QGS BIC, which applies the binary cascade (BIC) model for hadronic interactions at lower energies (between 70 MeV and 9.9 GeV for protons and neutrons). For hadronic interactions at higher energy, a quark-gluon string (QGS) model is applied. Geant4's precompound model is used for hadronic interactions at the lowest energies (below 70 MeV) and as a nuclear de-excitation model within the higher-energy QGS model. For neutron elastic and inelastic interactions below 20 MeV, a data-driven model is used (NEU-TRONHP). For neutron capture on Gd, Geant4's de-fault neutron capture library is modified to require energy conservation.
The full simulation is conducted without optical processes to increase the simulation speed. Without optical photons, reconstruction algorithms based on PMT hits cannot be used to determine the muon's trajectory and energy deposition. Therefore, the muon's path length and deposited energy are taken from the simulated path length and deposited energy.

Fluka
Fluka is another popular tool for simulations of particle transport and interactions with matter. Fluka version 2011.2b is used for this analysis. The hadron-nucleon interaction model in Fluka is based on resonance production and decay below a few GeV and on the dual parton model at higher energies. For hadron-nucleus interactions, a nuclear interaction model called Peanut [38] is used. For neutron interactions below 20 MeV, Fluka uses its own neutron cross section library containing more than 250 different materials.
The Daya Bay geometry is included in Fluka at the same level of detail used in the Geant4 simulation, with the exception that PMTs are not included in the former. Similar to the Geant4 simulation, the muon's path length and deposited energy are taken from the simulated path length and deposited energy.

A. Analysis Strategy
In this analysis, AD triggers following a detected muon are used to study neutrons produced by cosmic muons.
The neutron yield Y n can be expressed as where N n is the number of neutrons produced in association with N µ muons traversing the GdLS target, L avg is the average path length of muons in the GdLS from simulation, and ρ = 0.86 g/cm 3 [20] is the measured density of Daya Bay's GdLS. The following sections explain the details of the selection of muons traversing the target and the selection of neutrons produced by these muons.

B. Data set
The Daya Bay experiment began collecting data on 24 December 2011 with six ADs. Two ADs were located in EH1, one AD in EH2 and three ADs in EH3. In summer 2012, data taking was paused to install two new ADs, one in EH2 and one in EH3. Operation restarted on 19 October 2012. The results presented here are based on 404 days of data acquired with the full configuration of eight ADs and 217 days of data acquired with six ADs.

C. Muon event selection
In the water pool, muons are tagged by the PMT multiplicity, the number of PMTs with a signal above a threshold of 0.25 photoelectrons. The criterion for a water-pool tagged muon is more than 12 PMTs triggered in the inner or outer water pool. Muons passing through an AD are tagged by the amount of energy deposited in the AD. The criterion for an AD-tagged muon is an AD trigger with visible energy of at least 20 MeV. For this analysis, AD-tagged muons that fall within a [-2 µs, 2 µs] time window of a water-pool tagged muon are selected. (The negative time difference is allowed to account for time offsets between detectors.) Figure 6 shows the energy deposited in an AD for selected muon events. The peak around 800 MeV is due to muons with path lengths of 3-4 m, corresponding to the dimensions of the GdLS region, while the peak at low energy is dominated by muons losing energy as Cherenkov light in the mineral oil. With the criteria specified above, a sample of N µ,Obs muons is observed. Because this sample includes muons which have traversed the AD without entering the GdLS region, a purity correction P µ is applied to obtain N µ , where N µ is the number of muons passing through the GdLS region used in Eq. 1. The purity P µ is obtained from simulations as the ratio of the number of muons with non-zero path length in the GdLS that deposit at least 20 MeV in an AD to the total number of muons that deposit at least 20 MeV in an AD. Table II shows the values of N µ,Obs and P µ for each EH. P µ is approximately 62%; the remaining 38% of the muons deposit at least 20 MeV by passing through the LS region only.
Muons that reach the GdLS deposit a minimum of approximately 80 MeV of energy in the LS. Both P µ and the average muon path length in GdLS, L avg , are geometry-related parameters that depend largely on the muon angular distribution, and therefore the values are obtained from the muon simulation. The muon flux, energy, and angular distributions from the Music simulation of the Daya Bay site described in Section III A are input to the detector simulation. Figure 7 shows the distribution of muon path length through the GdLS from simulation. The average of this distribution is used as L avg . To verify the simulated muon angular distributions, the simulated muons are compared to a sample of muons with reconstructed tracks from data. By searching for coincident hits in the RPCs and telescope RPCs (shown in Fig. 3), a sample of muon tracks is obtained for which the muon direction can be precisely reconstructed. Using this method, the zenith angle (θ) and azimuthal angle (φ) distributions of muons for these RPC Telescope Coincidence (RTC) events are obtained. Figure 8 compares the angular distributions for RTC events in data and simulation. The RTC sample is approximately 1-2% of the total muon sample.
Because of the small angular acceptance of the telescope RPCs, the tracks in the RTC sample are a biased selection of muon tracks. However, this sample can be used to correct the simulated total underground muon distribution based on the ratio of the measured and simulated muon distributions in the RTC sample. rection is done in bins of θ and φ. This data-driven or tuned prediction for the total underground muon distribution is used as an input to the detector simulation, and the tuned muon simulation is used to estimate uncertainties in P µ and L avg . The angular distributions of the RTC events for the tuned muon simulation are nearly identical to the data distributions shown in Fig. 8. The muon-related parameters and the associated uncertainties are summarized in Table II. The uncertainty in the product P µ × L avg is evaluated to account for correlations between the parameters. The values from the nominal muon simulation are used as the central parameter values. The difference between the nominal values and the values obtained with the tuned muon simulation is included as an uncertainty ("tuned-nominal" in Table II). The maximum difference between the values calculated using Geant4 and Fluka is also included as an uncertainty ("Geant4-Fluka"). Uncertainty in the measured θ distribution, estimated conservatively at 10%, and the MC φ distribution, estimated conserva-tively at 5%, also introduce uncertainty in the parameter values from the tuned muon simulation ("θ uncertainty", "φ uncertainty"). An additional uncertainty in P µ is assigned to account for the inclusion of particles other than muons that are incorrectly tagged as muons by a 20 MeV energy deposit in the AD ("non-µ"). For example, neutrons produced from showering muons in the rock surrounding the hall or water may also deposit more than 20 MeV in the AD. Optical processes are turned off in the MC simulation, and therefore contributions to the deposited energy of the muon from Cherenkov light in the mineral oil (MO) region are not included in the simulation. Because the deposited energy is used in the calculation of P µ , a systematic uncertainty is assigned to P µ due to this effect ("MO dep-E"). A toy MC simulation is used to calculate P µ when the energy deposited in the mineral oil is included, and the difference from nominal value is used as the uncertainty. The muon deposited energy distribution in Fig. 6 has been corrected for this effect. The statistical uncertainty is also included. The total uncertainty is calculated as the square root of the sum of the squares of the individual uncertainties.

D. Neutron event selection
To determine the number of neutrons produced due to muons passing through the GdLS, neutron captures on Gd following a muon signal are selected. Simulations show that the average kinetic energy of a neutron produced in the GdLS by a muon is around 40 MeV, with a long tail that extends to around 1 GeV. Neutrons travel an average of ∼40 cm before capturing on Gd.
To select the Gd captures, AD triggers are chosen with energy between 6 and 12 MeV occurring in a signal time window, at least 10 µs and no more than 200 µs, after an AD-tagged muon. Triggers occurring <10 µs after a muon are not used due to afterpulsing and ringing in the PMTs following the passage of a muon. This criterion also vetoes decay electrons from stopping muons. Figure 9 compares the neutron multiplicity between data and MC, where the neutron multiplicity is defined as the number of AD triggers after an AD-tagged muon that satisfy the Gd capture criteria. This distribution has been corrected for readout window efficiency and the effect of blocked triggers, discussed later in this section. Studies indicate that cases where multiple muons pass through the detector before nGd capture candidates are rare and can be neglected in this analysis. The selected events in the signal time window include random backgrounds unrelated to the muon passage. These backgrounds are estimated by selecting AD triggers with energy between 6 and 12 MeV occurring long after the muon passage, in a time window between 1010 µs and 1200 µs after the muon. Because the neutron capture time is ∼30 µs, the contribution of neutron captures in this late time window is negligibly small. Given that both cosmic muons and background events are distributed randomly in time, Poisson statistics dictates that the distribution of background events measured in time since the last muon is an exponential with a time dependence on the muon rate, R µ . Therefore, the selected background events in the late time window slightly underestimates the background contribution to the selected events in the signal time window. A parameter α is used to apply a small correction to the number of events in the background window to take this into account. The number of selected neutron captures (N cap ) is given by where N 10−200 µs is the number of selected events in the signal time window, N 1010−1200 µs is the number of selected events in the background time window, and α is the correction factor. The correction factor α is given by where R µ is the measured rate of AD muons (1.21±0.12 Hz, 0.87±0.09 Hz, and 0.056±0.006 Hz in EH1, EH2, and EH3, respectively [21]). The number of selected events in the signal and background time windows and the value of α are shown in Table III. Figures 10  and 11 compare data and MC for distributions of candidate neutron captures. Figure 10 shows the distribution of time between the muon and delayed events for both data and MC. Figure 11 shows the energy of the delayed events. The background subtraction has been applied in both figures. A systematic uncertainty in the number of selected neutron captures (N cap ) is assigned due to blocked triggers. When the event rate is high, the electronics buffer can become saturated, and any trigger that occurs during this time will be blocked. Blocked triggers can occur when a muon suffers a large energy loss (> 4 GeV) in the AD, which generates many triggers, including neutrons. However, there is no way to determine if the blocked triggers are neutron captures. To be conservative, all blocked triggers are assumed to be neutron captures, and the systematic uncertainty is calculated from data as the number of selected neutrons that have blocked triggers in the time window since last muon relative to the total number of selected neutrons. Table IV gives the uncertainty in N cap due to this effect in each experimental hall. The number of selected neutron captures from Eq. 3 is further corrected for the efficiency of the selection criteria, including the energy selection efficiency ε E , the time window selection efficiency ε t , the fraction of neutrons captured on Gd ε Gd , and the electronics readout window efficiency ε ro . The values and uncertainties for ε E and ε Gd are taken from other Daya Bay analyses [39,40] and are the same for every AD. Because the time window to select neutron captures is different in this analysis, the value of ε t is calculated from the simulation of cosmogenically-produced neutrons. The uncertainty in the value of ε t is estimated by comparing the IBD neutron capture time distribution in data and MC. The electronics readout window efficiency ε ro corrects for the fact that within the 1.2-µs electronics readout window, only the first neutron capture will be read out by the electronics. For high multiplicity events, any subsequent neutron captures within that 1.2-µs readout time will be lost and not counted in the analysis. (The neutron multiplicity distribution in Fig. 9 has been corrected for this effect and the effect of blocked triggers.) The readout window efficiency factor is calculated from the simulation as the ratio of the number of neutron captures that would be selected (because their captures occur first in the time window) to the total number of neutron captures regardless of their timing. The uncertainty in the efficiency is calculated by comparing values of the neutron yield calculated using different-sized time windows to select the signal and background neutrons following a muon. The maximum fractional difference in the yield from nominal is taken as the relative uncertainty in the readout window efficiency. The values and uncertainties for all of the efficiency corrections are summarized in Tables V and VI. Neutrons that are captured in the stainless steel vessel (SSV) instead of Gd are included in the sample if the emitted gammas enter the LS or GdLS and produce a signal that satisfies the selection cuts. A correction is applied to account for the inclusion of these neutron captures in the sample. The SSV correction, f nSSV , is the ratio of the number of neutrons captured in the SSV to the total number of neutron captures, obtained from simulation. The values of this correction factor are shown in Table VII.
Another source of contamination is neutrons selected in time with a signal that is identified as a muon, but is actually some other type of particle. Secondary particles from a showering muon in the rock or water could deposit 20 MeV or more in the AD, causing the event to be incorrectly tagged as a muon, and the subsequent neutron captures are incorrectly included in the sample. To account for this effect, a correction f non-µ is applied, calculated as the ratio from simulation of the number of neutrons captured on Gd following a"non-µ" (≥ 20 MeV energy deposit, but not a muon) to the number of neutrons captured on Gd following any event tagged as a muon (any ≥ 20 MeV energy deposit).
Table VII summarizes the correction values and their uncertainties for each EH. The maximum difference between the values calculated using Geant4 and Fluka ("Geant4-Fluka") is used for the uncertainty in f nSSV .
The uncertainty in f non-µ is determined using the difference between the values calculated using Geant4 and Fluka, plus a smaller uncertainty due to neutron propagation ("n propagation") in the MC, which will be described below. Small statistical uncertainties in both parameters are also included. The total uncertainty for each parameter is the square root of the sum of the squares of the individual uncertainties. With the efficiency and purity corrections described above, the corrected number of neutron captures on Gd following a GdLS muon is The final step is to determine the number of neutrons produced due to muons passing through the GdLS based on the number of neutron captures on Gd given by Eqn. 5. Two parameters, R spill and R det , are used to determine this relationship. R spill accounts for the net effect of spill-in, where neutrons produced by muons outside the GdLS are detected via Gd capture, and spill-out, where neutrons produced by muons in the GdLS escape before capture on Gd. The value of R spill , obtained from simulation, is the ratio of the number of neutrons produced inside the GdLS to the number of neutrons captured on Gd. In EH1, the spill-out correction is approximately 26%, while the spill-in correction is around 20%, leading to a net correction of 6%. The values for R spill for each EH are shown in Table VIII.
The parameter R det accounts for the finite detector size. Some neutrons are neither produced nor captured inside the GdLS, but are indirectly produced by the passage of a muon through the GdLS. For example, a muon passing through the GdLS emits a gamma, which leaves the detector and subsequently produces a neutron that is not detected in the AD. An arbitrarily large GdLS detector would tag this neutron and associate it to the muon, but it goes undetected in this analysis due to the AD size. Therefore, this correction is necessary for consistency with the definition of neutron yield used in other experiments with different detector geometries. The value of R det , obtained from simulation, is the ratio of the number of neutrons produced inside the GdLS to the number of neutrons produced due to a muon's passage through the GdLS (regardless of the generation point of the neutron).
With these corrections, the number of neutrons produced due to the passage of a GdLS muon can be found with The values for R spill and R det and the associated uncertainties are summarized in Table VIII. The uncertainty in the ratio R spill /R det is evaluated to account for correlations between the parameters. The values R spill , R det , and f non-µ depend on the neutron propagation model in the simulation. Uncertainties in these parameters are estimated by comparing neutron propagation in data and simulation. Neutron captures associated with muon tracks from the RTC sample are selected. The neutron capture position is reconstructed with an uncertainty of approximately 20 cm using the method described in Ref. [18]. Figure 12 shows the distribution of the minimum distance between the reconstructed neutron capture position and muon track in semi-log scale for the EH1 data. A linear fit to the logarithm of the number of counts as a function of the distance performed on the data provides an allowed range for the slope from −1.45 to −1.25 m −1 . The neutron propagation in the simulation is modified by applying a scaling factor to the simulated neutron energy. In the simulation, slopes of −1.45 and −1.25 m −1 for the same distribution are obtained when the neutron energy is scaled by factors of 0.87 and 1.83, respectively. The parameters R spill , R det , and f non-µ are calculated in the simulation with the neutron energy scaling factor set to 0.87 and then set to 1.83. This provides a range in the values of R spill /R det and f non-µ which is used as the uncertainty. The same process is applied to the other EHs. This is the dominant source of uncertainty in the neutron yield. The maximum difference between the values of R spill and R det calculated using Geant4 and Fluka is also included as an uncertainty, as well as the statistical uncertainty.  FIG. 12. Semi-log distribution of the perpendicular distance between neutron capture position and the RTC event muon track for EH1 data. The best linear fit of the logarithm of the number of counts as a function of the distance for the data is shown, in addition to the lines drawn with the upper and lower limit slopes used to determine the neutron energy scaling in the MC.

V. RESULTS
The neutron yield calculated by Eq. 1 at each EH is shown in Table IX. E µ avg is the average muon energy for muons passing through the GdLS at each EH, calculated from the Music simulation. The uncertainty in the average muon energy predicted by Music is about 6%, dominated by uncertainties in the mountain profile and rock density. The neutron yields predicted from Geant4 and Fluka simulations at each EH are also shown, with statistical uncertainties only. TABLE IX. Measured and predicted neutron yield for each EH in units of ×10 −5 µ −1 g −1 cm 2 . The measured value is determined from Eq. 1 using the corrections described in previous sections. The predicted values from Geant4 and Fluka are obtained by counting neutrons produced due to simulated muons passing through the GdLS assuming a realistic muon flux and detector geometry. The average muon energy from the Music simulation in each EH is also given.

EH1
EH2 Comparisons of neutron yield measurements from different experiments are typically shown in terms of the average muon energy, despite the differences in muon energy distributions and angular distributions. Daya Bay's measured values for the neutron yield for all three experimental halls are shown in Fig. 13 as a function of average muon energy. The predicted yields at Daya Bay from Geant4 and Fluka are also shown. These results are compared with other experimental measurements over a wide range of average muon energies. For the references that quote an average depth instead of average muon energy, the depth is converted to an average muon energy based on an average rock density [35].
Previous studies [32][33][34][35] have shown that the yield as a function of average muon energy can be described by a power-law, A power-law fit applied to the measurements shown in Fig. 13 including the three points from Daya Bay yields coefficients a = (4.0 ± 0.6) × 10 −6 µ −1 g −1 cm 2 and b = 0.77 ± 0.03. Not all the references quote an uncertainty in the depth or muon energy. Therefore, zero uncertainty in the average muon energy is assumed for all points included in the global fit. Daya Bay has the unique capability to measure neutron production at three different underground sites with essentially identical detectors. A power-law fit applied to the three data points from Daya Bay (including the uncertainty in the average muon energy in the fit) yields coefficients a = (7.2 ± 3.8) × 10 −6 µ −1 g −1 cm 2 and b = 0.64 ± 0.12.
In Fig. 13, the global fit is compared to Fluka-based studies performed by Wang et al [32] and Kudryavtsey et al [33]. The Fluka predictions from this study and Refs. [32,33] are consistent with the measurement at EH3, but predict fewer neutrons for the shallower depths at EH1 and EH2. Other measurements shown in Fig. 13 show similar or even larger discrepancies with respect to the Fluka-based predictions. Geant4 has been shown to predict up to 30% fewer neutrons than Fluka at muon energies above 100 GeV [34], which is consistent with the MC predictions in this analysis. The points for Daya Bay EH1 and EH2, which differ in energy by less than 1 GeV, are shown in the inset. The predicted yields at Daya Bay from Geant4 and Fluka are also shown. Experimental data is shown from Hertenberger [6], Boehm [8], Aberdeen Tunnel [10], KamLAND [4], LVD [2] with corrections from [35], and Borexino [3]. The solid line shows the power-law fit to the global data set including Daya Bay. The dashed line and dash-dotted lines show Fluka-based predictions for the dependence of the neutron yield on muon energy from Wang et al [32] and Kudryavtsey et al [33].
dependence of the neutron yield on average muon energy is obtained by including the Daya Bay measurements with measurements from other experiments.