Understanding NaI(Tl) crystal background for dark matter searches

We have developed ultra-low-background NaI(Tl) crystals to reproduce the DAMA results with the ultimate goal of achieving purity levels that are comparable to or better than those of the DAMA/LIBRA crystals. Even though the achieved background level does not approach that of DAMA/LIBRA, it is crucial to have a quantitative understanding of the backgrounds. We have studied background simulations toward a deeper understanding of the backgrounds and developed background models for a 9.16-kg NaI(Tl) crystal used in the test arrangement. In this paper we describe the contributions of background sources quantitatively by performing Geant4 Monte Carlo simulations that are fitted to the measured data to quantify the unknown fractions of the background compositions. In the fitted results, the overall simulated background spectrum well describes the measured data with a 9.16-kg NaI(Tl) crystal and shows that the background sources are dominated by surface $^{210}$Pb and internal $^{40}$K in the 2 to 6-keV energy interval, which produce 2.4 counts/day/keV/kg (dru) and 0.5 dru, respectively.


Introduction
Numerous astronomical observations have led to the conclusion that the majority of the matter in our universe is invisible, exotic, and nonrelativistic dark matter [1,2]. However, it is still unknown what the dark matter is. Weakly interacting massive particles (WIMPs) are a Corresponding author: ejjeon@ibs.re.kr one of the most attractive dark matter particle candidates [3,4]. The lightest supersymmetric particle (LSP) hypothesized in theories beyond the standard model of particle physics is a suitable candidate for a dark matter WIMP. There have been numerous experiments that directly search for WIMPs in our galaxy by looking for nuclear recoils that are produced by WIMP-nucleus scattering [5,6].
To date, no other experiments, except for the DAMA experiment [7,8,9], have found an annual modulation signal interpreted as WIMP interactions with a significance of 9.2σ. However, this finding has spurred a continuing debate since the WIMP-nucleon cross sections inferred from the DAMA modulation are in conflict with limits from other experiments that directly measure the nuclear recoil signals, such as XENON100 [10], LUX [11], and SuperCDMS [12].
The Korea Invisible Mass Search (KIMS) is an experiment that aims at searching for dark matter at an underground laboratory located in Yangyang, South Korea (Y2L). We are performing a high-sensitivity search for WIMP interactions in an array of NaI(Tl) crystals in an attempt to reproduce the DAMA/LIBRA's observation of an annual modulation signal [24,25].
There are several groups, such as DM-Ice [13,14], ANAIS [16,18], and SABRE [19], developing ultra-lowbackground NaI(Tl) crystals with the goal of reproducing the DAMA/LIBRA results and currently, KIMS and DM-Ice have agreed to operate a single experiment, COSINE, at Y2L using NaI(Tl) crystals and a total mass of 106 kg is being used in the first-stage experiment, COSINE-100. As part of this program we have developed ultra-low-background NaI(Tl) crystals and studied their properties in a variety of test setups with the ultimate goal of achieving purity levels that are comparable to or better than those of the DAMA/LIBRA crystals. Even though current background levels achieved by the research and development are higher than those of DAMA/LIBRA, it is crucial to have a quantitative understanding of the backgrounds.
For further understanding of the backgrounds, we have performed Monte Carlo simulations based on Geant4 and compared their results with measured data (see Sect. 3.2). To build concrete background models for a 9.16-kg NaI(Tl) crystal used in one of test setups, we studied background simulations of internal radioactive contaminants, such as natural radioisotopes inside NaI(Tl), cosmogenic radionuclei, and surface contaminations in NaI(Tl) crystal (see Sect. 3.2.2 and 3.2.3), and external background sources from the exterior of crystals (see Sect. 3.2.1). We quantified their contributions by treating them as floating and/or constrained parameters in the data fitting (see Sect. 3.2.4). In addition, our evaluation of background prospects, based on this study, is described in section 4.

Detector shielding and configuration of the test arrangement
We have studied low-background NaI(Tl) crystals with various test setups at Y2L, as reported in references [24,25], and as shown in Fig. 1(a), the crystals being tested are enclosed by five different shield layers that were used for the KIMS-CsI experiments [21,22]. The outmost layer is a 30-cm-thick muon detector (MD) filled with mineral oil. The other layers include a sequence of 15-cm-thick lead, an iron sheet, 5-cm-thick polyethylene (PE), and a 10-cm-thick copper shield.
One test arrangement that was used for this study is shown in Fig. 1(b) with the CsI(Tl) crystal array. Three different-sized NaI(Tl) crystals, NaI-001, NaI-002, and NaI-005, are surrounded by ten CsI(Tl) crystals located inside the copper shield. Each end of the crystal was attached to a photomultiplier tube (PMT). NaI-001, NaI-002, and NaI-005 were produced by Alpa-Spectra (AS) and the first two crystals were grown from their own powders. For NaI-005, it was grown from AS WIMPScint-II grade (AS-WSII) powder that was aimed at reducing 210 Pb contamination in the powder and crystal growing process. As listed in Table 1 [25], the total α rate of NaI-005 was reduced by more than a factor of three compared to NaI-002.
In this paper, thus, we focused on the background model for the NaI-005 crystal and studied Geant4-based simulations of the background sources internal/external to NaI-005 in this test arrangement. Internal backgrounds of the two crystals, NaI-001 and NaI-002, were simulated for coincidence data of NaI-005 and backgrounds from the PMTs were also simulated for external backgrounds of NaI-005.
2.2 Background measurements in the NaI(Tl) crystal test setup NaI-002 and NaI-005 crystals have an identical cylindrical shape with a diameter of 4.2 inches, a length of 11 inches, and a mass of 9.16 kg. The light yield and the measured background rates from internal radioactive contaminants in the NaI(Tl) crystals are listed in Table 1. The details of the three crystals are discussed in reference [25].
Three different types of PMTs were used in the test arrangement: a metal-packed R11065SEL, a glasspacked R12669SEL, both manufactured by Hamamatsu Photonics, and 9269QA of Electron Tubes, Ltd.. R12669SEL PMTs were coupled with NaI-002 and NaI-005; R11065SEL and 9269QA were coupled with NaI-001 and CsI(Tl) crystals. The radioactivity levels of the PMTs were measured underground with a high-purity Ge (HPGe) detector and their measurements are listed in Table 2 [24].
We used the measured activities inside the crystals and from the PMTs for the simulation study of the NaI-005 backgrounds. In addition to natural radioisotopes measured, there are also backgrounds from cosmic excitation that are continuously decreasing as a function of time. We considered I, Te, and 22 Na isotopes from cosmic radiation as background sources and they are included in the simulations and their contribution described later.

Method of simulation
To understand the measured background of NaI-005 in the test arrangement described in section 2.1, we have performed Monte Carlo simulations with the Geant4 toolkit [27], version of 4.9.6.p02, including physics list classes of G4EmLivermorePhysics for low energy electromagnetic process and G4RadioactiveDecay for radioactive decay process.
In the simulation framework, radioactive sources such as full decay chains of 238 U and 232 Th were simulated, assuming that decay chains of 238 U and 232 Th are each in equilibrium, thus all related activities within the chains are simply equal to the 238 U and 232 Th activities multiplied by the branching ratios for decay of the   daughter isotopes. However, it needs to specify all longlived parts of 238 U and 232 Th daughters to consider broken decay chains when they are fitted to the measured data to quantify their unknown background fractions and, thus, we grouped daughter isotopes from the full decay chains of 238 U and 232 Th according to their halflives. We used five groups for 238 U and three groups for 232 Th, as listed in Table 3, in the data fitting. 40 K is treated as its own, additional group.
Each simulated event includes all energy deposited in the crystals within an event window of 10 µs from the time a decay is generated, to account for the conditions in the data acquisition system (DAQ) of the experimental setup [24]. Sometimes decays with relatively short half-life such as 212 Po decay (with a half-life of 300 ns) and the following decays will appear in the same event, called pileup events, and they were treated within one event in simulated event.
The simulated spectrum was smeared by the energy resolution, as a function of energy, that was measured during the calibration. Calibration points were measured using γ-ray sources, 241 Am and 60 Co; they contribute peaks at 59.5 keV( 241 Am), 1170 keV, and 1332 keV. Measured background spectra from several radioactivities were used for additional calibrations; there are peaks at 3.2 keV and 1460 keV from 40 K, 30.8 keV from 126 I, 67.2 keV from 125 I, and 609.3 keV and 1120.3 keV from 214 Bi, respectively. We also used the peak at 0.8 keV from 22 Na. To build the concrete background model for NaI-005 it needs to compare the measured background level in NaI-005 with the simulated background spectra that are each composed of different background compositions. Initially, we studied background simulations of natural radioisotopes in NaI(Tl) crystals and PMTs and compared the total of the simulations with the data. In this comparison, we found there remains significant backgrounds to be modeled and, thus, we considered additional background components to improve the model. The details of the simulation studies are described in subsections.
For the initial comparison, we simulated daughter isotopes from full decay chains of 238 U, 232 Th, and 40 K located inside NaI(Tl) crystals and 26 PMTs assuming a chain equilibrium. By using the measured activities of radioactive sources listed in Tables 1 and 2, event rates were normalized to the units of counts/day/keV/kg (dru). Fig. 2 shows the normalized background energy distributions in NaI-005 for single-hit events in which there is an energy deposit in NaI-005 only. The total of the simulations (solid red line) is compared with the measurement (open black circles). In the simulations, the background sources for energies below 10 keV are dominated by 40 K (solid blue line) and 210 Pb (solid cyan line) internal to the NaI(Tl) crystal, and for high energies above ∼100 keV external backgrounds from PMTs are dominant. The details are itemized in the following.
However, as shown in Fig. 2, the overall measured background level is higher than that of background simulations and, thus, there are significant remaining backgrounds. For the contributions of the background sources not included in Fig. 2, we studied background simulations of the external sources, cosmogenic radionuclei, and surface contaminations of the NaI(Tl) crystal. To quantify the contributions of all background sources in the simulation they are fitted to the measured data by floating and/or constrained the unknown fractions of the background spectra as fit parameters. In the fit, daughter isotopes from full decay chains of 238 U, 232 Th, and 40 K are grouped into 9 background spectra: 5 background spectra for 238 U, 3 background spectra for 232 Th, and a spectrum for 40 K. The fitted results for the contributions of those background sources to the total background level are described in sections 3.2.1, 3.2.2, and 3.2.3, quantitatively. All the fit results of activities of background sources and background events in the 2 to 6 keV energy region are listed in Tables 4 and 5 in section 3.2.4.
• Internal backgrounds of NaI(Tl) crystals To normalize the backgrounds from internal radioactive contaminants, we assumed a chain equilibrium. Therefore, all related activities within the chains are equal to 238 U, 232 Th, and 40 K activities, in Table 1, multiplied by the branching ratios for decay of the daughter isotopes. We also added the background simulation of internal 210 Pb by considering the measured α rate. The resultant background contributions, except for those from 40 K and 210 Pb, were negligible (< 10 −2 dru).
• External backgrounds from PMTs We used 26 3-inch PMTs of three different PMT types in the test arrangement. We thus considered measured activities in terms of different PMT types, as listed in Table 2. The normalized background contributions from the PMTs are represented in three radioactive sources : 238 U (solid magenta line), 232 Th (solid khaki line), and 40 K (dotted blue line).

Effects of external background sources
The majority of the external γ background comes from the radioactive decay of isotopes in the surrounding rocks. To block such an environmental background we installed several shielding layers, as described in section 2.1, and measured the background reduction with an ultra-low background 100% HPGe detector at Y2L. By using the full shielding structure and N 2 gas flowing into the inside of the copper shield to avoid backgrounds from 222 Rn in the air at Y2L, we could reduce the environmental background by a factor of 10,000 (measured to be 1.20 ± 0.49 pCi/L [23]), thus ensuring that those contributions would be negligible. However, there exist some background from radioactive sources in detector components inside the shielding. Typically, PMTs and materials-such as the copper shield, the copper structure to which the PMT is affixed to the crystal, the PMT base, connectors, and bunches of cables-would contain radioactive sources that can contribute to the background as an external background source inside the shielding. For PMTs, even though we measured the radioactivity levels, background contributions will be different for the main components in the PMTs, such as the PMT window, body, and stem, where radioactive sources were generated.
Although it is difficult to reproduce accurately background contributions from external sources without knowing well about the external background contamination, it is possible to estimate the effects of such radioactive sources, external to the NaI(Tl) crystal, by considering them as parameters in the fit. To take it into account we simulated the background spectra in such a way that radioisotopes contained in 238 U, 232 Th, and 40 K were generated randomly in the whole PMT body including the PMT window and stem, and we grouped the results into nine background spectra, instead of performing each of the simulations of all the external sources besides PMTs. It is because that the background energy distribution in NaI-005 is very similar to that from PMTs when we simulated external background sources from the inside space of the copper shield. We also separately considered background contributions from the two PMTs attached to NaI-005, apart from the other PMTs attached to NaI-001, NaI-002, and CsI(Tl) crystals, because the external background contributions to NaI-005 will be different for the distance from the external sources to the crystal. Consequently, we have two sets of nine background spectra for testing the effects of the external background sources. They were fitted to the measured data, as floating and/or constrained parameters, to estimate their fractions contributed to the total background and, in the fit, we used singlehit events and multiple-hit events, the latter of which have energy deposits in two or more crystals, simultaneously. We estimated the external background contributions by assuming that the fitted results for PMTs simulations include the backgrounds from all the external sources. As shown in Fig. 3, the fitted results well reproduce the measured data for both single-hit events and multiple-hit events for energies above ∼100 keV. The black circles represent the data and the solid red line represents the sum of simulated background spectra with the fitted fractions of groups 1 (dotted magenta line), 4 (solid magenta line), 7 (solid khaki line), and 8 (dotted khaki line) and 40 K (dotted blue line); the biggest contribution is from 226 Ra (group 4). Those external backgrounds are expected to be vetoed by an active veto detector in the new detector design.

Backgrounds from cosmogenic isotopes
We installed NaI-005 in the test setup in December 2014, soon after that crystal was delivered to Y2L, and we recorded data for about a month. Therefore, some backgrounds are expected from cosmogenic activations. We checked backgrounds over the specified time interval and observed a clear reduction of the peak at 30.8 keV from 126 I within the first 10 days. We thus did not use  the first 16 days data to exclude the background from 126 I in this study.We simulated backgrounds from cosmogenic isotopes, which are expected to be produced by cosmic ray exposure [15]. The simulated background spectra are used in the data fitting, by floating their unknown fractions, and the fitted results are shown in Fig. 4. The dominant isotopes are itemized in detail in the following. • 121 Te (solid gray line) The half-lives of 121 Te and 126 I are 19.17 and 12.93 days, respectively, and a clear reduction of the peak from 126 I within the first 10 days was observed. However, there is still peak at 573 keV identified with γ rays emitted by electron capture in the decay of 121 Te that is produced as a result of the decay of the longer-lived cosmogenic isotope 121m Te (halflife = 164.2 days), which is reported in reference to ANAIS [17].
• 125m Te (dotted brown line) and 123m Te (solid brown line) 125m Te and 123m Te are long-lived metastable states and their half-lives are 57.4 and 119.2 days, respectively. They contribute peaks at 145 and 248 keV.
• 22 Na (solid cyan line) 22 Na can be produced through the (n, 2n) reaction on 23 Na by energetic cosmic neutrons at sea level. It decays via positron emission (90%) and electron capture (10%), followed by 1270-keV γ-ray emission with a mean lifetime of 3.8 yr. The electron capture decay produces ∼0.8-keV X-rays. Therefore, ∼10% of the 22 Na decay will produce 0.8-keV X-rays and 1270-keV γ rays simultaneously. Meanwhile, the positron will be converted to two 511-keV annihilation γ rays. If one of the two 511 and 1270 keV γ rays escapes from a crystal, the energy deposited in the crystal will be 650-1000 keV. We looked for a coincident event that deposits 1270 keV of energy in NaI-002, resulting in a γ-ray hit in the 650 to 1000-keV energy interval in NaI-005. It is measured to be 0.8 ± 0.3 mBq/kg by considering the detection efficiency obtained from a Monte Carlo simulation and is used as a constrained parameter for the background fraction of 22 Na in the fit.

Background from NaI(Tl) surface
Energy (   The levels of 238 U and 232 Th contamination measured in NaI-005 are too low to account for the total observed α-particle rate, which suggests that they are due to decays of 210 Po nuclei that originate from 222 Rn contamination that occurred sometime during the powder and/or crystal processing stages. This is confirmed by the observation of a 46-keV γ peak that is characteristic of 210 Pb. There is also the possibility of surface contamination, which is expected to affect the lowenergy spectrum in a different way from bulk 210 Pb. To clarify the effect of surface 210 Pb contamination we simulated the background spectra at various surface depths from 0.1 to 10 µm by generating 210 Pb randomly within the depth. Fig. 5 shows a comparison of background spectra simulated in different surface depths. The low-energy background looks different from that due to bulk 210 Pb (cyan color), and the heights of the ∼50-keV peak are also different from that due to bulk 210 Pb.
In the simulation, we added a background contribution from surface 210 Pb generated randomly within a surface depth of 10 µm. Background fractions of bulk and surface 210 Pb are obtained in such a way that they are treated as floating parameters in the fit to the measured data. And they are found to be 0.05±0.76 and 0.81±0.25 mBq/kg, respectively, as shown in Fig. 6. The total amount of 0.86±0.80 mBq/kg is comparable to the total observed α activity from decays of 210 Po from crystals grown by AS-WSII powder (0.85 ± 0.04 mBq/kg).

Comparison of fitted Monte Carlo fractions with data
We simulated full decay chains of 238 U, 232 Th, and 40 K from NaI(Tl) crystals and 26 PMTs that are grouped into three sets of nine background spectra and additionally we simulated five cosmogenic isotopes and bulk and surface 210 Pb, including their progenies. Using all of the simulated background spectra, we fit the model to the measured data for both single-hit events and multiple-hit events of NaI-005, in the 2 to 1510-keV fitting range, to determine the unknown background fractions; we used a maximum likelihood fit using Poisson statistics [28]. In the fitting, nine of the internal background groups were constrained by the measured activity errors and the other groups for external backgrounds were treated as floating and/or constrained parameters. All the cosmogenic background spectra and bulk/surface Pb210 spectra were considered as floating parameters in the fit. The fitted results of activities of background sources in NaI-005 and PMTs are listed in Table 4. Fig. 7 shows the fitted results for all the simulated background spectra plotted as various lines with differ-   ent colors and styles for both single-hit events (a) and multiple-hit events (b) in the 0.5 to 2000-keV energy region. The overall energy spectrum summed over all simulations (solid red line) is well matched to the data (open black circles) not only for single-hit events but also for multiple-hit events. However, there are shown discrepancies between measurements and simulations for high energies above ∼1600 keV. It is because that the pulse shape above ∼1600 keV from the R12669SEL PMT is deformed and a significant fraction of the photoelectrons do not reach to the anode. Accordingly, the integrated charge signals of energetic events do not scale linearly with energy as one can see in Fig. 7. It is also the reason why we chose the fitting range up to 1510-keV.
The measured energy spectrum is compared with all the simulated spectra in the low-energy region below 20 keV in Fig. 8. As listed in Table 5, the main backgrounds in the 2 to 6-keV energy region are from surface 210 Pb and 40 K contaminations inside the crystal, the contributions of which are found to be 2.4 and 0.5 dru, respectively. Remnants of cosmogenic activation of Te (0.95 dru) still persist but will quickly be reduced with a lifetime of < 100 days. We also determined the external backgrounds, which amounted to 0.59 dru. In the low energy region below 0.5 keV the data is suppressed by the requirement for low energy noise rejection and there is shown less number of events in the data.

Conclusion
We have studied NaI(Tl) crystal backgrounds based on Monte Carlo simulation using the Geant4 toolkit. The studies show that all the simulated background spectra, normalized by measured activities and fitted fractions, describe the data well not only for single-hit events but also for multiple-hit events. According to the comparison between the Monte Carlo simulation and the data in the low-energy spectra, the background is found to be dominated by 210 Pb, mainly surface 210 Pb, which is due to 222 Rn exposure during crystal growing and/or handling procedures, and 40 K within the NaI(Tl) crystal, the background level of which is consistent with  ANAIS's expectation evaluated by assuming the activities of their characterized NaI(Tl) crystals [18]. External background contributions are expected to be tagged by an active veto detector with a liquid scintillator (LS) surrounding crystals in the new detector design. We already observed such a reduction, 0.76 ± 0.04 in the 6 to 20-keV energy region, with a prototype active veto system using an LS [26]. An LS veto detector can also reduce the contribution of the internal 40 K background, by tagging 3-keV X-ray events with required conditions for the LS veto signal, by a factor of 2 with the optimized thickness of the LS in the new detector design [26]. In addition, we are studying the suppression of the 210 Pb crystal-surface background to achieve a background level as low as bulk 210 Pb, which contributes 0.04 dru in the 2 to 6-keV energy region. Moreover, improving the purity of NaI(Tl) crystals with small concentrations of 210 Pb, ∼60 µBq/kg, is possible, as reported in reference to KamLAND-PICO [20].
As a result, we are expecting that we can reach a background level of < 0.5 dru in the 2 to 6-keV energy region.