Measurement of Muon-induced Neutron Production at the China Jinping Underground Laboratory

Solar, terrestrial, and supernova neutrino experiments are subject to muon-induced radioactive backgrounds. The China Jinping Underground Laboratory (CJPL), with its unique advantage of a 2400 m rock coverage and long distance from nuclear power plants, is ideal for MeV-scale neutrino experiments. Using a 1-ton prototype detector of the Jinping Neutrino Experiment (JNE), we detected 343 high-energy cosmic-ray muons and (7.86$ \pm $3.97) muon-induced neutrons from an 820.28-day dataset at the first phase of CJPL (CJPL-I). Based on the muon-induced neutrons, we measured the corresponding muon-induced neutron yield in a liquid scintillator to be $ (3.44 \pm 1.86_{\rm stat.}\pm 0.76_{\rm syst.})\times 10^{-4}\mu ^{-1}\rm g^{-1}cm^{2} $ at an average muon energy of \SI{340}{GeV}. We provided the first study for such neutron background at CJPL. A global fit including this measurement shows a power-law coefficient of (0.75$ \pm $0.02) for the dependence of the neutron yield at the liquid scintillator on muon energy.


I. INTRODUCTION
Cosmic-ray muons, produced by the interaction of primary high-energy cosmic rays with the Earth's atmosphere, have a powerful penetrating ability. When such muons pass through mountain rocks and detector materials, they create neutrons. The neutrons are the main source of environmental radiations for underground lowbackground experiments. They are produced in two ways. The direct production includes muon-nucleus spal- * Supported in part by the National Natural Science Foundation of China (11620101004, 11475093), the Key Laboratory of Particle & Radiation Imaging (TsinghuaUniversity), the CAS Center for Excellence in Particle Physics (CCEPP), and Guangdong Basic and Applied Basic Research Foundation (2019A1515012216). Portion of this work performed at Brookhaven National Laboratory is supported in part by the United States Department of Energy (DE-SC0012704) E-mail: orv@tsinghua.edu.cn lation, quasielastic scattering, and negative muon capture by the nucleus. The indirect one is from the muoninitiated electromagnetic and hadronic showers. The knowledge of cosmic-ray muon-induced neutrons is crucial for the search of rare events associated with neutrino interactions, neutrino-less double-beta decay, and weakly interacting massive particles (WIMPs) in underground low-background experiments.
Located in Sichuan Province, China, CJPL is one of arXiv:2108.04010v2 [hep-ex] 26 Jun 2022 the deepest underground laboratories in the world [14]. With its unique advantage of a 2400 m rock coverage [15] and ∼1000 km distance from commercial nuclear power plants, CJPL is ideal for rare-event experiments, such as dark matter searches [16,17] and neutrino physics. The proposed Jinping Neutrino experiment (JNE) [18] can serve as an ideal observatory for MeV-scale solar, terrestrial [19,20] and supernova [21] neutrinos. This paper reports on a measurement of the cosmic-ray muon-induced neutron yield in the LS with an average muon energy of 340 GeV [22] using a 1-ton prototype detector of the JNE. The dataset contains 820.28 live days from 31 July 2017 to 27 September 2020. Together with other experiments, this analysis and result contribute to the understanding of the relationship between muon energy and muon-induced neutron yield up to several hundred GeV.

II. PROTOTYPE DETECTOR AND PERFORMANCE
The 1-ton prototype detector was constructed in early 2017 to study new neutrino detection techniques and underground backgrounds [23]. Figure 1 shows a schematic diagram of the prototype detector, which measures 2 m in height and diameter. The outermost part is a 5 cm thick lead wall to reduce the environmental background. The innermost part is a 0.645 m radius acrylic spherical vessel containing a 1 t LS [24,25]. It consists of linear alkylbenzene (LAB, the molecular formula is C 18 H 30 ), 0.07 g/l of 2,5-diphenyloxazole (PPO), and 13 mg/l of the wavelength shifter 1,4-bis(2-methylstyryl)-benzene (bis-MSB). Thirty 8-inch Hamamatsu R5912 photomultiplier tubes (PMT) outside the acrylic vessel are used to detect light signals. The PMT gains are maintained at ∼ 10 7 by adjusting the high voltage power supplies. A pure water layer between the acrylic sphere and the stainless steel tank serves as passive shielding to suppress the environmental radioactive background.
The electronics system collects and processes the light signal detected by the PMTs. It consists of 4 CAEN V1751 flash ADC boards and one CAEN V1495 logical trigger module. The V1751 board has eight 10-bit ADC channels with 1 V dynamic ranges and 1 GHz sampling rates. All the PMT signals are digitized in the V1751 boards and then fed into the V1495 boards for logical judgment. We set the PMT-hit discriminator to 10 mV, and lowered it to 5 mV in phase E after a better understanding of PE charge distribution. A trigger is issued when more than 25 PMTs are fired simultaneously in a 125 ns window. We later lowered the threshold to 10 PMTs in phase G to efficiently detect the events affected by total internal reflections of the acrylic-water boundary. The data acquisition system records 1024 ns-sized voltage waveforms of the fired PMTs. For experimental studies, the trigger conditions changed several times. Table I provides the details to define each phase.
We performed a nitrogen purging and sealing experiment on 14 July 2019 to study the laboratory radon leakage into the detector. Meanwhile, we expected that this operation could also remove the oxygen in the LS to reduce the quenching effect, improving the light yield and consequently the energy scale in photo-electrons (p.e.) per MeV. Before the nitrogen purging and sealing, the light yield was measured to be 4010 photons per MeV [25]. This number later increased to 6445.
The 1-ton prototype detector began collecting data on 31 July 2017. Owing to a laboratory maintenance, the entire experiment eventually stopped running on 27 September 2020. We first determined a list of good runs the were neither for pedestal calibration nor detector maintenance. We discarded apparent noise by examining the event rate, single PMT trigger rate, PMT baseline, and baseline fluctuation. The live time after the data quality check was 820.28 days.

III. CALIBRATION
We performed several PMT gain and energy scale calibrations for the detector based on the PMT dark noise (DN) and the decay products of radioactive isotopes in the LS. No dedicated artificial isotopes were deployed because the use of characteristic gamma ray peaks from the contamination to calibrate the small prototype detector was sufficient.
The thermal electron emission from the PMT photocathode caused the major DN, and they mimicked actual single photoelectron emissions. Since an event time window was 1024 ns wide and the triggered position was between 150 to 600 ns, we selected the DN events from 0 to 150 ns for each PMT. We then fitted the DN charge distribution with a realistic PMT response function [26] to obtain the gain value for each PMT. Figure 2(a) shows an example of the DN charge spectrum.
With the gain value for each PMT of each run, we processed all the data and obtained the total charge for each triggered event. Figure 2(b) shows a charge distribution in one run. We extracted the monoenergetic gamma signals emitted by the radioactive isotopes 40 K and 208 Tl.
By scaling the MC simulations 40 K and 208 Tl gammas to the data, we estimated the light yields and energy scales in each DAQ phase Table I. They increased together owing to oxygen removal by nitrogen purging in July 2019. Three sources of energy scale uncertainty were identified: the inconsistency between the 40 K and 208 Tl scales and their fitting errors (1.5 %), background shape uncertainty as shown in Fig. 2(b) (2.34 %), and the PMT gain variation with time (2.84 %). All three contributions produced a light yield uncertainty of 4 %.

IV. SIMULATION
The muon-induced neutron yield refers to the number of neutrons produced by each muon per unit path length and density of the material. Some neutrons produced in the LS escape without being observed, while some detected neutrons are induced by muons traversing other detector components and the rock outside the detector. Considering the above processes, the corrections of the number of muon-induced neutrons in this analysis were estimated by a study based on Geant4 simulation [27,28]. We determined the average path length of muons through the LS using the Geant4 simulation with  the input of muon angular reconstruction result [22]. In the following, we will discuss several important aspects of the simulation.

A. Muon simulation in mountain
We used the Geant4-based simulation package by Guo et al. [22] to simulate the muon penetration in the mountain rock to predict various underground muon characteristic profiles. This simulation used Geant4's own standard electromagnetic and muon-nucleus processes. The Jinping mountain terrain data were obtained from the NASA SRTM3 dataset [29] with a rock density of 2.8 g/cm 3 [30]. The composition of the rock in the simulation utilized the values from the abundance of elements in Earth's crust (percentage by weight) [31]: oxygen (46.1 %), silicon(28.2 %), aluminum (8.2 %), and iron (5.6 %). Gaisser's formula [32,33] reasonably describes the muon flux at sea level. We adopted the modified version from the Daya Bay Collaboration [34] for this simulation, as it has a better flux description for low energy. It features large zenith angle ranges by parameterizing the kinetic energy E and zenith angle θ of cosmic-ray muons. Figure 3 shows the simulated angular distribution of underground muons.

B. Detector and neutron simulation
We used Geant4 to generate the incident muon sample and neutrons along the muon's track. The details of the simulated neutrons and their interactions with matter inside the detector are described in this section.
The detector simulation package was a Geant4-based framework. We adopted the same underground muon energy spectrum and angular distribution as [22] to estimate the neutron detection efficiency, and we used a set of empirical formulas from [35] to generate the neutrons induced by muons. The induced neutron spectrum as a function of the mean energy of muon E µ in GeV is given by where both E µ and E n are in GeV. The neutron zenith angular cos θ distribution with respect to the muon is ( We assumed that the induced neutron azimuth angular φ distribution relative to the parent muon is uniform. Figure 3 shows the angular distributions of muons and neutrons according to Eq. (2).

A. Principle
In this analysis, we selected muon events to search for the associated neutrons. The neutron yield (Y n ) in the LS can be expressed as where N n,LS is the number of neutrons produced by the muons that traverse the LS target, N µ is the number of these muons, L is the track length of muons in the LS, and ρ is the density, which was measured to be 0.86 g/cm 3 .
Out of N n,LS neutrons, N n,LS neutrons are triggered with the efficiency ε LS : N n,LS includes two contributions: neutrons eventually captured in the LS and those escaping to other detector parts but still producing sufficient photons to become triggered.
Because the detected neutrons could emanate from all the detector components and the rock outside the detector, we define the total number of detected neutrons as, where i is the index of the materials, including the LS, water, iron tank, lead, and rock. N n,LS is given by: where ξ LS is the fraction of detected neutrons generated in the LS.
According to Eqs. (3), (4) and (6), the muon-induced neutron yield is finally expressed as The following sections present the selection criteria for the muon-induced neutrons and the determination of all the quantities shown in Eq. (7).

B. Muon event selection and flux
We required a minimum energy deposit of 98 MeV for each muon event. We cut out the electronics noise and flasher events by requiring the ratio of the maximum p.e. among the PMTs to the total p.e. in each event, denoted by r max , to be less than 0.15. Figure 4 shows a twodimensional distribution of the deposited energy E dep and r max . There were 343 high-energy cosmic-ray muons, i.e., N µ in Eq. (7). Using the same method as Ref. [22] and an enlarged exposure, we updated the muon flux of CJPL-I to (3.61 ± 0.19 stat. ± 0.10 sys. ) × 10 −10 cm −2 s −1 .

C. Neutron event selection
We selected events for the neutron capture on hydrogen following a tagged parent muon. The neutron capture time t distribution follows where τ n = 216 µs is the mean neutron capture time. It was calculated using the LS's elemental composition and the thermal neutron capture cross-sections from [36] at the detector temperature.
To avoid electronics-induced baseline distortions, we define the signal window to start from 20 µs after the muon's passage and end at 1020 µs to cover neutron capture time, shown in Fig. 5 as "data". The "backgroud" in Fig. 5 is scaled from the 1020 µs < ∆T < 5 001 020 µs time sideband window following the parent muons. A 2.2 MeV peak was evident in the energy spectrum, but it did not occur in the background. Three peaks were observed in the background energy spectrum. The two at approximately 0 MeV and 1.5 MeV were due to the trigger condition changes. In phases A-F, the 25-PMT threshold was just below 1.5 MeV. The one at approximately 2.6 MeV was owing to the radioactive 208 Tl background. We evaluated the number of muon-induced neutrons N n with an unbinned maximum likelihood fit incorporating both the neutron-like event energy and time to the parent muon in this analysis.
We describe the energy spectrum using the following probability density function: where ω is the ratio of the background rate to the total event rate in the signal time window 20 µs < ∆T < 1020 µs, f bkg is the background probability, and f sig is the signal probability [37]: in which where µ and σ are the mean and standard deviation of the expected 2.2 MeV neutron-hydrogen capture peak, respectively, α is the signal fraction, λ is the slope of the exponential tail caused by energy leak and edge effect, and erf is the Gaussian error function: erf(E) = 2 √ π E 0 e −t 2 dt. Equation (10) describes the MC energy spectrum very well, and the parameters were determined from the MC energy spectrum. Since the detector ran in several different trigger conditions, the final MC energy spectrum was the average of those weighted by the live time for each DAQ phase.
The time spectrum of neutron capture is described by where N n is the total number of neutron captures associated with the selected signal window, and N bkg is the background rate. Considering the selection of signal time window, we have the probability density function where, and which are the normalized signal and background time spectra, respectively.
We consider that the number of neutron-capture events follows a Poisson distribution with the expected value ν. Since the energy and time of neutron captures are independent, the 2-dimensional likelihood function is where m = 39 is the number of neutron candidates in the signal window.
We divided the DAQ phases into two parts: the high threshold period and the low one, namely DAQ A-F and DAQ G-J. Analyzed separately, the numbers of muoninduced neutrons were observed to be 6.37 ± 3.46 stat. and 1.49 ± 1.93 stat. , respectively. The total number of muoninduced neutrons N n was calculated as 7.86 ± 3.97 stat. , which shows the observed neutrons with a statistical significance at the 2σ level. Figure 5 shows the fit of energy (a) and time (b) spectra of neutron capture against the combined dataset of the two periods.
We estimated the systematic effect on the number of neutrons due to the uncertainties of the parameters obtained from MC simulation and the energy scale to be 0.64 % and 0.81 %. They provided a combined systematic uncertainty on N n of 1.0 %.

D. ε LS , L, and ξ LS
The MC simulation determined the neutron detection efficiency ε LS in Eq. (7) for each DAQ phase. Owing to the small volume of the LS region in this detector and the relatively large mean free path of muon-induced neutrons, ε LS is dependent on the neutrons' distribution positions. We compared the neutrons produced along their parent muons track and those uniformly generated in the LS region and evaluated a 6.65 % contribution to the systematic uncertainty of ε LS . Table II shows the result for each DAQ phase. We obtained the weighted average by the live time for ε LS used in Eq. (7).
The underground muon energy and angular spectra and muon selection criteria determined the distribution of L and its uncertainty. The average track length L of a muon passing the LS was estimated to be 108.01 ± 14.10 cm from the MC simulation.
We calculated ξ LS by simulating the detected number of muon-induced neutrons in each component N n,i defined by Eqs. (3), (4) and (5). Table II shows the ξ LS value in each DAQ phase. We obtained the weighted average of ξ LS for use in Eq. (7) using the live time.
As shown in Table II, ε LS and ξ LS varied with the trigger conditions and light yields of the detector. Similar to ε LS and ξ LS , we defined the detection efficiency ε Other and fraction of detected neutrons ξ Other generated in other components, which include the water, acrylic vessel, iron tank, lead wall of the detector, and rock in the Jinping Mountain. They all can contribute to the detected neutrons. ε Other and ξ Other are used to understand the contribution of detected neutrons from other components and the detector performance of neutron detection. Table III lists the ε Other and ξ Other values for DAQ phase A-D, as an example. The sum of all the ξ values equals 100 %.
Previous studies [35,[38][39][40] have shown that the yield follows a power law of average muon energy Figure 6 shows the result of a global power-law fit to the experimental measurements solid red line, including the result from this study. The fit yields coefficients a = (4.52 ± 0.55) × 10 −6 µ −1 g −1 cm 2 and b = (0.75 ± 0.02). Compared with the power-law fit from Daya Bay, this fit was consistent within the permitted error range, but it did not contribute significantly to the global fit owing to a large uncertainty.
We noted that the MC simulation [35] underestimates the data from LSD [11] at a high energy. Mei and Hime [40] attributed the deviation to the lack of experimental data for high-energy muon interactions with nuclei. However, the cross-section of muons calculated using the Bezrukov-Buagaev model [41] and used in the simulation [35] was consistent with the measurement by MACRO [42] and ATLAS [43]. We expect the next hundred-ton JNE detector under construction to increase the exposure by at least 100 times, which will measure   [3], and the black dashed line shows FLUKA-based predictions for the dependence of the neutron yield on muon energy from Ref. [35]. the neutron yield precisely to provide insights to the tension.

VII. SUMMARY
We studied the cosmic-ray data collected by the 1ton prototype detector of JNE at CJPL-I. By analyzing the neutrons associated with the muons, we measured the neutron yield (3.44 ± 1.86 stat. ± 0.76 syst. ) × 10 −4 µ −1 g −1 cm 2 at an average muon energy of 340 GeV, corresponding to an upper limit of Y n < 7.78 × 10 −4 µ −1 g −1 cm 2 at 99% C.L. A power-law fit of the neutron yield against average muon energy shows that our result is consistent with those of previous studies but lacks statistics to contribute significantly. With the next hundred-ton JNE detector under construction, we foresee significantly a much more precise muon-induced neutron study in the near future.