A measurement of the scintillation decay time constant of nuclear recoils in liquid xenon with the XMASS-I detector

We report an in-situ measurement of the nuclear recoil (NR) scintillation decay time constant in liquid xenon (LXe) using the XMASS-I detector at the Kamioka underground laboratory in Japan. XMASS-I is a large single-phase LXe scintillation detector whose purpose is the direct detection of dark matter via NR which can be induced by collisions between Weakly Interacting Massive Particles (WIMPs) and a xenon nucleus. The inner detector volume contains 832 kg of LXe. $^{252}$Cf was used as an external neutron source for irradiating the detector. The scintillation decay time constant of the resulting neutron induced NR was evaluated by comparing the observed photon detection times with Monte Carlo simulations. Fits to the decay time prefer two decay time components, one for each of the Xe$_{2}^{*}$ singlet and triplet states, with $\tau_{S}$ = 4.3$\pm$0.6 ns taken from prior research, $\tau_{T}$ was measured to be 26.9$^{+0.7}_{-1.1}$ ns with a singlet state fraction F$_{S}$ of 0.252$^{+0.027}_{-0.019}$.We also evaluated the performance of pulse shape discrimination between NR and electron recoil (ER) with the aim of reducing the electromagnetic background in WIMP searches. For a 50\% NR acceptance, the ER acceptance was 13.7${\pm}$1.0\% and 4.1${\pm}$0.7\% in the energy ranges of 5--10 keV$_{\rm ee}$ and 10--15 keV$_{\rm ee}$, respectively.

A : We report an in-situ measurement of the nuclear recoil (NR) scintillation decay time constant in liquid xenon (LXe) using the XMASS-I detector at the Kamioka underground laboratory in Japan. XMASS-I is a large single-phase LXe scintillation detector whose purpose is the direct detection of dark matter via the NR resulting from collisions with Weakly Interacting Massive Particles (WIMPs). The inner detector volume contains 832 kg of LXe. 252 Cf was used as a source of irradiating neutrons outside of the detector. The decay time constant of the resulting neutron induced NR was evaluated by comparing the observed photon detection times with Monte Carlo simulations for various decay time constants. Fits to the decay time prefer two decay time components, the singlet and the triplet state τ S and τ T respectively; with τ S = 4.3±0.6 ns taken from prior research, τ T was measured to be 26.9 +0.8 −1.2 ns with a singlet state fraction F S of 0.252 +0.028 −0.020 . We also evaluated the performance of pulse shape discrimination between NR and electron recoil (ER) with the aim of reducing the electromagnetic background in WIMP searches. For a 50% NR acceptance the ER acceptance was 13.7±1.0% and 4.1±0.7% in the energy ranges of 5-10 keV ee and 10-15 keV ee , respectively.

Introduction
Liquid xenon (LXe) has been used in many modern experiments such as dark matter and neutrinoless double beta decay searches [1][2][3][4][5]. The scintillation timing information can be used for position reconstruction of an event in the detector [6] as well as for particle identification and discrimination between nuclear recoil (NR), gamma-rays, and beta-rays [7].

The other process involves the recombination process between electrons and ions
Xe + + Xe → Xe + 2 , Xe + 2 + e − → Xe * * + Xe , Xe * * → Xe * + heat , The Xe * 2 dimer has a singlet and triplet state, each with its own decay time constant. The decay time constants of the singlet and triplet states were reported to be 4.3±0.5 ns and 21±2 ns, respectively using 252 Cf fission fragments [10], and the recombination process with its longer decay time constant of more than 30 ns was measured with 1 MeV electrons from a 207 Bi source [10,11]. In the case of neutron induced NR events, the decay time constant of the triplet state was reported to be ∼20 ns for both under the zero-electric field [13,17] and the non-zero electric field (0.1-0.5 kV/cm) [18,19]. While the decay time constants of the singlet and triplet states depend weakly on the density of the excited species, the ratio of singlet to triplet state as well as the recombination time depends on the deposited energy density [20]. Such pulse shape information may allow for the discrimination between NR and ER events [21].
XMASS-I is a large single phase LXe detector built primarily for dark matter searches, and previously reported the decay time constant of ER events using low energy gamma-rays calibration sources [8]. In this work, we measured the decay time constant of NR events using an external 252 Cf neutron source irradiating the XMASS-I detector and also evaluated the usefulness of pulse shape discrimination (PSD) between NR and ER in the energy region of interest for dark matter searches.

The XMASS-I detector
The XMASS-I detector is located in the Kamioka mine under 1,000 m of rock (2,700 meter water equivalent). The inner detector (ID) contains 832 kg of LXe inside a spherical, oxygen free high conductivity (OFHC) copper structure with an 80 cm diameter. Scintillation light from the LXe is detected by 630 hexagonal R10789 photomultiplier tubes (PMTs) and 12 cylindrical R10789Mod PMTs with a total photocathode coverage of 62.4%. The inner containment vessel contains the LXe and the PMT holder, while the outer containment vessel holds vacuum for thermal insulation. In order to reduce external gamma-rays and neutrons from the surrounding rock, the ID is placed at the center of the outer detector (OD). The OD is a cylindrical tank 10 m in diameter and 11 m in height filled with ultrapure water. 72 Hamamatsu 20-inch R3600 PMTs are mounted on the inner surface of the water tank to provide an active muon veto. More details can be found in Ref. [1].
Signals from the 642 ID PMTs were recorded by CAEN V1751 waveform digitizers with a 1 GHz sampling rate and 10-bit resolution. Analog-timing-modules (ATMs) that were previously used in the Super-Kamiokande experiment [22,23] worked for generating a trigger. The threshold for an ID PMT to register a hit in the ATMs is set at 0.2 photoelectron (PE). When 4 or more ID PMTs hit in a 200 ns coincidence timing window, a global trigger is issued to both the ATMs and the waveform digitizers. For each triggered event, the waveform is recorded with a width of 10 µs. As for OD PMTs, if 8 or more PMTs hit, an OD trigger is issued.

LED calibration
The individual PMT gains are monitored by a blue LED embedded in the inner surface of the PMT holder. This LED is flashed every second using the one-pulse-per-second signal from the global positioning system. LED calibration data is taken continually during the physics runs and identified by the trigger information.

Energy calibration and light yield
To check the stability of the detector's light yield, inner calibration data using a 57 Co source are taken on a weekly or every two weeks basis. Deploying the 57 Co 122 keV gamma-ray source at the center of the detector, we obtain the quoted photoelectron yield as ∼15 PEs / keV and also trace the timing offsets of the PMT channels. In addition, the 55 Fe, 109 Cd, 241 Am, and 137 Cs sources are also used for the position-dependent energy calibration of the detector.

The neutron source and its deployment
252 Cf undergoes spontaneous fission with a branching ratio of 3.11%. An average fission event emits 8 gamma-rays with a total energy of 7 MeV and 3.75 neutrons [24].
This gamma-ray emission is used to tag such fission events. Figure 1 shows the calibration setup. To detect the gamma-rays, the 252 Cf source was deployed in the cylindrical Bicron BC400 plastic scintillator 40 mm in diameter and 50 mm long. It has a central hole 9 mm in diameter and 12 mm deep to place the source. The scintillator was read out by a Hamamatsu R580 1.5 inch PMT. Hereafter this scintillator and PMT were referred to as neutron tagging assembly. A timing calibration between the plastic scintillator and XMASS-I detector electronics was performed with a 60 Co source instead of the 252 Cf source. The signal of the plastic scintillator was recorded by the same waveform digitizer as ID signals. A cylindrical polyethylene pipe with a 45 mm outer diameter, 10 mm inner diameter, and a 300 mm length was installed in front of the plastic scintillator and worked as a support structure. A 10 mm long lead plug in front of the source prevents 252 Cf gamma-rays from reaching the detector. For deployment at the detector, this assembly was inserted into a 45 mm inner diameter stainless steel calibration pipe that reached through the water tank toward the outer containment vessel of the inner detector. A stainless steel rod allowed to insert and position this neutron tagging assembly in the calibration pipe. The position of the neutron tagging assembly was adjusted at 150 mm from the tip of the stainless steel calibration pipe to keep the trigger rate less than 100 Hz from the requirement of the capability of the data acquisition.

Monte Carlo simulations
The XMASS Monte Carlo (MC) simulation is based on Geant4 [25]. It includes a detailed detector geometry, particle tracking, the scintillation process, photon tracking, the PMT response, and electronics responses. As for the input parameters to the simulation, optical parameters such as the absorption and scattering lengths of the LXe are extracted from the comparison between data and simulation for the position-dependent inner calibration with the 57 Co source. The copper reflection was deduced from the comparison of gamma-rays from the surface of the detector between data and MC simulations. Table 1 summarizes the input parameters of the MC simulation.
In the simulation of the 252 Cf calibration, the Brunson model [24] and the Watt spectral model [28] were used for the input energy spectrum of the gamma-rays and the neutrons, respectively. For the cross section of xenon-neutron elastic scattering, inelastic scattering and neutron capture, we used both the ENDF/B-VII.0 library following the instructions in Ref. [29,30] and the G4NDL3.13 library based on the ENDF/B-VI library, and compared the results to evaluate the systematic uncertainty of the cross section. For NR events, we considered the relative scintillation efficiency, named L eff [31], which is defined as the scintillation yield of xenon for NR relative to the zero-field scintillation yield for 122 keV gamma-rays from 57 Co. As for ER events, the non-linearity of scintillation light yield was taken into account using a model from Ref. [32] tuned with our gammaray calibration data. Since the rate for multiple neutrons or gamma-rays from the same fission event to enter the ID simultaneously was found to be negligible in this measurement, we generated a single neutron or a gamma-ray for each MC event in the 252 Cf simulation while considering their intensity. The detection time of the i th photon T i after a 252 Cf spontaneous fission (T = 0) was defined as Here t Edep is the time when the incident particle deposits its energy in the LXe. Scintillation time t i scinti follows the scintillation decay time profile parameterized as We assumed that the scintillation decay time constant has two components τ S and τ T , corresponding to the decay constants of the singlet and the triplet state with the effect of recombination, respectively. The recombination process contributes at most 10% of total scintillation light for nuclear recoil [33]. F S denotes the fraction of photons following τ S . t i TOF is the time of flight (TOF) of the scintillation photon. Here, the group velocity of the scintillation light was calculated from the refractive index of LXe. t i TT is the transit time in the PMT, which we assume to be the same for all PMTs. The transit time spread (TTS) of σ = 2.4 ns for our PMT [34] was included in the timing calculation. t i jitter is a smearing parameter accounting for the jitter in the electronic channel of PMT and extracted from the 57 Co calibration data. It follows a Gaussian distribution with a standard deviation of 0.93 ns. After calculating T i for all photons, waveform for each PMT was simulated using the one PE pulse shape extracted from LED calibration data.

Event selection
We took 1.5 hours of 252 Cf source data with an ID trigger rate of roughly 80 Hz. The signal of the plastic scintillator was searched for by offline analysis. The energy threshold of the neutron tagging assembly was found to be about 100 keV, based on the comparison between 60 Co data and simulation. Events related to neutrons from the 252 Cf source were selected using the following four criteria. (1) The event has a signal from the R580 1.5 inch PMT within +2 −1 µs from the timing of the

Figure 2.
Timing distributions of 252 Cf source data (solid black) and simulation (solid red). ∆T is defined in section 3.2. Timing of each event is aligned so that the events induced by the 252 Cf gamma-rays peak at ∆T = 0 ns. Magenta dotted and green dash-dotted histograms show the timing distribution from the 252 Cf neutron simulation (single site events only) and the 252 Cf gamma-ray simulation, respectively. The accidental coincidence contribution was derived from event rate in -1900 < ∆T < -100 ns time window. This accidental rate was about 3 × 10 −4 events/s and was subtracted from data. Events in 30 < ∆T < 100 ns were used for NR events analysis. ID trigger. (2) Only the ID trigger is issued. (3) The time difference from the previous ID event is longer than 500 µs and the root mean square of the timings of all hits in the event is less than 100 ns. This removes noise events that often follow particularly high energy events. (4) The time difference (∆T) between the time when the neutron tagging assembly detects a signal and the time when the ID trigger issued is 30 < ∆T < 100 ns. As the response of the neutron tagging assembly is not implemented in our MC simulation, ∆T is calculated from the time difference between the time of the particle generation and the time that the ID trigger issued. Figure 2 shows the ∆T distributions of events which have less than 500 detected PE in the ID. In the figure, the definition of single-site (multi-site) in the simulation is that all photons detected by PMTs originate from single (multiple) energy deposition in the LXe. Figure 3 shows the energy spectrum after all cuts. It includes the systematic uncertainties related to the detection efficiency of the plastic scintillator (±10%), the cross section difference between the G4NDL3.13 library and the ENDF-B/VII.0 library (25% at most), the scintillation efficiency for NR (±1σ in Ref. [31], 10% at most), and the position dependence of energy scale (5%), which accounts for the discrepancy of the detected PE distribution between data and simulation. The detection efficiency of the plastic scintillator was estimated to be (70±10)% by comparing the count rate for the 137 Cs and 60 Co source in data and MC simulation using a small setup. The ∆T distribution of the simulation agreed with data as shown in Fig. 2. There was about 25% count rate difference below 20 detected PEs in the energy spectra in Fig. 3, we discuss the impact on scintillation decay time constant in the section 4.1.
The contribution of the multi-site events is about 75% as deduced from MC simulation. Events with PE counts between 10 and 100 PEs, corresponding to the energy range from 1.5 keV ee (6.3 keV nr ) to 8.3 keV ee (40 keV nr ), were used to evaluate the NR decay constant described in the following section.

Evaluation of the nuclear recoil decay time constant
The decay time constant of NR events was evaluated by comparing the raw timings of the detected PEs over all PMTs between data and MC simulation with various timing parameters. To analyze waveforms of individual PMTs, we developed a peak finding algorithm based on a Savitzky-Golay filter [35] to obtain individual photon hit timings. Each peak was fitted with a single PE waveform template obtained from the LED calibration data. We obtained the timing and PE information for each PE detected at each PMT. Figure 4 shows a NR event waveform from a single PMT with the fitting result. Due to fluctuations of the baseline and electronic noise, the peak-finding algorithm sometimes misidentifies the tail of the single PE distribution as a peak. Such misidentified peaks typically have PE smaller than 0.5 PE. In this study, only peaks that have more than 0.5 PE are used. For each event, all peaks from all PMTs are sorted in order of detected timings, and the timing of the fourth earliest peak is set to T = 0 ns, reflecting the trigger implementation in DAQ. All peak timings in an event are shifted relative to the time of the fourth earliest peak, to sum up all events.
To obtain the decay time constant for NR, we performed a χ 2 fit defined as

The scintillation decay time constant of the nuclear recoil
The simulation using (τ S ,τ T ,F S ) = (4.3 ns, 27.0 ns, 0.25) gives the smallest χ 2 / ndf = 114.6 / 115. Data overlaid with the best fit simulation is shown in Fig. 5. Figure 6 shows the χ 2 map in the F S -τ T plane. Since the parameters were scanned in discrete steps, a parabolic function fit as defined in Eq. 4.1 was performed to obtain the decay time constant.  were discovered to be (4.3 ns, 26.9±0.5 ns (stat), 0.252±0.013 (stat)). All systematic uncertainties are listed in Table 2. They were evaluated as follows: Table 2. Summary of systematic uncertainties on τ T and F S .
(2)Neutron cross section: The event rate, including the fraction of multi-site events, depends on the neutron cross-section. We found an event rate difference of about 10% between the simulation using the ENDF-B/VII.0 and the simulation using the G4NDL3.13 library. Parameter scans of F S and τ T with the G4NDL3.13 library were conducted and the difference of the respective best fit value was used as the systematic uncertainty due to the neutron cross section. The effect of the count rate difference mentioned in the section 3.2 was also evaluated. We evaluated the impact on the scintillation decay time constant by lowering the weight of events in the data below 20 PEs and it turned out to be a negligible effect.
(3)L eff : Following Ref. [31]'s error estimates, we ran simulations also with L eff ±1 σ. Reanalyzing using these simulations we used the differences of the respective best fit value as the systematic uncertainty due to L eff .
(4)Jitter: Timing jitter affects the determination of the rising edge of timing distributions. The uncertainty was evaluated by comparing the timing distribution of data and simulated samples with different assumptions for the amount of jitter. t jitter was changed from σ = 0.93 ns to 0.0, 0.5, 1.5 ns and the differences of the respective best fit value were assigned as systematic uncertainty.
(5)Afterand pre-pulses: While studying our PMTs we found single-PE level after-pulses about 40 ns after the main pulse and pre-pulses about 15 ns before. The rate and timing of both the after-and pre-pulses were measured independently in a laboratory setup. These afterand pre-pulses are not included in our standard MC simulation. To study the impact of these missing pulses, we added the appropriate rates of after-and pre-pulses with Gaussian distributions in the timing offset: 40 ns offset, 5 ns width and 0.65% for the after-pulses, and -15 ns offset, 2 ns width and 0.10% for the pre-pulses for each peak in our standard 252 Cf simulation. The difference between the two evaluation on the standard and the after-and pre-pulse enhanced simulation is used as systematic uncertainty.  Figure 7. The best-fit parameters τ T and F S with various measurements. Filled markers and solid lines correspond to NR measurements. Open markers and dotted lines correspond to ER measurements. Results from Akimov et al. [13], Teymourian et al. [15], the LUX experiment [18], Hogenbirk et al. (0 V/cm, 0.1 kV/cm and 0.5 kV/cm) [19], and XMASS measurement (This work and [8]) are indicated by the black triangle, cyan cross, yellow diamond, green cross, blue circle, magenta triangle, and red square, respectively.
The obtained τ T is close to that for ER (27.8 +1.5 −1.1 ns using 55 Fe 5.9 keV gamma-rays ) reported in Ref. [8], although the obtained F S is larger than that of ER (0.145 +0.022 −0.020 using 55 Fe 5.9 keV gammaray). Figure 7 shows the τ T and F S with various measurements including both ER measurements and NR measurements. This measurement had the lowest energy threshold in all experiments conducted without an external electric field. D. Akimov et al. reported decay constants using a single component exponential fit [13]. Our single component fit value of τ = 22.5 ns for 1.5 < E < 8.3 keV ee is close to their reported value, although the simulation does not reproduce the data well ( χ 2 / ndf = 368.1 / 116 ). Our evaluated singlet fraction is consistent with the results of Ref. [18,19]. On the other hand, our τ T is about 5 ns longer than the values measured under the electric field. This might be an effect of the recombination process.

Performance of the pulse shape discrimination
For Weakly Interacting Massive Particle (WIMP) searches in data from a single phase LXe detector, the possibility of PSD between NR and ER is of significant interest. We evaluated the performance of PSD in XMASS-I based on our decay time constant measurement. To obtain the relevant timing distributions, we first simulated the ER events with uniform energy from 0 to 20 keV and NR events which followed the energy distribution of 100 GeV/c 2 WIMPs elastically scattering in the LXe target at the geometrical center of the detector. Figure 8 shows the peak timing distributions of those simulated ER and NR events that had energy deposits from 5 to 10 keV ee . In Fig. 8, TOF was subtracted using the velocity of 110 mm/ns and the timing of the fourth earliest peak in each event was set to T = 0 ns again to reflect the trigger implementation in DAQ. As mentioned in section 4.1, F S in NR is larger than in ER. Therefore a difference in the timing distributions can clearly be seen. These histograms in Fig. 8 were used as the probability density function ( f ER,NR (t i )) of the PMTs hit timings and we evaluated the following log likelihood ratio This log likelihood ratio was calculated using PMT hit timing with T > 0 ns, after TOF subtraction and T 0 determination. TOF subtraction uses the reconstructed event position. The performance of this PSD method was evaluated using another MC simulation. Electron events with uniform energy from 0 to 20 keV and the NR events which followed the energy distribution of 100 GeV/c 2 WIMPs elastic scattering were now generated uniformly throughout the detector. Figure 9 shows the log likelihood ratio distribution of the simulated ER and NR events with energy deposition from 5 to 10 keV ee . The ER acceptances when requiring a 50% NR acceptance for energies between 5 to 10 keV ee and between 10 to 15 keV ee were estimated to be 13.7±1.0% and 4.1±0.7%, respectively. This corresponds to a S/ √ N ratio of 1.4 and 2.5, respectively. Events which were reconstructed within 20 cm from the detector center were used for the evaluation of the performance of PSD. In Fig.  10 the blue graph shows the ER acceptance as a function of energy. The performance of this PSD method when evaluated using the XMASS-I detector simulation is consistent with the performance that we reported previously using a small chamber [21]. We also evaluated the performance of this PSD in an ideal case where the measurement is not affected by the jitter or TTS (red graph in Fig.  10). In this ideal case, the PSD performance improves by about a factor of 2 between 5 and 7 keV ee , and by about one order of magnitude between 15 and 20 keV ee .

Conclusions
We evaluated the time profile of NR scintillation emission in LXe with the XMASS-I detector using the 252 Cf sources. Two decay components are needed to reproduce the timing distribution of the  Figure 8. Pulse timing distribution of the simulated ER and NR events with energy deposition from 5 to 10 keV ee . In this figure, TOF was subtracted and the timing of the fourth earliest peak in each event was shifted to T=0 ns to reflect the trigger implementation. Areas are normalized to 1. NR data. We obtained the decay time constant of triplet state τ T = 26.9±0.5 (stat) +0.6 −1.1 (sys) ns and the singlet fraction F S = 0.252±0.013 (stat) +0.025 −0.015 (sys) with a decay time constant of singlet state τ S = 4.3±0.6 ns taken from a prior research. This measurement had the lowest energy threshold without applying an electric field. We also developed a PSD method based on a log likelihood ratio. The ER acceptances with a 50% NR acceptance at energies between 5 and 10 keV ee and between 10 and 15 keV ee were estimated to be 13.7±1.0% and 4.1±0.7%, respectively.  Figure 10. ER acceptance as a function of energy. NR equivalent energy is indicated on the top scales. WIMPs with a mass of 100 GeV/c 2 were used in the NR simulation. Electrons with energy distributed uniformly between 0 and 20 keV were used in the ER simulation. Both simulations generated events uniformly distributed in the XMASS-I detector and events which were reconstructed within 20 cm from the detector center were used for the evaluation. The two graphs correspond to our log likelihood performance evaluated for a 50% NR acceptance. Error bars show the quadratic sum of systematic uncertainty and statistical uncertainty. The blue graph shows the result of our MC evaluation, and the red one represents an ideal case where measurements are not affected by time jitter in the electronics and TTS in the PMTs.
Tokyo, and partially by the National Research Foundation of Korea Grant funded by the Korean Government (NRF-2011-220-C00006).