Measurement of muon-induced neutron yield at the China Jinping Underground Laboratory

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 low-background experiments. They are produced in two ways. The direct production includes muon-nucleus spallation, quasielastic scattering, and negative muon capture by the nucleus. The indirect one is from the muon-initiated 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.

In 2000, the CERN Super Proton Synchrotron (SPS) studied the production of the muon-induced neutrons using an artificial muon beam [1]. The cosmic-ray muon-induced neutrons generated in a liquid scintillator (LS) in underground detectors have also been studied by CUBE [2], Daya Bay [3], Aberdeen Tunnel [4], ASD [5–7], KamLAND [8], LVD [9], Borexino [10], LSD [11], and many others [5, 12, 13]. A power-law is used to describe the dependence of a neutron yield in an LS on the average muon energy, but measurements up to a few hundreds of GeV are rare.

Located in Sichuan Province, China, CJPL is one of 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.

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 $ \rm 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 $ \sim 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.

Figure 1.  (color online) Schematic diagram of the 1-ton prototype detector for the Jinping Neutrino Experiment.

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 1 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 that 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.

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.

Figure 2.  (color online) Fit example of DN charge spectrum of a PMT (a) and p.e. charge spectrum in one run (b), in which we used three exponential functions to describe the background in three regions.

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}\rm K $ and $ ^{208}\rm Tl $.

By scaling the MC simulation $ ^{40}\rm K $ and $ ^{208}\rm Tl $ gammas to the data, we estimated the light yields and energy scales in each DAQ phase (Table 1). 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}\rm K $ and $ ^{208}\rm 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%.

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 discuss several important aspects of the simulation.

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³ [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.

Figure 3.  (color online) Angular distribution of the muon and muon-induced neutron. These two figures show the azimuth angular distributions of the muon ($ \phi_{\mu} $) and neutron ($ \phi_{n} $), the zenith angular distribution of the muon ($ \theta_{\mu} $) and neutron ($ \theta_{n} $), and the angular distribution of the neutron relative to the muon ($ \theta_{n}' $) according to Eq. (2). All variables are in the laboratory frame.

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_{\mu} $ in GeV is given by

where both $ E_\mu $ and $ E_n $ are in GeV. The neutron zenith angular $ \cos\theta $ 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).

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, \rm LS} $ is the number of neutrons produced by the muons that traverse the LS target, $ N_{\mu} $ is the number of these muons, L is the average track length of muons in the LS, and ρ is the density, which was measured to be 0.86 $ \rm g/cm^{3} $.

Out of $ N_{n, \rm LS} $ neutrons, $ N'_{n, \rm LS} $ neutrons are triggered with the efficiency $ \varepsilon_{\rm LS} $:

$ N'_{n, \rm 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, \rm LS} $ is given by

where $ \xi_{\rm 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).

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_{\rm max} $, to be less than 0.15. Figure 4 shows a two-dimensional distribution of the deposited energy $ E_{\rm dep} $ and $ r_{\rm max} $. There were 343 high-energy cosmic-ray muons, i.e., $ N_{\mu} $ in Eq. (7). Using the same method as Ref. [22] and an enlarged exposure, we updated the muon flux of CJPL-I to $ \rm (3.61\pm 0.19_{stat.}\pm 0.10_{sys.})\times 10^{-10}cm^{-2}s^{-1} $.

Figure 4.  (color online) Scatter plot of deposited energy and $ r_{\rm max} $; the black dots represent the muon events from the MC simulation, and the red dots represent the muon candidates selected from the experimental data. Three hundred forty-three selected muon events spread in the region of $E_{\rm dep} > 98 \;{\rm MeV}$ and $ r_{\rm max}<0.15 $, while flasher and electronic noise events had larger $ r_{\rm max} $ and relatively smaller $ E_{\rm dep} $ values. We also removed low energy events that may contain indistinguishable radioactive background, shower, and noise events.

We selected events for the neutron capture on hydrogen following a tagged parent muon. The neutron capture time t distribution follows

where $ \tau_{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<\Delta {T}<5001020 $ µ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 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 $ {\rm ^{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.

Figure 5.  (color online) Fit of energy (a) and time (b) spectra of neutron capture. The data points were captured from events in the $ 20<\Delta {T}<1020 \;$μs window after the muon. Background distributions are scaled from time sideband data. The neutron signal energy spectrum was estimated using Monte Carlo detector simulation. The dataset was a simple addition of the low and high threshold periods.

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 <\Delta {T}<1020 $ µs, $ f_{\rm bkg} $ is the background probability, and $ f_{\rm 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 $ \rm erf $ is the Gaussian error function: $ {\rm erf}(E)={2}/{\sqrt{\pi}}\int_{0}^{E} {\rm e}^{-t^{2}}{\rm d}t $. 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-capture events associated with the selected signal window, and $ N_{\rm 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-capture events 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 muon-induced neutrons were observed to be $ 6.37 \pm3.46_{\rm stat.} $ and $ 1.49\pm1.93_{\rm stat.} $, respectively. The total number of muon-induced neutrons $ N'_{n} $ was calculated as $ 7.86\pm3.97_{\rm 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 events 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%.

The MC simulation determined the neutron detection efficiency $ \varepsilon_{\rm 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, $ \varepsilon_{\rm LS} $ is dependent on the neutrons' distribution positions. We compared the neutrons produced along their parent muon's track and those uniformly generated in the LS region and evaluated a 6.65% contribution to the systematic uncertainty of $ \varepsilon_{\rm LS} $. Table 2 shows the result for each DAQ phase. We obtained the weighted average by the live time for $ \varepsilon_{\rm LS} $ used in Eq. (7).

DAQ phase	$\varepsilon_{\rm LS }$(%)	$ \xi_{\rm LS} $(%)
A–D	40.65±2.70	74.01±12.53
E–F	40.74±2.71	72.67±12.24
G	53.54±3.56	67.86±11.25
H	55.17±3.67	61.79±10.09
I, J	56.42±3.75	61.36±9.86

Table 2.  Detection efficiency $ \varepsilon_{\rm LS} $ and fraction factor $ \xi_{\rm LS} $ values in Eq. (7) for different DAQ phases.

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\pm 14.10 \rm \ cm $ from the MC simulation.

We calculated $ \xi_{\rm LS} $ by simulating the detected number of muon-induced neutrons in each component $ N'_{n,i} $ defined by Eqs. (3), (4), and (5). Table 2 shows the $ \xi_{\rm LS} $ value in each DAQ phase. We obtained the weighted average of $ \xi_{\rm LS} $ for use in Eq. (7) using the live time.

As shown in Table 2, $ \varepsilon_{\rm LS} $ and $ \xi_{\rm LS} $ varied with the trigger conditions and light yields of the detector.

Similar to $ \varepsilon_{\rm LS} $ and $ \xi_{\rm LS} $, we defined the detection efficiency $ \varepsilon_{\rm Other} $ and fraction of detected neutrons $ \xi_{\rm 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. $ \varepsilon_{\rm Other} $ and $ \xi_{\rm Other} $ are used to understand the contribution of detected neutrons from other components and the detector performance of neutron detection. Table 3 lists the $ \varepsilon_{\rm Other} $ and $ \xi_{\rm Other} $ values for DAQ phases A–D, as an example. The sum of all the ξ values equals 100%.

	Water	Acrylic	Iron	Lead	Rock
$ \varepsilon_{\rm Other} $(%)	4.89±0.33	18.79±1.25	0.58±0.04	0.50±0.03	0.073±0.005
$ \xi_{\rm Other} $(%)	9.55±2.89	1.00±0.27	2.02±0.47	13.16±2.17	0.26±0.07

Table 3.  Efficiency $ \varepsilon_{_{\rm Other}} $ and fraction factor $ \xi_{_{\rm Other}} $ value of all components for DAQ phases A-D.

Table 4 lists all the variables and uncertainties used in Eq. (7). We calculated the neutron yield in the LS of DAQ phases A–F to be $ Y_{n}=(4.84 \pm 2.63_{\rm stat.}\pm 1.08_{\rm syst.})\times 10^{-4} $ μ⁻¹g⁻¹cm², and DAQ phases G–J to be $ Y_{n}= (2.03 \pm 2.63_{\rm stat.}\pm 0.44_{\rm syst.})\times 10^{-4} $ μ⁻¹g⁻¹cm². These two measurements were statistically independent, but their systematical uncertainties were fully correlated. The combined neutron yield is $ Y_{n}=(3.44 \pm 1.86_{\rm stat.}\pm 0.76_{\rm syst.})\times 10^{-4} $ μ⁻¹g⁻¹cm². Owing to maintenance of the laboratory, the experiment has been shut down since September 2020, preventing an opportunity to accumulate more data to improve the statistical precision as desired. A corresponding upper limit of $ 7.78\times 10^{-4} $ μ⁻¹g⁻¹cm² at 99% C.L. was also indicated.

DAQ phase	Parameter	Value	Statistical uncertainty	Systematic uncertainty
A–F	$ N'_{n} $	6.37	3.46	0.06
	$ N_{\mu} $	256	−	−
	$\varepsilon_{\rm LS } (\%)$	40.68	−	2.70
	$\xi_{\rm LS } (\%)$	73.50	−	12.42
	$\rho /(\rm g/cm^{3})$	0.86	−	−
	$L /\rm cm$	108.01	−	14.10
G–J	$ N'_{n} $	1.49	1.93	0.01
	$ N_{\mu} $	87	−	−
	$\varepsilon_{\rm LS } (\%)$	56.19	−	3.73
	$\xi_{\rm LS} (\%)$	61.81	−	9.96
	$ \rho /(\rm g/cm^{3}) $	0.86	−	−
	$L /\rm cm$	108.01	−	14.10

Table 4.  Parameters used in Eq. (7) to calculate the neutron yield in the LS.

We compare the neutron yield measurements from different experiments in terms of the average muon energy, with various muon energies and angular distributions. Figure 6 shows the measurement at CJPL and the results from CUBE [2], Bezrukov et al. [12], Palo Verde [13], Daya Bay EH1 [3], Daya Bay EH2 [3], Aberdeen Tunnel [4], Enikeev et al. [5], ASD [5–7], Daya Bay EH3 [3], KamLAND [8], LVD2011 [9], Borexino [10], and LSD [11] as a function of average muon energy.

Figure 6.  (color online) Muon-induced neutron yield as a function of average muon energy from Jinping compared with other experiments. The solid red line shows the power-law fit to all the experiments, the blue dashed line is the fit result from Daya Bay [3], and the black dashed line shows FLUKA-based predictions for the dependence of the neutron yield on muon energy from Ref. [35].

Previous studies [35, 38–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\pm0.55)\times 10^{-6} $ μ⁻¹mg⁻¹cm² and $ b= (0.75\pm 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 the neutron yield precisely to provide insights into the tension.

We studied the cosmic-ray data collected by the 1-ton prototype detector of JNE at CJPL-I. By analyzing the neutrons associated with the muons, we measured the neutron yield $ (3.44 \pm 1.86_{\rm stat.}\pm 0.76_{\rm syst.})\times 10^{-4} $ μ⁻¹g⁻¹cm² at an average muon energy of 340 GeV, corresponding to an upper limit of $ Y_{n}<7.78\times 10^{-4} $ μ⁻¹g⁻¹cm² 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 a significantly more precise muon-induced neutron study in the near future.

We thank the anonymous referees for their critical comments and insightful suggestions. We acknowledge Orrin Science Technology, Jingyifan Co., Ltd, and Donchamp Acrylic Co., Ltd, for their efforts in the engineering design and fabrication of the stainless steel and acrylic vessels. Many thanks to the CJPL administration and the Yalong River Hydropower Development Co., Ltd. for logistics and support.