Flux measurement of fast neutrons in the Gran Sasso underground laboratory

The neutron energy spectrum for energies between 5 and 15 MeV has been measured in the INFN Laboratori Nazionali del Gran Sasso in Italy, by using a $1.5~m^3$ stainless-steel tank filled with 1.2 ton Gd-doped liquid scintillator. A two-pulse signature, composed of a prompt and a delayed signal, due to proton recoils and gammas from neutron capture on Gd respectively, has been used for events selection. The neutron energy spectrum has been obtained from the detected proton-recoil spectrum, after unfolding the detector response, which had been simulated with GEANT4. Although this spectrum has been measured previously, the present work represents an improvement in the energy resolution with respect to previously published results. The results obtained with this technique are compared with other measurements existing in literature.


Introduction
The INFN Laboratori Nazionali del Gran Sasso (LNGS) is the largest facility in the world for underground astroparticle physics [1]. It is located in the center of Italy at an average depth of 3600 m w. e., where the cosmic ray flux is suppressed with respect to sea level, and is a suitable place for hosting experiments looking for rare processes, such as neutrinos from gravitational stellar collapse, neutrino-less double beta decay and dark matter interactions. Nevertheless, radiation backgrounds are still present, primarily originating from the natural radioactivity of the surrounding rock. The knowledge of background is therefore essential for a careful design of detectors and their shielding [2]. Neutron transport in matter is mainly computed by simulation packages and both the composition and the geometry of the shielding are calculated on the bases of the impinging neutron flux, known from measurements, and the attenuation factor required to reach certain physics goal.
Aim of the present work is the measurement of the neutron flux inside one of the experimental hall of the LNGS which could be used as input in simulation works for designing future large-scale detectors. Exploiting the existence of a 1.5 m 3 Gd-loaded scintillator detector, partially exposed to the cavern rock radioactivity, a cost effective measurement of the neutron flux was possible. The big volume of the Gd-doped scintillator detector allows the identification of neutron events by their capture on Gd, which is characterized by the emission of a γcascade with enough energy (∼8 MeV) to emerge from background. During the thermalization process preceding its capture, the neutron hits protons of the scintillator and by studying the signals induced by proton recoils the original neutron spectrum can be reconstructed. This paper describes the detector and its location inside the experimental hall A in section 2. Calibration procedures for both, electron and nuclear recoils, are described in section 3. Finally, event selection and unfolding method are discussed in section 4 and 5, respectively. * Email address: gianmarco.bruno@nyu.edu

Experimental Setup
This work has been carried out in the framework of the LVD (Large Volume Detector) experiment, which is a neutrino observatory monitoring the galaxy since June 1992 to detect neutrinos from core collapse supernovae [3]. LVD is placed in the Hall A of the Gran Sasso Laboratory and is made by an array of 840 scintillator detectors 1.5 m 3 each. One of these, highlighted in figure 1, has been adopted to perform the measurement reported in this paper. It is placed on the southern region of the array and has several differences compared with the typical LVD detector, regarding scintillator, photomultipliers, front-end electronics and data acquisition.
The liquid scintillator of the standard LVD detector is a mixture of aliphatic and aromatic hydrocarbons (C n H 2n with n = 9.6), also known as White Spirit or Ragia Minerale. The liquid scintillator of the counter we used for these measurements contains about 8% of aromatics, and is doped with 1 g/l of PPO as activator and 0.03 g/l of POPOP which acts as wavelength shifter. This scintillator has been developed several years ago by the component of the LVD Collaboration from the Institute for Nuclear Research of the Russian Academy of Sciences [4]. Years later, in the frame of a collaboration between the LVD and MetaLS (Metal and Rare Earth doped Liquid Scintillator) groups [5], this counter was doped with [Gd] = 0.93 g/l (about 0.115% in weight). Optical properties and characteristics of this scintillator can be found in [5]. In particular the attenuation length at 425 nm, corresponding to the maximum sensitivity of the photocathode of the used photomultipliers, remains greater than 10 m after Gd doping, well above the detector dimensions (1x1x1.5) m 3 .
Organic scintillators contain a relatively large number of hydrogen atoms per atom of carbon. They are therefore efficient moderators and are potentially suitable for fast neutron spectroscopy as well. On the other hand, doping the scintillator with Gd enhances the capability of tagging neutron captures, with respect to n-captures on protons, resulting in a shorter average capture time and a higher energy emitted by the capture process, both contributing to significantly increase the signal to noise ratio. At a Gd concentration of 0.93 g/l, the detector average thermal n-capture time is ≈ 27µs [5], about 7 times shorter than that of a standard LVD counter [6]. At the same time the total energy emitted by n-capture on Gd is about 8 MeV, instead of the 2.2 MeV γ quantum for the p(n, γ)D reaction.
The detector has been equipped with 3 5" photomultipliers (PMTs) Photonis XP3550, that have a good gain stability over time. The PMTs have been tested and calibrated by measuring the single photoelectron spectrum (SPE) with a light emitting diode (LED) and their gain has been monitored for years. Each PMTs is connected to a 14 bit, 40 MHz bandwidth and 100 MS/s waveform digitizer (CAEN V1724). The trigger condition is the three-fold coincidence among the PMTs. The trigger threshold has been set low enough to have a good efficiency for the detection of gammas from neutron captures and will be discussed in more details in section 4. Finally, to tag muon events, and then to be able to disentangle the muon induced neutron background, the logical OR of all triggers coming from the LVD counters surrounding the detector is also acquired.

Detector Response
The detector response has been studied for both electron and nuclear recoils. The former by using a 60 Co gamma source and the latter by using an AmBe(α, n) neutron source to characterize the light output response to the different incident radiations.

Electron recoils
Energy calibration of the detector has been performed at the beginning of the data taking and repeated regularly during the whole period of data taking of about two years, using a 60 Co source. The 60 Co decays with a half life of t 1/2 = 5.27 y emitting two gammas of energy E 1 = 1.17 and E 2 = 1.33 MeV. The source has been inserted inside the detector using a stainless steel rod and was kept in the center of the scintillator volume during the calibration run. This operation was made under a constant Argon stream to avoid any contact of the scintillator with the oxygen of the atmosphere. The spectrum obtained, after background subtraction, summing up the contribution of all PMTs, is shown in figure 2 black dots.
The mean free path of 1.25 MeV photons in liquid scintillator is about 20 cm, which corresponds to a probability for a gammaray generated in the center of the detector to be completely contained in the active volume of about 95%. The remaining 5% of events contributes to the bump visible in figure 2 from 1 to 1.8 MeV. This calibration procedure allows us to establish the  Figure 2: 60 Co calibration: black and red dots represent experimental data and Monte Carlo simulation respectively. The energy scale has been established with 3% accuracy. In addition to the full absorption peak, corresponding to events completely contained in the active volume of the detector, we can distinguish a bump of events corresponding to one of the two gammas escaping the detector. energy scale with a 3% uncertainty, dominated by systematics, ensuring that the detector gain remained constant, within uncertainty, during the whole period of data taking. Moreover, the main crucial parameters, i.e., PMT high voltage, Radon activity, temperature and humidity of the environment were also continuously monitored. The MC simulation of the detector response based on the Geant4 toolkit [7] has been developed including the position of nearby LVD counters (details can be found in [8]). Data and MC are shown in figure 2, with black and red dots, respectively.

Nuclear recoils
To calibrate and monitor during data taking the detector response to nuclear recoils, an AmBe(α, n) source was used, inserted inside the detector using a stainless steel rod and kept in the center of the scintillator volume as the 60 Co one.
This source emits about 10 neutrons per second, and the emitted neutrons are thermalized by elastic scattering on the hydrogen atoms of the organic scintillator. During the slowing down process, a fraction of the neutron kinetic energy is converted to photons which can be detected by the PMTs giving a prompt signal. After being thermalized, neutrons are eventually captured by Gd and H (and with lower probability by C and Fe nuclei) producing a delayed signal. Neutron captures The mean n-capture time can be measured by studying the distribution of the time intervals between the prompt signal, due to the neutron elastic scattering with H, possibly accompanied by 12 C * de-excitation gamma-ray of energy 4.44 MeV, and the n-capture signal. The mean n-capture time, strongly depends on the characteristics of the detector, mostly the Gd concentration (n-capture cross section), in our case 0.115% in weight, and the detector geometry. The simulation predicts a mean n-capture time t MC = 26.7 ± 0.2 µs, which is in agreement with the experimental value t exp = 26.9 ± 0.2 µs, and also in agreement with previous measurements and simulations [5]. This quantity allowed us to monitor the scintillator stability during data taking in terms of Gd concentration (in solution) and light transmittance.
Finally, from the energy spectrum of the prompt signals, and the known energy spectrum of the neutrons emitted by the AmBe source [9], we were able to measure the light output response to nuclear recoils of the detector. Details of this measurement can be found in [8] where the quenching factor was found to be in good agreement with Birks semi-empirical model with a Birks factor equal to 0.0140 ± 0.0007 g cm −2 MeV −1 .
In figure 4 is shown the energy spectrum of the detected prompt signals (after background subtraction) from the AmBe source (black dots), compared with the MC simulation (red dots). Note that the energy scale in abscissa is expressed in terms of MeV (electron equivalent), the same unit as in figure  2. The excellent agreement between the measured and simulated spectra validate the MC simulation and give us the confidence to use it for the reconstruction of the neutron flux from the cavern rock.

Experimental Method
The key idea of this work is identifying neutron events by triggering on the neutron capture on Gd, which has enough energy to emerge from the background, then, by looking backward in time, to possibly find the proton recoil signal induced by the same neutron. As a result, the neutron energy spectrum will be obtained by unfolding the proton recoil spectrum.
The three PMTs waveforms are continuously written in a circular memory buffer. When a trigger occurs, High Threshold trigger (HT from now on), a time window of ≈ 640µs around the trigger signal (2 16 samples for each PMT) is frozen. The DAQ software supports the Zero Length Encoding (ZLE) algorithm, providing data transfer in compressed mode. The ZLE threshold, which we will refer to as Low Threshold (LT) since it is lower than the trigger threshold, can be set independently on each input channel. The lower is LT the less is the energy of the neutrons we can detect. On the other hand HT has consequences on the detection efficiency but does not affect the energy range of the neutron spectrum.
HT has been set to have 100% efficiency at 4.0 MeV, and LT to have 100% efficiency at 1.0 MeV. In these conditions the trigger rate is approximately 0.8 Hz, therefore the detector runs with almost zero dead-time.
In figure 5 the energy spectrum of trigger signals is shown. The red histogram reports the energy distribution of all HT trigger pulses, the black histogram has been obtained rejecting "muon" candidates from the red one, i.e., events in coincidence with at least one of the 280 LVD counters which belongs to the same LVD tower [3]. The bump at about 8 MeV corresponds to the maximum of the gamma emission due to neutrons captured by Gd.
It can also be that the recoil proton has enough energy to trigger the data acquisition by itself, in this case the signal due to the n-capture will be following the trigger. To distinguish between the two signals of the pair the whole story of the event is needed. Therefore an acquisition window of 640µs is used with the trigger in the center, in order to have 320µs in both pre-trigger and post-trigger regions. Since the mean n-capture time is about 27µs, and considering a time delay of about 10µs for thermalizing, more than 98% of neutrons are expected to be captured in the 120µs closest to the trigger, that is the region we use for the proton recoil search. The remaining part of the waveform is used to measure the background.  In figure 6 the arrival time of pulses in the acquisition window is shown for four different trigger thresholds ranging from 3 to 6 MeV, from top to bottom in the figure. Triggers are not shown to focus the attention on the shape of the distributions, hence the missing points at 320µs. The distribution of the arrival time of the pulses is expected to have a flat component due to background signals (being uncorrelated to the trigger), while an exponential distribution is a distinguishing feature of a neutron event in which one pulse is due to the proton recoil and the other one to the n-capture. The exponential distribution emerges from the background when the trigger threshold is greater than 4 MeV. Its slope is in agreement with the expected one, then demonstrating that, at this energy threshold, the detector can discriminate between neutron events and gamma background.

Neutron Energy Spectrum
In this work, we use data from June 2011 to October 2012: over this period the detector was active for about 485 days, corresponding to 97% livetime. Neutron induced nuclear recoils are selected requiring (at least) 1 pulse of energy greater than 1.0 MeV (electron equivalent) occurring within ± 120 µs from the trigger signal, which has to be ≥ 5 MeV. Whenever the time difference exceeds this limit, the pair of pulses is considered as uncorrelated and therefore due to background.
Muon induced neutron events, which occur with a probability depending on the muon energy and the mass of the target nucleus, can be recognized because they are in time coincidence with the passage of the parent muon 1 . At the same time, in our detector, it is impossible to disentangle the two energy releases. In order to avoid overestimating of the high-energy component of the neutron spectrum, events tagged as muon, because in coincidence with a signal in at least one of the 279 surrounding LVD counters, were discarded from the analysis. Due to the low rate of the LVD detectors the dead time introduced is negligible. In this way we reject most of the neutrons produced by muons inside the experiment.
During 485 live days, after muon rejection and background subtraction, we found 360 nuclear recoil events. The energy distribution of these events is shown in figure 7 in black, compared with the background in magenta. The excess of events from 4 to 5 MeV is attributed to 4.44 MeV de-excitation γ-rays from the inelastic scattering 12 C(n, n ) whose contribution is anyway taken into account by the MC simulation and included in the response functions of the detector.
The method used to evaluate the neutron flux Φ n , from the nuclear recoil energy distribution, is based on the MC simulation of the detector response to mono-energetic neutrons. This simulation describes accurately both electron and nuclear recoil interactions in this same detector, see figure 2, 3, 4, and [8], when radioactive sources are placed in its center. Now the point source is replaced in the simulation by an isotropic neutron field, i.e., monochromatic neutrons are produced randomly on a spherical surface surrounding our detector together with the nearest LVD counters and including the GERDA water tank as well, which both contribute to reduce the absolute efficiency (see eq. 2) of the detector.
Let N MC gen be the total number of neutrons generated by the simulation, N MC imp the number of neutrons reaching the active volume, and N MC det the number of neutrons which satisfy the selection criteria (giving at least 1 MeV electron equivalent prompt signal and at least 5 MeV from the following capture) intrinsic and absolute efficiencies can be defined as: The intrinsic efficiency, int (E n ) drops to less than 1% for energies lower than 5 MeV, while it reaches ∼17% in the range from 15 to 20 MeV. Higher energies have not been studied in this work. Clearly, neutrons of energy equal to 5 MeV or less have a very low chance to give a prompt signal greater than 1 MeV electron equivalent, because of the scintillator light yield (quenching effect). This is the reason why, in the final result, the neutron spectrum starts at 5 MeV.
Considering the proton recoil spectrum of figure 7, the number of counts C i detected in the i-th energy bin is related to the incoming neutron spectrum Φ(E) by a Fredholm integral equation of the first kind: where S i (E) is the probability that a neutron of energy E deposits an energy between E and E + ∆E in the detector, corresponding to the i-th energy bin. This equation can be written in discrete form: where Φ n is the neutron flux in the n-th energy interval and S in is a squared matrix whose rows are the visible-energy distribution for mono-energetic neutrons, i.e., the response functions. The vector Φ can be (uniquely) determined by minimization only if its length is equal to the rank of the matrix S in . Since this is usually not the case for experimental data, which are affected by uncertainties, equation 4 can not be solved exactly. To avoid negative solutions and easily handle errors, the matrix inversion is performed using the least-squares approach, finding the vector Φ that minimizes S Φ − C .
The result of the deconvolution (or unfolding) process is the spectrum Φ shown, with purely statistical error bars, in figure  8 in black, compared with other measurements and simulations existing in literature, in the same figure.
To double-check that the matrix inversion algorithm worked correctly we folded again the neutron spectrum with the response functions finding good agreement with the experimental data C i .
At energies higher than 10 MeV the fraction of neutrons induced by muons that are not detected by other LVD counters and then untagged, becomes dominant. This component was not included by Wulandari, et al. [11], who simulated the flux of neutron induced by radioactivity only, that's why their spectra (red and blue empty circles in figure 8) end at about 8 MeV. Measurements of the neutron flux, integrated over a broad range of energies, have been published along the years and can be found in: [14], [15], [16], [17] and [18]. They are not shown in figure 8.

Conclusions
A 1.5 cubic meter detector located in the INFN Laboratori Nazionali del Gran Sasso (LNGS) filled with Gd-doped liquid scintillator has been used to measure the energy spectrum of neutron from rock radioactivity. The dimensions of the detector, allows the confinement of the entire event, and the strong signature of the neutron events, provided by n-capture on Gadolinium, allows a direct measurement of the energy spectrum of scattered protons. Previous studies of the quenching effect, made on this same detector [8], allow a correct interpretation of this energy spectrum. Deconvolution of kinematics and detection effects leads to a final reconstructed neutron spectrum obtained in the absence of any hypothesis on its shape.
The obtained results are limited to the neutron energy range between 5 and 15 MeV and take already into account the existence of the GERDA water tank in hall A, which reduces the detected neutron flux of about 30%. At energies lower than about 5 MeV, mainly because of quenching, neutrons have a very low chance to give a prompt signal greater than 1 MeV electron equivalent that is our energy detection threshold for recoils. For energies greater than about 10 MeV, on the other side, the contamination due to untagged muon induced neutrons becomes dominant. In this energy region our results can not be compared with MC simulations from [11], because this component was not considered by the authors being strongly dependent on how the experimental hall itself is filled.  [13] and orange from [12]. Empty circles, in red and blue, represent the results of the simulations from [11] of Hall A dry and wet concrete respectively. Each point in the figure shows the integral flux in a 1.0 MeV energy bin.