Measurements of the neutron background for the Colorado Underground Research Institute (CURIE)

Introduction

Cosmic-ray interactions in the upper atmosphere generate cascades of secondary particles, including pions and kaons [1, 2]. Muons produced primarily through charged-pion and kaon decay constitute the most penetrating and abundant charged component of cosmic rays at sea level [3]. Their flux of approximately $170\ \mathrm{m}^{-2}\mathrm{s}^{-1}$ at sea level and mean energy of approximately 4 GeV [3, 4] produce important backgrounds for experiments searching for rare processes. Rock overburden provides the most effective passive suppression of this component. As muons propagate from the surface to laboratory depth, they lose energy through ionization and radiative processes, and their surviving flux decreases strongly with depth [3, 5]. Underground shielding is conventionally reported in meters water equivalent (m.w.e.), permitting comparisons among facilities with different rock compositions and overburden geometries. The deepest laboratories host flagship rare-event programs, including direct dark-matter searches [6, 7] and searches for neutrinoless double-beta decay [8], while shallower facilities support detector development, material screening and radioassay [9, 10], copper electroforming [11], and cryogenic and quantum-sensor characterization [12].

Reducing the cosmogenic muon flux does not eliminate it as a background source. Surviving muons can often be identified through their characteristic ionization signals and rejected using coincidence vetoes [13, 14, 15], but they also produce secondary radiation in the material surrounding an underground laboratory. Neutrons are particularly problematic because they can penetrate substantial shielding and interact differently from the $\gamma$ rays for which high-$Z$ shielding is designed, requiring dedicated suppression with hydrogen-rich moderators or neutron-absorbing materials. Their signatures can enter directly into experimental regions of interest, e.g., $\gamma$ rays following inelastic scattering, $(n,n'\gamma)$, can deposit energy near a $0\nu\beta\beta$ endpoint [16], neutron capture can activate detector materials and target isotopes [17, 18], and elastic scattering can produce nuclear recoils that mimic dark-matter interactions [19, 20]. The underground neutron field has two physically distinct origins with different energy scales and dependencies on site conditions, requiring both components to be characterized.

The first origin is cosmogenic. Residual muons produce neutrons directly, through spallation on nuclei and $\mu$-capture, and indirectly, through photons and hadrons generated in the electromagnetic and hadronic cascades they induce along their track [3]. The resulting spectrum is hard, extending to GeV scale [21], and these neutrons can traverse substantial shielding before reaching a detector. They are also among the least constrained of all underground backgrounds. Predictions rely on muon photonuclear cross sections carrying large uncertainties [22, 23, 24], and general-purpose transport codes are known both to disagree with one another and to fall below measured yields, in some cases by factors of two or more [25, 26, 27]. Direct measurements are scarce at any depth and scarcer still at shallow sites, so each additional measurement serves the dual purpose of quantifying the background at a given facility and constraining the codes used to predict it elsewhere.

The second origin is radiogenic. The ${}^{238}\mathrm{U}$ and ${}^{232}\mathrm{Th}$ decay chains in the surrounding rock emit $\alpha$ particles of 4 MeV to 9 MeV, which produce neutrons through $(\alpha,n)$ reactions on light nuclei; the spontaneous fission of ${}^{238}\mathrm{U}$ contributes a smaller additional yield [28, 29, 30]. This spectrum is soft, and after moderation in the host rock it is confined to energies below $\approx10$ MeV, where it substantially exceeds the muon-induced contribution at depths beyond approximately 100 m.w.e. [31, 32]. It is this low-energy component that sets passive shielding and moderator requirements. Its magnitude is governed by the local uranium and thorium content and by the abundance of light-$Z$ target nuclei, and measured fluxes vary by an order of magnitude between facilities of comparable depth [33, 34, 35, 36]. Therefore, the radiogenic neutron flux cannot be reliably predicted by depth-based parameterizations, and it must instead be measured in situ at each site.

The Colorado Underground Research Institute (CURIE) is a shallow-underground laboratory operated by the Colorado School of Mines within the Edgar Experimental Mine in Idaho Springs, Colorado [37]. Excavated in stable gneiss and reached by a horizontal, rail-driven tunnel, the facility sits beneath an effective overburden of $\approx400$ m.w.e., attenuating the surface muon flux by a factor of 600–700 [38]. Its primary research spaces are the Subatomic Particle Hideout (SPH) and Cryolab I. The cosmogenic muon flux at CURIE has been measured and the underlying simulation framework validated in Ref. [38], and the muon-induced secondary backgrounds have been predicted through coupled MUTE [39, 40] and GEANT4 [41] simulations in Ref. [42]; the reader is referred to those works for a detailed description of the facility and its cosmogenic environment.

In this work, we present the first direct measurements of the neutron background at CURIE, performed in the SPH. The radiogenic component below approximately 10 MeV was measured with a custom Bonner Cylinder Spectrometer comprising a ${}^{3}\mathrm{He}$ proportional counter and cylindrical high-density polyethylene moderators, and was unfolded into an energy spectrum and integrated flux. The high-energy muon-induced component was measured with a BC-501A organic liquid scintillator, which detects fast neutrons through nuclear recoils and separates them from the underground $\gamma$-ray background using pulse shape discrimination. Together, these measurements characterize the radiogenic and high-energy cosmogenic regimes and provide an experimental benchmark for the simulation framework of Ref. [42].

Radiogenic Neutrons

The radiogenic neutron flux at CURIE was measured in the Subatomic Particle Hideout using a custom Bonner Cylinder Spectrometer (BCS), consisting of a ${}^{3}\mathrm{He}$ proportional counter and six cylindrical high-density polyethylene (HDPE) moderators. The detector response was simulated in GEANT4 and empirically calibrated against a ${}^{239}\mathrm{Pu}$-Be source using the shadow-cone method, described in Section II C.1. Because the calibration measurements used parallel-beam geometries whereas the underground field is approximately isotropic, an angular correction was derived before the measured count rates were unfolded into the neutron spectrum and integrated flux.

Bonner Cylinder Spectrometer

The BCS exploits the 5330 b thermal-neutron capture cross section of ${}^{3}\mathrm{He}$ [43]. The cylindrical proportional counter used here [44] was filled with ${}^{3}\mathrm{He}$ at a pressure of 2 atm and detects neutrons through the reaction

which releases 764 keV of kinetic energy, shared between a 573 keV proton and a 191 keV triton. The resulting pulse-height signature provides efficient neutron identification with low sensitivity to the ambient $\gamma$-ray background.

Sensitivity above thermal energies is conventionally obtained by enclosing the counter in spherical HDPE moderators of increasing diameter [45]. In the low neutron flux at CURIE, however, the counter size permitted by a spherical moderator would require prohibitively long acquisition times. We therefore use cylindrical moderators, which accommodate a larger active volume at the cost of an angularly dependent response. The correction for this dependence is described in Sec. II D.

Six HDPE cylinders of common height 12 in and diameters from 2 to 11 in were selected using GEANT4 response simulations. Together with the bare counter, they provide seven configurations spanning thermal energies to approximately 40 MeV, as shown in Figure 1 and summarized in Table I. Each cylinder was capped with an HDPE plug, which fully covered the exposed ${}^{3}\mathrm{He}$ proportional counter.

Image not extracted

Response Functions

The response function, $R(E)$, relates an incident neutron field to the resulting detector signal, incorporating the detector geometry, material composition, interaction cross sections, and intrinsic detection efficiency. For a monoenergetic neutron fluence, it is defined as

where $N_{\mathrm{d}}$ is the number of detected ${}^{3}\mathrm{He}(n,p){}^{3}\mathrm{H}$ events and $\Phi(E)$ is the incident fluence. For isotropic incidence, this is the local scalar fluence, integrated over all arrival directions. The response therefore has units of area and can be interpreted as an effective response area that combines the geometric acceptance and detection probability at neutron energy $E$. Increasing the moderator thickness shifts the detector sensitivity toward higher neutron energies.

All response functions were simulated in GEANT4 [41] version 11.2.1 using the SHIELDING_HP physics list [46] and G4NDL neutron data. Thermal scattering for hydrogen bound in polyethylene was enabled. The ${}^{3}\mathrm{He}$ counter was modeled according to the manufacturer specifications [44], and the moderators were modeled as HDPE with a density of $0.96\,\mathrm{g\,cm^{-3}}$ and mass fractions of $14.37\%$ hydrogen and $85.63\%$ carbon. Each response function was evaluated at 140 logarithmically spaced neutron energies from $1 \times10^{-9}$ to $40\,\mathrm{MeV}$, with $1 \times10^{8}$ neutrons simulated at each energy.

For a parallel beam, the incident fluence is defined by the number of generated neutrons, $N_{\mathrm{gen}}$, and the area of the generation plane, $A_{\mathrm{gen}}$. (2) then becomes

Configuration	$D \times H \times T$ (in)	$E_{\mathrm{peak}}$ (MeV)
Bare	$1.3 \times8.9 \times-$	Thermal
Cylinder 1 (C1)	$2.0 \times12.0 \times0.35$	$10^{-7}$
Cylinder 2 (C2)	$3.0 \times12.0 \times0.85$	$10^{-6}$
Cylinder 3 (C3)	$5.0 \times12.0 \times1.85$	$1$
Cylinder 4 (C4)	$7.0 \times12.0 \times2.85$	$3$
Cylinder 5 (C5)	$9.0 \times12.0 \times3.85$	$5$
Cylinder 6 (C6)	$11.0 \times12.0 \times4.85$	$10$

For incidence parallel to the cylinder axis, the generation plane is circular, whereas for side-on incidence it is rectangular, with $A_{\mathrm{gen}} = DH$, where $D$ and $H$ are the cylinder diameter and height. At an intermediate incidence angle $\theta_{\mathrm{inc}}$, measured relative to the cylinder axis, the projected area is

This expression recovers the face-on and side-on projected areas at $\theta_{\mathrm{inc}} = 0^\circ$ and $\theta_{\mathrm{inc}} = 90^\circ$, respectively.

Previous cylindrical neutron spectrometers have generally calculated their response functions using parallel-beam simulations [47, 48, 49]. This treatment is appropriate for the source characterization measurements described in Sec. II C, but not for the approximately isotropic neutron field encountered underground. Applying a response calculated for a single incident direction to an isotropic field can therefore introduce a systematic bias [47].

For an isotropic field, no unique generation plane exists. The appropriate geometric quantity is instead the mean projected area presented by the detector over all incident directions, which we denote the effective acceptance area, $A_{\mathrm{eff}}$. For a convex body, Cauchy’s mean-projected-area theorem gives

where $S$ is the total surface area [50].

The effective acceptance areas were also evaluated directly in GEANT4. Each detector configuration was placed at the center of a spherical source surface, with neutron starting positions sampled uniformly over the surface and incident directions sampled from a Lambertian distribution. This produces a spatially uniform, angularly isotropic fluence within the sphere. If $R_{\mathrm{s}}$ is the radius of the source sphere and $N_{\mathrm{hit}}$ is the number of generated neutrons whose trajectories intersect the detector assembly, then

The isotropic response function is consequently

The simulated effective areas, listed in Table II, agree with the corresponding analytic values of $S/4$, confirming the normalization of the isotropic source geometry.

Name	$A_{\mathrm{eff}}$ (cm$^2$)
Bare	65.94
C1	131.74
C2	205.22
C3	367.36
C4	549.78
C5	752.46
C6	975.41

The isotropic response functions retain the general energy dependence of their side-on counterparts but are lower in magnitude, as shown in Figure 2. This difference reflects the angular dependence of the cylindrical geometry and demonstrates the importance of calculating the response for the incident field being measured.

Image not extracted

Detector Characterization

As mentioned earlier, the ${}^{3}\mathrm{He}(n,p){}^{3}\mathrm{H}$ reaction releases kinetic energy shared between a 573 keV proton and a 191 keV triton. When both reaction products stop within the active gas, their full combined energy is recorded. Events in which one product reaches the counter wall before depositing its full energy form the proton and triton shelves below the full-energy peak [43]. These characteristic features are visible in the bare-counter pulse-height spectrum shown in Figure 3.

Image not extracted

The full-energy peak was used to determine the channel-to-energy calibration. Neutron events were then selected using the region of interest (ROI) from channel 245 to channel 1085. The lower boundary corresponds to the 191 keV triton edge, while the upper boundary lies above the 764 keV full-energy peak to ensure its complete inclusion. This ROI retains the ${}^{3}\mathrm{He}(n,p){}^{3}\mathrm{H}$ events while rejecting the lower-energy region dominated by ambient $\gamma$ rays and electronic noise.

Potential background contributions from intrinsic alpha activity in the ${}^{3}\mathrm{He}$ proportional counter were also investigated. A six-hour background measurement was performed in the CoorsTek basement with the bare counter surrounded by 6 in of HDPE to suppress the ambient neutron background. No events were observed within the neutron ROI, indicating no measurable intrinsic alpha contribution under these conditions.

The simulated response functions were characterized experimentally using the shadow-cone method [51, 52]. This method separates neutrons arriving directly from a calibration source from those scattered by the air, floor, walls, and surrounding materials.

Shadow Cone Method

The shadow cone consists of eight 1 in thick steel slabs, with an Fe content greater than 95%, followed by twelve 1 inch thick slabs of 5% borated polyethylene. The assembly approximates a pyramidal frustum with a total length of 20 in, following International Organization for Standardization guidelines for shadow-cone measurements [52, 53]. The smaller, source-facing end measures 2.75 in $\times$ 1.38 in, and the larger, detector-facing end measures 6.23 in $\times$ 5.07 in. The steel section scatters and degrades the energy of neutrons in the direct beam, while the borated polyethylene moderates and absorbs the remaining neutrons.

Two measurements are required for each detector configuration. A measurement without the cone contains both direct and scattered neutrons together with the ambient background. With the cone, the direct beam is blocked, leaving the scattered contribution and ambient background. The ambient background is assumed to remain unchanged between the two measurements and therefore cancels in their subtraction. The number of direct neutron events is therefore

where $F_A(d)$ corrects for neutrons scattered out of the direct beam by the air over the source-to-detector distance $d$. This correction is

where $\overline{\Sigma}$ is the spectrum-averaged macroscopic neutron cross section in air. Subtracting the shadow-cone measurement and applying this correction isolates the direct-beam response needed for comparison with the simulated response functions.

PuBe Neutron Source

The characterization measurements used a cylindrical, doubly encapsulated ${}^{239}\mathrm{Pu}$-${}^{9}\mathrm{Be}$ source (NUMEC-B1, 32 g Pu, 2 Ci) with a reported neutron emission rate of approximately $3.6 \times10^{6}$ n/s. The long ${}^{239}\mathrm{Pu}$ half-life of approximately 24 000 yr makes the neutron emission rate effectively constant over the timescale of the measurements. Neutrons are produced through

The source was modeled in GEANT4 using its reported plutonium activity and the double-encapsulation geometry of Ref. [54]. The simulated neutron emission rate was 2.1% higher than the reported value. The simulated mean neutron energies are $(4.06^{+0.08}_{-0.04})$ MeV over 0 MeV to 10.5 MeV and $(5.15^{+0.04}_{-0.01})$ MeV over 1 MeV to 14.1 MeV, where the uncertainties arise from the energy binning. These mean energies agree with measurements of comparable PuBe sources [55, 56].

Measurements

The characterization measurements were performed in the CoorsTek Atrium at the Colorado School of Mines, a large open space selected to reduce the contribution from room-return neutrons. The detector was initially oriented side-on to the source, as shown in Figure 4. The cone-to-detector distance was fixed at 100 cm following shadow-cone measurement guidelines [53, 57].

Image not extracted

The source-to-cone distance was varied from 4.0 to 39.5 cm using the C4 configuration to identify a plateau region in the measured cone-in count rate [52, 57]. A source-to-cone distance of 19.9 cm was selected within this region to fully shadow the detector while limiting suppression of the scattered field. This distance was used for all seven detector configurations.

Each detector configuration was measured for 300 s with and without the shadow cone. The air-outscatter correction of Eq. (9) was evaluated over the PuBe neutron spectrum, giving $\overline{\Sigma} = 8.35 \times10^{-5}\ \mathrm{cm}^{-1}$. Corresponding GEANT4 simulations were performed using the PuBe spectrum and the same source-detector geometry, with air and room-return contributions omitted, providing only the response to a direct neutron flux. Each simulation corresponded to an equivalent live time of 1200 s.

For each configuration, the ratio of the measured to simulated direct count rates defines an empirical scaling factor,

These factors empirically normalize the simulated response functions to the PuBe measurements [47, 58, 59, 60, 61]. Each factor is applied uniformly over neutron energy, retaining the simulated energy dependence of the response. This procedure assumes that the residual discrepancy can be represented by an overall normalization correction. Because the calibration measures the response integrated over the PuBe spectrum, it does not independently constrain energy-dependent discrepancies. The measured and simulated rates and the resulting side-on scaling factors are listed in Table III.

Name	$R_{\mathrm{meas}}$ (cps)	$R_{\mathrm{sim}}$ (cps)	$S$
Bare	$0.06^{+0.02}_{-0.05}$	$0.07^{+0.01}_{-0.02}$	$0.83^{+0.36}_{-0.70}$
C1	$1.84^{+0.10}_{-0.12}$	$3.66^{+0.16}_{-0.18}$	$0.50^{+0.04}_{-0.04}$
C2	$31.73^{+1.08}_{-1.07}$	$32.09^{+1.02}_{-0.61}$	$0.99^{+0.04}_{-0.05}$
C3	$127.19^{+4.18}_{-4.04}$	$113.11^{+2.86}_{-3.36}$	$1.12^{+0.05}_{-0.04}$
C4	$155.03^{+5.07}_{-5.01}$	$146.95^{+4.71}_{-4.77}$	$1.05^{+0.04}_{-0.05}$
C5	$139.29^{+4.51}_{-4.58}$	$140.13^{+4.99}_{-4.12}$	$0.99^{+0.04}_{-0.05}$
C6	$112.89^{+3.67}_{-3.65}$	$116.69^{+2.56}_{-5.38}$	$0.97^{+0.05}_{-0.04}$

Spectrum Unfolding

Each BCS configuration measures the incident neutron field through a broad, energy-dependent response.

Expressed in terms of count rate and differential neutron flux, the measured rate for configuration $j$ is

where $R_j(E)$ is the response function and $\mathrm{d}\Phi/\mathrm{d}E$ is the differential neutron flux. The equivalent count–fluence formulation is obtained by multiplying the count rate and flux by the measurement live time. Discretizing Eq. (12) over the response-function energy grid gives

where $\phi_i = (\mathrm{d}\Phi/\mathrm{d}E)_i$ and $\Delta E_i$ is the energy bin width.

The seven measured responses cannot uniquely determine the neutron flux in 140 energy bins. Moreover, the response functions are broad and strongly overlapping, so additional moderator configurations would provide only partially independent information. The measurement therefore constrains broad thermal, epithermal, and fast components more reliably than fine bin-to-bin structure. Spectrum unfolding regularizes this underdetermined inversion by imposing physical constraints and incorporating an input spectrum that guides the solution where the data provide limited information.

Three established unfolding algorithms are used in this work: MAXED, GRAVEL, and MLEM. MAXED maximizes the entropy of the solution relative to a supplied default spectrum while constraining the agreement between the predicted and measured count rates through a $\chi^2$ condition [62, 63]. The GRAVEL and maximum-likelihood expectation-maximization (MLEM) algorithms instead iteratively update an initial spectrum until its response-folded count rates converge toward the measurements [64, 65].

The PuBe data provide an internal closure test of the calibrated response matrix and unfolding procedure. The simulated PuBe spectrum described in Sec. II C 2 was normalized and supplied as the initial spectrum on the same 140-bin energy grid as the response functions. Both GRAVEL and MLEM reproduce its principal features, with GRAVEL providing better agreement below approximately 0.1 MeV. Because the measured rates also determine the response scaling factors and the source spectrum supplies the initial shape, this agreement establishes internal consistency but does not independently validate the energy dependence of the response functions. For the underground analysis, MAXED is additionally included to provide a reconstruction based on a methodologically distinct regularization procedure.

Each unfolded spectrum is evaluated by refolding it through the response matrix,

Image not extracted

and comparing the predicted and measured count rates [66, 67]. For the GRAVEL reconstruction, the predicted rates agree with all seven measured rates within their uncertainties, as shown in Figure 6. Integrating the unfolded spectrum gives a source emission rate of

Image not extracted

consistent with the reported PuBe emission rate. The recovered spectral shape, source strength, and response-folded count rates demonstrate closure of the calibrated response and unfolding procedure before it is applied to the underground measurements.

Isotropic Corrections

The scaling factors derived in Sec. II C were measured in a side-on parallel-beam geometry, whereas the underground neutron field is assumed to be approximately isotropic. Ideally, the corrections applied to the isotropic response functions would themselves be determined using a known isotropic neutron field. No calibration source, however, produces a well-characterized and truly isotropic field over the required energy range. This limitation has generally led experiments either to use spherical moderators, whose response is intrinsically near-isotropic, or to approximate the isotropic response using parallel-beam measurements [33, 47, 48, 49, 68].

To account for the angular response of the cylindrical moderators, additional shadow-cone measurements were performed with the detector rotated by 45° and 90° relative to the original side-on orientation (Figure 7). Together with the 0° measurements of Sec. II C, these measurements provide scaling factors $S(0^\circ)$, $S(45^\circ)$, and $S(90^\circ)$ for each configuration, spanning the side-on, intermediate, and face-on geometries.

Image not extracted

The cone-to-detector distance was selected independently for each orientation so that the near face of the cone fully shadowed the projected area of the largest moderator, C6. Distances of 406 cm and 100 cm were used for the 45° and 90° measurements, respectively. The larger distance at 45° was required by the greater projected area in that orientation. These differences do not bias the extracted scaling factors because each measurement was compared with a GEANT4 simulation performed using the same source-detector geometry. Consequently, each ratio

is determined at a matched source-detector distance and can be combined with those from the other orientations.

Each measured scaling factor describes the detector response for a single incident orientation, whereas an isotropic field illuminates the detector from all directions. The correction required for the isotropic response functions

tions is therefore the response-weighted average

where $w(\theta)$ accounts for both the angular phase space of an isotropic field and the detector response at each orientation.

For an isotropic neutron field, the angular weighting includes the solid-angle element $\sin\theta_{\mathrm{inc}}\,\mathrm{d}\theta_{\mathrm{inc}}$, where $\theta_{\mathrm{inc}}$ is measured relative to the cylinder axis. The detector’s angular sensitivity is represented by $R_{\mathrm{sim}}(\theta)$, the simulated count rate integrated over the PuBe neutron spectrum. These rates were calculated using a common source-detector geometry and source normalization, so their angular variation incorporates the projected area and intrinsic detection efficiency averaged over the PuBe spectrum. The common normalization cancels in (17). The angular weight is therefore

The laboratory rotation angle $\theta$ and the incidence angle $\theta_{\mathrm{inc}}$ are complementary, as illustrated in Figure 8. The side-on orientation ($\theta= 0^\circ$) corresponds to incidence perpendicular to the cylinder axis ($\theta_{\mathrm{inc}} = 90^\circ$), while the face-on orientation corresponds to axial incidence. Thus,

Image not extracted

and $\sin\theta_{\mathrm{inc}} = \cos\theta$. This weighting accounts simultaneously for the directional response of the detector and the distribution of arrival directions in an isotropic field.

The detector response was simulated at four additional orientations, 23°, 34°, 56°, and 68°, using a common source-detector geometry. These simulations provide $R_{\mathrm{sim}}(\theta)$, and therefore $w(\theta)$, over the full angular

(end_mid=1 range. The three measured scaling factors were interpolated onto a fine angular grid, and Eq. (17) was evaluated as a discrete sum in 5° intervals.

Because only three measured scaling factors are available for each configuration, the interpolation was performed using both piecewise-linear and quadratic functions. The mean of the two resulting values of $S_{\mathrm{iso}}$ was adopted as the central scaling factor. Rather than using half their separation, the full absolute difference between the two interpolation results was taken as a conservative systematic uncertainty. The statistical uncertainty propagates the Poisson uncertainties in the measured and simulated count rates, while the total systematic uncertainty combines the interpolation uncertainty with that arising from the simulated geometry.

The resulting isotropic scaling factors are listed in Table IV. These uncertainties are propagated into the underground unfolding by scaling the response functions according to the central, statistical, and systematic variations of $S_{\mathrm{iso}}$ and unfolding the resulting response matrices separately.

Configuration	$S_{\mathrm{iso}}$
Bare	$1.000 \pm0.238_{\mathrm{stat.}} \pm0.055_{\mathrm{sys.}}$
C1	$0.897 \pm0.033_{\mathrm{stat.}} \pm0.020_{\mathrm{sys.}}$
C2	$1.028 \pm0.010_{\mathrm{stat.}} \pm0.029_{\mathrm{sys.}}$
C3	$1.147 \pm0.006_{\mathrm{stat.}} \pm0.024_{\mathrm{sys.}}$
C4	$1.019 \pm0.004_{\mathrm{stat.}} \pm0.023_{\mathrm{sys.}}$
C5	$0.917 \pm0.004_{\mathrm{stat.}} \pm0.032_{\mathrm{sys.}}$
C6	$0.835 \pm0.003_{\mathrm{stat.}} \pm0.021_{\mathrm{sys.}}$

Figure 9 compares the isotropic scaling factors obtained using piecewise-linear and quadratic interpolation of the three measured angular scaling factors. The two methods give similar results across all seven configurations. Their full absolute difference is included in the systematic uncertainty to account for sensitivity to the interpolation choice.

Underground Measurements

The BCS measurements were performed in the Subatomic Particle Hideout with the detector deployed near the south wall, approximately 0.5 m above the floor and at least 0.5 m from each wall. The $^{3}\mathrm{He}$ counter was oriented vertically, and the detector assembly was kept at the same position throughout the campaign. The unfolding uses the isotropic response functions, assuming an approximately isotropic local neutron field. Each configuration was measured in a single continuous live-time window, with run lengths chosen to bring the statistical uncertainty on the integrated count rate below approximately 5%.

Image not extracted

The channel-to-energy mapping was fixed by fitting the full-energy thermal-capture peak in a high-statistics PuBe dataset and anchoring it to the known 764 keV $Q$-value. The neutron ROI discussed earlier is applied consistently to every measurement. Its lower bound corresponds to the 191 keV triton edge, with an uncertainty of $\pm5$ channels reflecting the ambiguity in locating the edge. This uncertainty is propagated as a systematic uncertainty on the extracted counts. The upper bound lies several channels above the full-energy peak to ensure its complete inclusion. The resulting ROI retains the $^{3}\mathrm{He}(n,p)^{3}\mathrm{H}$ events while rejecting the lower-energy region dominated by ambient $\gamma$ rays and electronic noise.

During the campaign an unexpected source of electronic noise appeared in the low-energy region of the pulse-height spectrum. To suppress it, the $^{3}\mathrm{He}$ tube and charge-sensitive preamplifier were enclosed in a Faraday cage. The noise was confined to channels below $\approx200$, well separated from the ROI whose lower edge lies at channel 245, and because the noise population lies entirely below the ROI it makes a negligible contribution to the integrated capture counts (Figure 10). The extracted count rates are therefore insensitive to any residual noise not removed by the cage. The measured data are summarized in Table V, where the statistical errors are from counting and the systematic errors are dominated by varying the lower ROI boundary by $\pm5$ channels.

Image not extracted

Name	Counts	$T_{\mathrm{live}}$ (h)	$R_{\mathrm{exp}}$ (cps)
Bare	823	164.44	$(1.390 \pm0.049_{\mathrm{stat.}} \pm0.012_{\mathrm{sys.}}) \times10^{-3}$
C1	1162	242.29	$(1.332 \pm0.039_{\mathrm{stat.}} \pm0.010_{\mathrm{sys.}}) \times10^{-3}$
C2	957	159.22	$(1.670 \pm0.054_{\mathrm{stat.}} \pm0.011_{\mathrm{sys.}}) \times10^{-3}$
C3	740	165.21	$(1.244 \pm0.045_{\mathrm{stat.}} \pm0.017_{\mathrm{sys.}}) \times10^{-3}$
C4	868	264.22	$(9.125 \pm0.315_{\mathrm{stat.}} \pm0.063_{\mathrm{sys.}}) \times10^{-4}$
C5	619	310.60	$(5.536 \pm0.224_{\mathrm{stat.}} \pm0.045_{\mathrm{sys.}}) \times10^{-4}$
C6	693	498.32	$(3.863 \pm0.150_{\mathrm{stat.}} \pm0.061_{\mathrm{sys.}}) \times10^{-4}$

A direct comparison between the underground count-rate distribution (Table V) and that measured with the PuBe source (Table III) already reveals the qualitative character of the underground field, before any unfolding is performed. For the PuBe source, whose mean neutron energy is several MeV, the response peaks at the intermediate-to-large moderators, reaching a maximum at C4 and falling off sharply for the bare and thin configurations that cannot thermalize such a hard spectrum. The underground distribution is markedly differThe response is largest for the bare and thin moderators (Bare–C2) and decreases monotonically toward the largest configurations. This shift of the response peak toward thinner moderators is the signature of a much softer incident spectrum, in which a substantial fraction of neutrons are already thermal or epithermal and are detected efficiently with little or no added moderation, while thick moderators serve only to increase absorption. This is precisely the behavior expected of a radiogenic field dominated by $(\alpha,n)$ and spontaneous-fission neutrons that have been moderated through scattering in the surrounding rock, and it anticipates the thermal- and epithermal-dominated spectrum recovered by the unfolding below.

Input Spectrum

Unlike the PuBe characterization, for which the source spectrum is well known, the underground neutron spectrum is itself the quantity being measured. Its reconstruction therefore requires an assumed input spectrum. Radiogenic neutron production can be estimated from the elemental composition and uranium and thorium activities of the surrounding rock [31, 32, 69]. Although the major-element composition of the host rock was characterized by scanning electron microscopy for the simulations in Ref. [42], the uranium and thorium activities in the rock immediately surrounding the SPH have not been measured. We therefore adopt a physically motivated default spectrum.

The default spectrum contains three components: a Maxwellian distribution for thermal neutrons, a $1/E$ distribution for the epithermal region produced by repeated downscattering, and a Watt spectrum describing neutrons from the spontaneous fission of ${}^{238}\mathrm{U}$ in the surrounding rock [20, 30, 70]. An example is shown in Figure 11, with the three components identified separately. Only the relative spectral shape is supplied to the unfolding algorithms; the absolute normalization is determined by the measured count rates. The input and response functions use the same 140 logarithmically spaced energy bins.

Image not extracted

Dependence on the assumed input was tested using two default spectra constructed from the same Maxwellian, $1/E$, and Watt components, but with substantially different relative weights: one dominated by the thermal component and the other by the fast component. The measured count rates were unfolded with each default spectrum using MAXED, GRAVEL, and MLEM, yielding six reconstructions. Their spread is included in the systematic uncertainty on $\Phi_n$ as an estimate of sensitivity to the selected algorithms and component weightings.

Results

Figure 12 compares the six reconstructions, grouped by default spectrum with the corresponding input overlaid. The reconstructed spectra exhibit consistent broad thermal, epithermal, and fast components despite the different input weightings. This agreement indicates stability within the tested set of algorithms and default spectra. It does not independently establish the detailed shapes of the components, which remain influenced by the input spectra where the response functions provide limited information.

Image not extracted

The final spectrum (Figure 13) is taken as the arithmetic mean of the six reconstructions. Integrating this spectrum gives a total local scalar neutron flux at the detector for $E_n \lesssim10\ \mathrm{MeV}$ of

Image not extracted

The integrated-flux uncertainty includes contributions from the measured counts, response-function scaling factors, and unfolding procedure. The counting contribution is evaluated by varying the input count data by the uncertainties in Table V. For the scaling contribution, the response functions are scaled using the central $S_{\mathrm{iso}}$ values and their upper and lower statistical and systematic bounds. The five resulting response matrices are unfolded separately, and deviations from the nominal result determine the propagated uncertainties. The statistical part originates in the Poisson uncertainties of the measured and simulated calibration counts; the systematic part combines the interpolation spread and simulated-geometry uncertainty (Sec. II D). The unfolding contribution is systematic and is assigned as the maximum deviation of the six reconstructions from their mean. Statistical and systematic contributions are summed separately in quadrature.

As a final check, the averaged spectrum was forward-folded with the response functions to predict the count rate for each configuration. The ratio $C_{\mathrm{pred}}/C_{\mathrm{meas}}$, shown in Figure 14, quantifies the agreement between the reconstructed and measured rates. All seven configurations are consistent with the measurements, with a maximum relative difference of approximately 8 %. The bare configuration has the largest uncertainty because of its low count rate during the PuBe characterization and the resulting statistical uncertainty in its response-function scaling factor. This calibration uncertainty, rather than the counting statistics of the underground measurement, is the dominant statistical limitation associated with the bare response function and propagates into the uncertainty on the integrated flux.

Image not extracted

The integrated flux was also evaluated over three energy intervals: thermal ($E_n \le0.5\ \mathrm{eV}$), epithermal ($0.5\ \mathrm{eV} < E_n \le0.1\ \mathrm{MeV}$), and fast ($E_n > 0.1\ \mathrm{MeV}$). The contribution reconstructed above 10 MeV is negligible relative to the total flux, so no upper bound is imposed on the fast interval. Because the BCS sensitivity decreases at these energies, no detailed interpretation of the spectrum above 10 MeV is attempted.

The integrated flux was also evaluated over three energy intervals: thermal ($E_n \le0.5\ \mathrm{eV}$), epithermal ($0.5\ \mathrm{eV} < E_n \le0.1\ \mathrm{MeV}$), and fast ($E_n > 0.1\ \mathrm{MeV}$). The contribution reconstructed above 10 MeV is negligible relative to the total flux, so no upper bound is imposed on the fast interval. Because the BCS sensitivity decreases at these energies, no detailed interpretation of the spectrum above 10 MeV is attempted.

Table VI compares the resulting integrated fluxes with measurements from other underground facilities. The CURIE values are consistent with previous measurements, particularly those from facilities with similar rock compositions [33, 34]. At CURIE, the epithermal component accounts for approximately 57% of the total flux, consistent with radiogenic neutrons undergoing repeated scattering in the surrounding rock and experimental hall before reaching the detector. The order-of-magnitude variation among facilities with overburdens ranging from 140 to 6230 m.w.e. further demonstrates that radiogenic neutron production depends more strongly on the local uranium and thorium content and host-rock composition than on facility depth.

Muon-induced Neutrons

The BCS measurement of Sec. II constrains the neutron field up to approximately 10 MeV, above which the detector sensitivity decreases substantially. In this higher-energy region, the radiogenic contribution is strongly suppressed and muon-induced neutrons be come the dominant component of the underground neutron background. Below this energy, separating the two populations would require muon (anti)coincidence tagging [27, 76, 77], while the muon-induced contribution remains subdominant to the radiogenic flux [15, 20]. The present measurement therefore focuses on the high-energy neutron field, where the cosmogenic component can be isolated without direct muon tagging.

Underground Facility	Depth (m.w.e.)	Thermal $E_n \le0.5\,\mathrm{eV}$	Epithermal $0.5\,\mathrm{eV} < E_n \le0.1\,\mathrm{MeV}$	Fast $0.1\,\mathrm{MeV} < E_n$	Total $(\mathrm{cm}^{-2}\mathrm{s}^{-1}) \times10^{-6}$
Felsenkeller [71]$^{\mathrm{a}}$	140	—	—	—	$235 \pm43$
CURIE (this work)	406	$9.09 \pm1.55$	$14.40 \pm4.63$	$1.63 \pm0.87$	$25.10 \pm6.99$
JUSL [31]	1604	$34.0 \pm1.3$	$27.2 \pm0.6$	$47.0 \pm1.6$	$107.0 \pm7.0$
Y2L [33]	2000	$14.4 \pm1.5$	—	—	$44.6 \pm6.6$
Soudan [72]	2032	$0.7 \pm0.16$	—	—	—
Yemilab [73]$^{\mathrm{a}}$	2350	$14.1 \pm0.9$	—	$3.1 \pm0.7$	$35.0 \pm1.7$
Canfranc (Hall A) [36]	2450	$3.4 \pm0.2$	—	—	$14.8 \pm0.2$
Kamioka [34]	2700	$7.88$	—	—	$23.5 \pm0.7^{+1.9}_{-2.1}$
LNGS [35]	3100	$2.05 \pm0.06$	$1.28 \pm0.31$	$2.56 \pm0.27$	$5.89 \pm0.64$
SURF [72]	3650	$9.9 \pm1.2$	—	—	—
Modane [74]	4080	$3.57 \pm0.32$	—	—	$9.6$
SURF [72]	4230	$1.7 \pm0.3$	—	—	—
CJPL [75]	6230	$7.03 \pm1.81$	—	—	$26.9 \pm10.2$

Muon-induced neutrons were measured using a Bicron BC-501A organic liquid scintillator (OLS) detector [78]. The cylindrical active volume has a diameter and height of 5 in, corresponding to approximately 1.61 L, and is enclosed in a sealed aluminum housing coupled to a 5 in photomultiplier tube (PMT). Fast neutrons produce nuclear recoils within the scintillator, while differences in the temporal response to nuclear and electronic recoils enable neutron–$\gamma$ discrimination through pulse-shape analysis.

ences in the temporal response to nuclear and electronic recoils enable neutron–$\gamma$ discrimination through pulse-shape analysis. The measured count is converted to an integral neutron flux at the rock–cavern boundary using the simulated transport and detector response. The reported flux follows the first-entry, boundary-area convention of Ref. [42], allowing a comparison with the prediction reported there.

Fast Neutron Interactions

Fast neutrons deposit energy in an organic liquid scintillator primarily by scattering from hydrogen and carbon nuclei [79, 80]. In elastic scattering, part of the neutron kinetic energy is transferred to a recoiling nucleus, which subsequently ionizes the scintillator. A neutron can transfer up to its full kinetic energy to hydrogen in a single collision, compared with approximately 28 % for carbon. Combined with the greater hydrogen abundance and larger elastic-scattering cross section, this makes proton recoils the dominant energy-deposition mechanism at neutron energies below a few MeV.

At higher energies, inelastic reactions on carbon become increasingly important. These include excitation of the 4.44 MeV first excited state of ${}^{12}\mathrm{C}$ and, above approximately 8 MeV, breakup reactions such as ${}^{12}\mathrm{C}(n,n^\prime)3\alpha$ [81, 82]. Energy deposited by the resulting charged particles, including protons, $\alpha$ particles, and the recoiling nuclear fragments, contributes to the scintillation signal and must be included when reconstructing the detector response.

Figure 15 shows the simulated partition of energy deposited in BC-501A by neutrons sampled from the muon-induced spectrum predicted for the SPH in Ref. [42]. Hydrogen recoils account for approximately 80 % of the deposited energy below 10 MeV. The contribution from carbon-inelastic reactions increases above the corresponding reaction thresholds, becomes comparable to that from $n$–$p$ elastic scattering near 80 MeV to 100 MeV, and accounts for approximately two-thirds of the deposited energy at 150 MeV. When integrated over the incident spectrum above 1 MeV, proton recoils contribute 62 % of the deposited energy, while carbon interactions contribute 36 %, including 29 % from inelastic and 7 % from elastic scattering. Although proton recoils dominate the spectrum-integrated response, the increasing contribution from carbon reactions must be included when relating the visible energy to the incident neutron energy.

Image not extracted

Pulse Shape Discrimination

BC-501A is sensitive to both nuclear and electronic recoils. In the underground environment, nuclear recoils are produced predominantly by neutron interactions, while electronic recoils are induced by $\gamma$ rays and residual cosmogenic muons. The greater ionization density of proton and carbon recoils produces a larger slow component in the scintillation light, giving neutron-induced pulses a longer decay tail than electronic-recoil pulses of comparable amplitude [83, 84]. This difference provides the basis for pulse shape discrimination (PSD).

Several PSD methods are available, including charge integration, zero-crossing, rise-time, and frequency-domain techniques [85, 86, 87]. We employ an analog decay-time method. The PMT signal is passed through a delay-line amplifier (DLA), whose amplified unipolar output is sent simultaneously to the digitizer and a pulse shape analyzer (PSA). The PSA generates timing signals when the trailing edge crosses 90 % and 10 % of the pulse amplitude. These signals provide the start and stop inputs to a time-to-amplitude converter (TAC), whose output amplitude is proportional to the interval between the two crossings. A pulse with a longer decay tail therefore produces a larger TAC amplitude, as illustrated in Figure 16. The complete acquisition chain is shown in Figure 17.

Image not extracted

Plotting the TAC amplitude against the DLA pulse amplitude separates events according to their decay time. The surface-level data in Figure 18 exhibit a neutron-like band at larger TAC amplitudes and an electronic-recoil band at smaller TAC amplitudes. The populations are well separated over most of the measured range but begin to overlap at low pulse amplitudes. In this region, the reduced number of detected photoelectrons increases statistical fluctuations in the measured pulse shape, while quenching compresses the visible-energy response of proton and carbon recoils relative to electron recoils [88]. The loss of discrimination at low visible energy ultimately determines the analysis threshold applied to the underground neutron measurement.

Image not extracted

Image not extracted

Energy Calibration

Unlike the $^{3}\mathrm{He}$ measurement, which requires only the identification of events within a fixed ROI, the BC-501A measurement requires the pulse amplitude to be related

to the deposited visible energy. Because Compton scattering dominates the response of the low-$Z$ scintillator to $\gamma$ rays, standard calibration sources produce continua rather than full-energy peaks [89, 90, 91, 92]. The calibration is therefore anchored to the corresponding Compton edges.

Conventional $\gamma$-ray sources provide calibration points only up to approximately 3 MeV. The calibration range for this work was extended using the 4.439 MeV de-excitation $\gamma$ ray from $^{12}\mathrm{C}$ produced by the PuBe source and the 6.128 MeV de-excitation $\gamma$ ray from $^{16}\mathrm{O}$ produced by a $^{244}\mathrm{Cm}$-$^{13}\mathrm{C}$ source. A final high-energy calibration point was obtained from cosmic-ray muons. Unlike the photon sources, muons deposit energy directly in the scintillator and produce a minimum-ionizing peak rather than a Compton continuum. The complete set of calibration energies is listed in Table VII.

Source	$E$ (MeV)	$E_{\mathrm{CE}}$ (MeV)
$^{137}\mathrm{Cs}$	0.662	0.477
$^{40}\mathrm{K}$	1.461	1.243
$^{208}\mathrm{Tl}$	2.614	2.382
$^{12}\mathrm{C}$	4.439	4.197
$^{16}\mathrm{O}$	6.128	5.883
$\mu$	$21.485 \pm0.501$	—

Finite detector resolution broadens each Compton continuum and obscures its otherwise sharp endpoint. To determine the edge consistently, a low-pass smoothing filter was applied to each measured $\gamma$-source spectrum from Table VII and the minimum of the derivative of the smoothed distribution was identified [94]. An example for the 4.439 MeV $\gamma$ ray from the PuBe source is shown in Figure 19.

Image not extracted

The scintillation light output is not strictly proportional to deposited energy across the full dynamic range because the relative contribution of ionization quenching varies with particle stopping power [76, 88]. The response below the high-energy muon region was therefore modeled empirically as a second-order polynomial [76],

constrained by the Compton-edge measurements in Table VII and the muon calibration point. The fit was performed using orthogonal-distance regression [95] so that uncertainties in both the measured edge positions and the associated energies were included. Energies are reported in MeV electron equivalent (MeVee), defined as the energy an electron would deposit to produce the same detector signal as the interacting particle. The muon calibration point was determined by simulating the energy deposited by cosmic-ray muons in the BC-501A detector using geant4, following Ref. [76], and using the energy and angular distributions from Ref. [42]. The simulated distribution has the Landau-like shape expected for ionization energy loss in a finite detector. After accounting for the finite detector resolution, the simulated and measured pulse-height spectra were compared by matching both the minimum-ionizing peak and its high-energy tail, as shown in Figure 20. This procedure gives a most-probable voltage signal of 2.423 V and constrains the high-energy calibration, where the scintillator response is approximately linear.

Image not extracted

The resulting calibration is described by a second-order polynomial below the muon calibration point and

a linear relation at higher energies:

where $a = (0.379 \pm0.162)\ \mathrm{MeVee}/\mathrm{V}^2$, $b = (7.832 \pm0.151)\ \mathrm{MeVee}/\mathrm{V}$, and $c = -0.008 \pm0.025\ \mathrm{MeVee}$. The quoted parameter uncertainties are obtained from the diagonal elements of the fit covariance matrix using the absolute input uncertainties. At the transition voltage $V_{\mu} = 2.423\ \mathrm{V}$, the fitted polynomial gives $E(V_{\mu}) \approx21.191\ \mathrm{MeVee}$. The high-energy linear slope is defined as $E(V_{\mu})/V_{\mu}$, ensuring continuity between the two branches. The simulated muon energy listed in Table VII is a fit input, through which the polynomial is not required to pass exactly. The calibration curve is shown in Figure 21.

Image not extracted

Underground Measurements

The muon-induced neutron measurement was performed in the Subatomic Particle Hideout with the BC-501A detector resting on its side, with its cylinder axis horizontal. The detector was positioned approximately 20 cm above the floor, 50 cm from the south-facing wall, and 140 cm from the laboratory door. Data were acquired for a total live time of $T_{\mathrm{live}} = 1431\ \mathrm{h}$ using a hardware trigger threshold of 300 mV. The waveform of each triggered event was recorded for offline analysis.

The pulse-height spectrum of all recorded events is shown in Figure 22, with the measured voltage converted to visible energy using Eq. (22). The minimum-ionizing muon peak is visible near 21 MeVee. Applying the pulse shape analysis described in Sec. III B separates the data into the populations shown in Figure 23. The dense electronic-recoil band contains events induced primarily by $\gamma$ rays and muons, while the less populated nuclear-recoil band contains the neutron candidates. The two populations are well separated at high visible energy but converge toward the low-energy region, where the discrimination deteriorates.

Image not extracted

Image not extracted

Neutron candidates were selected using an empirical, energy-dependent boundary in TAC amplitude. The boundary was chosen manually to follow the separation between the observed electronic- and nuclear-recoil populations while suppressing electronic-recoil contamination. Its functional form is

where $E$ is the visible energy, $a = 0.70$, $E_0 = 5.75$ MeVee, and $c = 4.37$, with $a$ and $c$ expressed in the TAC-amplitude units used in Figure 23. The exponential term

raises the selection boundary at low visible energy, where the populations approach one another, while the constant term sets the high-energy boundary. Events above this boundary and satisfying the pulse-amplitude requirement are classified as neutron candidates.

A vertical selection boundary at a measured pulse amplitude of 0.577 V excludes the low-energy region where the nuclear- and electronic-recoil populations cannot be reliably separated. Using Eq. (22), this boundary corresponds to a visible-energy threshold of $(4.63 \pm0.17)$ MeVee. The calibration uncertainty is assigned conservatively by varying each coefficient within its quoted uncertainty interval and taking the largest absolute deviation from the nominal calibrated energy over all parameter combinations. Neutron candidates must satisfy both this pulse-amplitude requirement and the PSD boundary. The selection is applied to individual events using their recorded pulse amplitudes, without histogram binning.

The sensitivity of the selected count to the PSD boundary was evaluated by varying $c$ within an offset-variation range $\Delta c = 0.02$, while holding the shape parameters and vertical pulse-amplitude boundary fixed. This variation includes or excludes events near the nominal boundary and gives

neutron candidates over the 1431 h live time. The statistical uncertainty follows from Poisson counting statistics, while the quoted systematic range describes the dependence on the adopted PSD boundary.

PSD Event Corrections

The pulse-shape discrimination selection introduces two potential sources of bias in the measured neutron count: electronic-recoil events may be misidentified as nuclear recoils, while genuine neutron-induced recoil events may be excluded due to the finite neutron acceptance of the PSD selection. Corrections for both effects were therefore studied for the measured neutron flux.

To estimate ER contamination, the DLA pulse-height distribution recorded during the underground measurement was used. The low-energy portion of this spectrum is dominated by the radiogenic ER background, providing an ER-enhanced sample from which the probability of misidentification can be evaluated. Events satisfying a horizontal TAC threshold of 4.37 MeVee were selected. This threshold corresponds to the constant offset of the nominal PSD boundary and conservatively includes events that would otherwise be rejected by the energy-dependent component.

The resulting distribution was fitted with a Crystal Ball function [96, 97] and an additional continuum component. The Crystal Ball function describes the dominant ER peak, while the continuum accounts for events outside this peak, which may include genuine neutron-induced recoils. Both components were included in the total fit, but only the Crystal Ball contribution was used to estimate ER leakage, since the continuum cannot be uniquely attributed to ER events. The fit, shown in Figure 24, was performed by minimizing the Poisson negative log-likelihood, yielding a reduced Poisson deviance of $D/n_{\mathrm{dof}} = 1.1$, where $D$ is the Poisson deviance.

Image not extracted

The expected number of events from the fitted ER peak above the DLA threshold was estimated by integrating the Crystal Ball component,

where $N_{\mathrm{CB}}$ is the fitted peak normalization, $p_{\mathrm{CB}}(E)$ is the normalized Crystal Ball distribution, and $E_{\mathrm{cut}} = 4.63$ MeVee is the DLA analysis threshold.

This yields an estimated ER contribution of $N_{\mathrm{ER}} = 1.01 \pm0.32_{\mathrm{stat.}}{}^{+1.08}_{-0.47}{}_{\mathrm{sys.}}$ events, where the quoted systematic uncertainty was obtained by varying the lower bound of the fit interval, and the statistical uncertainty is calculated from the correlations in the fitting parameters. The corresponding multiplicative correction to the measured neutron count is

with $N_{\mathrm{meas}}$ being the 55 neutron candidates.

To estimate the neutron acceptance of the PSD selection, a second data set was acquired at the surface in a laboratory adjacent to the room housing the plutonium-beryllium source described in Sec. II C 2. The measurement used the same detector and acquisition settings as the underground measurement, including the amplifier gain and acquisition threshold. The proximity to the PuBe source provided a substantially larger population of neutron-induced recoil events than the underground measurement, particularly at low visible energies where the ER and neutron populations approach one another and PSD losses are most likely to occur.

The nominal PSD boundary was applied to this data set, as shown in Figure 25. To suppress ER contamination, an additional horizontal TAC threshold of 4.42 MeVee was imposed, together with the 4.63 MeVee DLA threshold. The TAC threshold was selected to sufficiently suppress the ER population while remaining representative of the low-energy region of the underground data, where the exponential PSD boundary approaches its constant offset of 4.37 MeVee. The upper DLA energy limit, approximately 14.7 MeVee, was defined by the intersection of the horizontal threshold with the nominal PSD boundary. This restricts the analysis to the energy interval in which events passing the horizontal threshold may still be rejected by the nominal PSD selection, thereby defining a neutron-enriched control region.

Image not extracted

Within this region, 199 candidate neutron events were identified, of which six fell below the nominal PSD boundary. The corresponding neutron acceptance is

where the statistical uncertainty is obtained from binomial counting statistics.

The sensitivity to the control-region definition was evaluated using two alternative selections chosen to probe limiting regions of the PSD boundary. The first restricted the analysis to the low-energy interval $4.63 < E_{\mathrm{DLA}} < 10.6$ MeVee, with a more restrictive TAC threshold of 4.48 MeVee, emphasizing the region where the exponential component of the PSD boundary is most significant. The second considered $E_{\mathrm{DLA}}>14.7$ MeVee without an upper energy limit, using a TAC threshold of 4.37 MeVee, corresponding to the constant offset approached by the PSD boundary at high energies. These selections yielded acceptances of $(96.73 \pm1.44_{\mathrm{stat}})\%$ and $(96.70 \pm0.98_{\mathrm{stat}})\%$, respectively. The maximum deviation from the nominal acceptance was 0.28 percentage points, substantially smaller than the statistical uncertainty, and was therefore neglected as an additional systematic contribution.

The resulting correction for neutron events excluded by the PSD selection is

This correction accounts for the finite neutron acceptance of the PSD selection, which is not included in the simulated detector efficiency (Sec. III F) defined by the visible-energy threshold alone.

Detector Efficiency

The BC-501A detection efficiency depends on both neutron energy and arrival direction. An incident neutron must interact in the scintillator and deposit sufficient visible energy to exceed the 4.63 MeVee analysis threshold. The efficiency must therefore be averaged over the energy and angular distributions of the incident neutron field. We determine this flux-weighted efficiency using a dedicated GEANT4 simulation of the detector in the SPH, with neutrons sampled from the muon-induced spectrum and angular distribution predicted in Ref. [42] at the rock–cavern boundary.

The visible energy produced by a nuclear recoil is reduced relative to that produced by an electron with the same deposited energy. The deposited energies of the individual recoil species were therefore converted to MeVee using species-dependent light-output functions for protons, $\alpha$ particles, and carbon, beryllium, and boron recoils [98, 99]:

Here, the deposited recoil energies are expressed in MeV and the corresponding light outputs in MeVee. The visible contributions from all recoil products in an event were summed, while energy deposited through electromagnetic interactions was treated on the electron-equivalent scale.

Near the analysis threshold, the deposited energy is dominated by proton and carbon elastic recoils. Carbon-inelastic reactions, including ${}^{12}\mathrm{C}(n,\alpha){}^{9}\mathrm{Be}$ and ${}^{12}\mathrm{C}(n,p){}^{12}\mathrm{B}$, nevertheless contribute at the few-percent level. Including the light-output response of all relevant recoil species therefore provides a more complete reconstruction of the detector signal.

The operational incident-neutron threshold was defined using simulated events in the lowest visible-energy band, $4.63\le E_{\mathrm{vis}<8.0}$ MeVee. The fifth percentile of the corresponding incident-neutron energy distribution was adopted as the reporting threshold, giving

Approximately 95% of the simulated events in this band therefore originate from neutrons above the reporting threshold. Its systematic uncertainty was evaluated by varying the visible-energy threshold over $(4.63\pm0.17)$ MeVee and varying the percentile criterion from the third to the seventh percentile. The contributions from these two variations were combined in quadrature.

Across all simulated events satisfying $E_{\mathrm{vis}}\ge4.63$ MeVee, a fraction $f_{\le}=0.01574$ originates from neutrons with $E_n<10.5$ MeV. The reporting threshold is therefore an operational definition based on the assumed incident spectrum and detector response, rather than a strict kinematic boundary. This below-threshold contribution is accounted for when converting the measured candidate count into the integral flux above 10.5 MeV, as described in Eq. (40).

Image not extracted

The flux-weighted efficiency was calculated for neutrons with boundary energy $E_n\ge E_{\mathrm{thr}}$ that enter the detector. In each boundary-energy bin, $N_{\mathrm{inc}}(E_k)$ is the number of these neutrons entering the detector and $N_{\mathrm{det}}(E_k)$ is the number producing a reconstructed visible signal at or above 4.63 MeVee. Because the source neutrons are sampled from the predicted boundary spectrum and transported to the detector, the flux-weighted efficiency is

This gives

The quoted uncertainty is systematic. It was evaluated by varying the visible-energy threshold over $4.63 \pm0.17$ MeVee and repeating the simulated mapping to determine the corresponding variation in $E_{\mathrm{thr}}$. The efficiency was then recalculated using each paired visible-energy and incident-neutron threshold. This procedure preserves the correlation between the two quantities and avoids treating the same threshold uncertainty as two independent contributions.

Effective Area and Predicted Rate

Converting the measured neutron count into a flux at the rock–cavern boundary requires the detector acceptance and intrinsic efficiency. These quantities were evaluated using a downstream GEANT4 version 11.2.1 simulation with the SHIELDING_HP physics list and G4NDL neutron data. Each source neutron was initialized with the position, energy, and direction of its corresponding recorded crossing of the rock–cavern boundary in Ref. [42], preserving the correlations among these quantities. The surrounding rock geometry and air-filled SPH cavern were retained, allowing subsequent neutron propagation and scattering within the modeled geometry. A detailed BC-501A detector model, constructed from the measured dimensions and manufacturer-specified scintillator composition [78], was placed at the experimental measurement position. The resulting detector response therefore includes transport from the cavern boundary to the scintillator.

The effective area characterizes the fraction of neutrons entering the laboratory that intersect the detector. It is distinct from the intrinsic efficiency calculated in Sec. III F, which gives the probability that a neutron entering the detector produces a signal above the analysis threshold. Neutrons with $E_n \ge E_{\mathrm{thr}}$ were generated at the laboratory boundary over an area $A_{\mathrm{gen}} = 66.185\ \mathrm{m}^2$, corresponding to the area used to define the neutron flux in Ref. [42]. The effective area relates the boundary flux to the rate of neutrons entering the detector. For source neutrons with boundary energy $E_n \ge E_{\mathrm{thr}}$,

where $N_{\mathrm{gen}}$ is the number of neutrons generated at the laboratory boundary and $N_{\mathrm{hit}}$ is the number entering the BC-501A detector volume. This gives

The uncertainty in $A_{\mathrm{eff}}$ includes contributions from the incident-neutron threshold, the detector dimensions, and the reconstructed detector position. The threshold contribution was evaluated using the correlated visible-energy and incident-neutron threshold variations described in Sec. III F. The detector height and diameter were each varied by $\pm3\ \mathrm{mm}$, and its position within the SPH was varied by $\pm5\ \mathrm{cm}$ in each coordinate. The position dependence was evaluated using all eight combinations of these coordinate displacements, with the largest deviations from the nominal result defining the asymmetric bounds. The position uncertainty dominates the total uncertainty in $A_{\mathrm{eff}}$.

The same detector simulation was used to predict the number of selected events originating from neutrons with boundary energy $E_n \ge10.5\ \mathrm{MeV}$. Neutrons were generated according to the flux predicted in Ref. [42] for an exposure equivalent to the 1431 h experimental live time. Applying the 4.63 MeVee analysis threshold gives

predicted events. The statistical uncertainty is the Poisson uncertainty in the simulated event count, while the systematic uncertainty arises from the normalization uncertainty of the predicted neutron flux.

The measured integral neutron flux at the rock–cavern boundary is obtained from

where $N_{\mathrm{meas}}$ is the observed neutron-candidate count, $T_{\mathrm{live}}$ is the measurement live time, $f_{<}$ is the simulated fraction of selected events originating below the incident-neutron reporting threshold, and $C_{\mathrm{ER}}$ and $C_{\mathrm{NR}}$ are the electronic and nuclear recoil corrections, respectively. The product $A_{\mathrm{eff}}\epsilon$ is the response area calculated for neutrons above that threshold. The factor $(1-f_{<})$ accounts for the below-threshold contribution to the selected sample without changing the reported number of observed candidates.

The visible-energy and incident-neutron thresholds are varied coherently within each detector-response calculation. For the final flux uncertainty, all four limiting combinations of the resulting $A_{\mathrm{eff}}$ and $\epsilon$ values are evaluated, and the largest upward and downward deviations are retained as a deliberately conservative uncertainty. The independent geometry-only contribution is then added in quadrature.

Results

Applying Eq. (40) to the $N_{\mathrm{meas}} = 55 \pm7_{\mathrm{stat}}{}^{+5}_{-8\ \mathrm{sys}}$ neutron candidates recorded over 1431 h gives the integral rock–cavern boundary flux of

The statistical uncertainty follows from the Poisson uncertainty on $N_{\mathrm{meas}}$, as well as from the NR correction factor. The systematic uncertainty combines the variation in the number of candidates selected by the PSD boundary, the geometry-related uncertainty in $A_{\mathrm{eff}}$, the conservative threshold envelope described above, and the ER and NR corrections. These contributions are added in quadrature. The resulting systematic bounds describe a conservative uncertainty range.

The radiogenic contribution to the selected sample was estimated using the unfolded neutron flux above approximately 7 MeV and the BC-501A response over the 1431 h exposure. This gives approximately 0.1 events before any additional rejection by the PSD boundary, corresponding to about 0.2 % of the 55 selected candidates. This estimate is small compared with the statistical uncertainty of the candidate count, and no radiogenic-background subtraction is applied.

The primary comparison of the measurement is with the integral muon-induced neutron flux predicted in Ref. [42],

The measured central value is approximately 10 % higher than the prediction but is consistent within the reported uncertainties. The simulated neutron energy and angular distributions enter the effective area, flux-weighted efficiency, and correction for the below-threshold contribution. This comparison therefore tests the predicted flux normalization under the assumed spectral shape and detector-response model.

Propagating the predicted boundary field through the detector-level simulation and converting the count back to the same boundary-flux convention gives

This value agrees with the integrated prediction to within approximately 1 %. However, the detector-level result is not statistically independent of either the integrated prediction or the measurement: its normalization originates from Ref. [42], while its conversion to flux uses the same effective-area and efficiency calculations applied to the data. It is therefore interpreted as a closure test demonstrating the internal consistency of the detector-response and flux-extraction procedure. The three flux determinations are summarized in Table VIII.

Determination	$\Phi(E_n > 10.5\ \mathrm{MeV})\ (\mathrm{cm}^{-2}\ \mathrm{s}^{-1})$
Integrated prediction [42]	$(1.86 \pm0.29_{\mathrm{sys}}) \times10^{-7}$
Direct simulation	$(1.84 \pm0.27_{\mathrm{stat}}{}^{+0.73}_{-0.46}{}_{\mathrm{sys}}) \times10^{-7}$
Measurement	$(2.09 \pm0.29_{\mathrm{stat}}{}^{+0.79}_{-0.53}{}_{\mathrm{sys}}) \times10^{-7}$

Several independently constrained inputs support this agreement. The MUTE calculation reproduces the measured underground muon flux and its angular dependence at CURIE [38]; the rock composition used in the

GEANT4 model was measured by scanning electron microscopy [42]; and the simulated muon-induced neutron yield reproduces its measured dependence on muon energy after application of the empirical multiplicity correction developed in Ref. [42]. These checks constrain the incident muon field, local material composition, and total neutron yield.

The present measurement provides a single integral constraint on the rock–cavern boundary flux above 10.5 MeV and cannot distinguish a difference in the predicted normalization from a difference in the neutron spectral shape. We therefore retain the simulated differential spectrum from Ref. [42] and adjust its normalization [76] using

Here, $N_{\mathrm{sim}} = 48$ includes only selected events originating from incident neutrons with $E_n \ge10.5\ \mathrm{MeV}$, while $(1-f_{\mathrm{c}})C_{\mathrm{NR}}C_{\mathrm{ER}}N_{\mathrm{meas}}$ is the inferred contribution from the same incident-energy range in the measured sample. The resulting spectrum, shown in Figure 27, retains the simulated boundary spectral shape, with its normalization constrained by the measured integral flux above 10.5 MeV. Below this energy, the scaled spectrum represents an extrapolation of the model.

Image not extracted

Direct measurements of the muon-induced neutron flux remain scarce, particularly near an overburden of 400 m.w.e. Table IX places the present result alongside measurements at other underground facilities, with their respective depths and neutron-energy thresholds. Differences in detector response and energy threshold limit direct quantitative comparisons. The result presented here is, however, qualitatively consistent with these existing measurements, in that the neutron flux decreases approximately with increasing facility depth.

Facility	Depth (m.w.e.)	$\Phi(E_n > E_{\mathrm{thr}})$ (cm$^{-2}$ s$^{-1}$)	$E_{\mathrm{thr}}$ (MeV)
Fermilab [100]	265	$(1.5 \pm1.4) \times10^{-6}$	15
CURIE (this work)	406	$(2.09 \pm0.29_{\mathrm{stat}}{}^{+0.79}_{-0.53}{}_{\mathrm{sys}}) \times10^{-7}$	10.5
BUST [101]	850	$(2.1 \pm0.5) \times10^{-8}$	28.6
Soudan [76]	2032	$(2.23 \pm0.52_{\mathrm{stat}} \pm0.99_{\mathrm{sys}}) \times10^{-9}$	20
Boulby [102]	2860	$\approx4.5 \times10^{-9}$	0.55

Discussion and Conclusions

This work presents the first direct measurements of the radiogenic and muon-induced neutron backgrounds at the Colorado Underground Research Institute. The two components were measured in the Subatomic Particle Hideout using complementary detector systems matched to their respective energy ranges.

The radiogenic component was measured using a custom Bonner Cylinder Spectrometer comprising a ${}^{3}\mathrm{He}$ proportional counter and six cylindrical HDPE moderators. The response functions were calculated in geant4 and calibrated against a ${}^{239}\mathrm{Pu}$-${}^{9}\mathrm{Be}$ source using the shadow-cone method. Refolding and unfolding the PuBe data provided a closure test of the calibrated response matrix. To account for the angular dependence of the cylindrical geometry, measurements at $0^\circ$, $45^\circ$, and $90^\circ$ detector rotation were interpolated and weighted by the simulated angular response and the angular phase space of an isotropic field. This correction permits the larger active volume of a cylindrical spectrometer to be used without leaving its directional response unaccounted for.

Applying the corrected response matrix to the seven underground count rates gives a total low-energy ambient scalar neutron flux at the detector for $E_n \lesssim10$ MeV of

This field is expected to be dominated by radiogenic production; no cosmogenic contribution has been subtracted. The epithermal component accounts for approximately 57 % of the integrated flux. The six reconstructions obtained using three unfolding algorithms and two default spectra give consistent broad thermal, epithermal, and fast components. This consistency supports stability against the tested changes in algorithm and input component weights, while the detailed component shapes remain less constrained. The variation among underground measurements in Table VI, despite facility depths spanning 140 to 6230 m.w.e., underscores the importance of site-specific measurements: radiogenic neutron production depends on the local uranium and thorium content and the composition of the surrounding rock [31, 32].

The muon-induced component was measured with a BC-501A organic liquid scintillator operated for 1431 h. Pulse shape discrimination identified 55 neutron candidates above the 4.63 MeVee analysis threshold, corresponding to an operational neutron energy threshold of 10.5 MeV. After correcting for the effective area, flux-weighted efficiency, the below-threshold contribution, and the electronic and nuclear recoil corrections the measured integral rock–cavern flux is

This central value is approximately 12% higher than the $(1.86 \pm0.29_{\mathrm{sys}}) \times10^{-7}\ \mathrm{cm}^{-2}\mathrm{s}^{-1}$ prediction of Ref. [42] but is consistent within the reported uncertainties. Because the predicted neutron energy and angular distributions enter the calculations of $A_{\mathrm{eff}}$ and $\epsilon$, the measurement tests the normalization of the predicted flux under the simulated response model; it does not independently determine the neutron spectral shape. The detector-level simulation provides an additional closure test but shares inputs and uncertainties with the flux extraction and is therefore not treated as an independent comparison.

The agreement with Ref. [42] is supported by several separately constrained inputs. The MUTE framework reproduces the measured underground muon flux and its angular dependence at CURIE [38]; the rock composition used in the GEANT4 model was determined experimentally; and the simulated dependence of neutron yield on muon energy was tested against existing measurements before applying the empirical multiplicity correction of Ref. [42]. The present result therefore provides an experimental benchmark for the corrected simulation framework at shallow depth. Such measurements remain valuable because general-purpose transport calculations can underpredict muon-induced neutron yields by factors of two or more in some materials and energy regimes [25, 26, 27].

Together, the measurements characterize complementary aspects of the neutron environment at CURIE. The BCS determines the local low-energy scalar flux, which is expected to be dominated by radiogenic production, while the BC-501A measurement constrains the high-energy muon-induced flux at the rock–cavern boundary. Because these quantities use different flux definitions, their numerical values do not directly establish the relative local contributions of the two components. Such a comparison would require the predicted muon-induced field to be evaluated as a local scalar flux at the BCS position. The measurements are nevertheless consistent with the expected transition from a low-energy background governed primarily by local radioactivity to a high-energy background dominated by muon-induced neutrons [20, 31].

The measurements also have important limitations. The radiogenic unfolding is intrinsically underdetermined because the seven response functions are broad and strongly overlapping. Consequently, the measurement constrains the relative strengths of the thermal, epithermal, and fast components more reliably than fine spectral structure. The largest statistical contribution to the radiogenic flux uncertainty arises from the response-function scaling of the bare configuration, whose PuBe calibration rate was low.

function scaling of the bare configuration, whose PuBe calibration rate was low. For the muon-induced measurement, the loss of neutron–$\gamma$ discrimination at low visible energy restricts the analysis to $E_n > 10.5$ MeV. Furthermore, the 55 selected events provide an integral constraint but are insufficient to reconstruct the differential muon-induced spectrum. Resolving that spectrum would require a substantially longer exposure, a larger detector volume, or both.

The methods developed here are applicable to other low-background facilities. The isotropic correction provides a framework for using cylindrical moderator geometries in low-flux environments, where increasing the active counter volume can reduce the required acquisition time. For fast-neutron measurements, the combination of species-dependent light-output reconstruction, a simulation-derived incident-energy threshold, and a flux-weighted efficiency provides a consistent method for converting a PSD-selected event count into an integral neutron flux. Together with the underground muon measurements of Ref. [38] and the muon-induced background calculations of Ref. [42], these results provide a combined characterization of the radiogenic and cosmogenic backgrounds at CURIE and the quantitative inputs needed to design and evaluate future low-background experiments at the facility.

Acknowledgments

This material is based upon work supported by a National Science Foundation Graduate Research Fellowship under Grant No. DGE-213‌7099, a U.S. Department of Energy Office of Science Grant No. DE-FG02-93ER40789, and the Colorado School of Mines via faculty startup funds and the ARCS Foundation.

References