Improved search for two-neutrino double electron capture on $^{124}$Xe and $^{126}$Xe using particle identification in XMASS-I

We conducted an improved search for the simultaneous capture of two $K$-shell electrons on the $^{124}$Xe and $^{126}$Xe nuclei with emission of two neutrinos using 800.0 days of data from the XMASS-I detector. A novel method to discriminate $\gamma$-ray/$X$-ray or double electron capture signals from $\beta$-ray background using scintillation time profiles was developed for this search. No significant signal was found when fitting the observed energy spectra with the expected signal and background. Therefore, we set the most stringent lower limits on the half-lives at $2.1 \times 10^{22}$ and $1.9 \times 10^{22}$ years for $^{124}$Xe and $^{126}$Xe, respectively, with 90% confidence level. These limits improve upon previously reported values by a factor of 4.5.


Introduction
Double electron capture (ECEC) is a rare nuclear decay process where a nucleus captures two orbital electrons simultaneously. There might be two modes of the process: where Z and A are the atomic number and atomic mass number of the nucleus, respectively. Detecting the neutrinoless mode of this process (0νECEC) would provide evidence for lepton number violation and the Majorana nature of the neutrino if observed. To release the decay energy in 0νECEC, there are two proposed mechanisms: the radiative and the resonant mechanisms. In the case of the radiative mechanism, the decay energy is carried away by emitting, for example, an internal Bremsstrahlung photon [1,2]. This process is, however, expected to have a much longer life-time than neutrinoless double beta decay. On the other hand, an enhancement of the capture rate by a factor as large as 10 6 is possible if the initial and final (excited) masses of the nucleus are degenerate [3][4][5][6][7][8]. Therefore, experimental searches for 0νECEC have been recently performed for a variety of candidate nuclei [9][10][11][12][13][14][15][16][17][18][19].
Although two-neutrino double electron capture (2νECEC) is allowed within the Standard Model of particle physics, only a few positive experimental results for 2νECEC have been reported: geochemical measurements of 130 Ba [22,23] and a direct measurement of 78 Kr [19,24] with half-lives of the order of 10 21 -10 22 years. Despite the nuclear matrix element for the two-neutrino mode differs from that for the neutrinoless mode, they are related to each other through the relevant parameters in a chosen nuclear model [25]. For instance, the nucleus' axial current coupling constant g A and the strength of the particle-particle interaction g pp in the quasiparticle random-phase approximation (QRPA) model are obtained from single βdecay and two-neutrino double beta decay measurements [26]. Measurements of the 2νECEC half-lives with various nuclei would shed new light on constraining these parameters.
Natural xenon contains 124 Xe (abundance 0.095%) and 126 Xe (0.089%), in which ECEC can be observed. 124 Xe has the highest Q-value among all the known candidate nuclei for ECEC at 2864 keV [27]. This Q-value is sufficiently large to open the β + EC and β + β + channels. The predictions in the literature for the half-lives of 124 Xe 2νECEC are spread over a wide range between 10 20 and 10 24 years [21,[28][29][30][31][32] depending on the models used for calculating the corresponding nuclear matrix element and the effective value of the nucleus' g A . Although 126 Xe can also undergo 2νECEC, the life-time of this process for 126 Xe is expected to be much longer than that for 124 Xe since its Q-value is smaller at 920 keV [27].
Previous experimental searches for 2νECEC on 124 Xe have sought the simultaneous capture of two K-shell electrons (2ν2K) using a gas proportional counter with enriched xenon and large-volume liquid xenon (LXe) detectors with natural xenon as the target. In this paper, we report the results of an improved search for 124 Xe and 126 Xe 2ν2K events, using data from the XMASS-I detector. We analyze a new data set taken between November 2013 and July 2016. The total live time amounts to 800.0 days and the fiducial xenon mass was enlarged to 327 kg (containing about 311 g of 124 Xe). We developed a novel method for discriminating the 2ν2K signal from the β-ray background using LXe scintillation time profiles.

The XMASS-I detector
XMASS-I is a large single-phase LXe detector located underground (2700 m water equivalent) at the Kamioka Observatory in Japan [39]. An active target of 832 kg of LXe is held inside a pentakis-dodecahedral copper structure that hosts 642 inward-looking 2-inch Hamamatsu R10789 photomultiplier tubes (PMTs) on its approximately spherical inner surface at a radius of about 40 cm. The photocathode coverage of the inner surface is 62.4%. Signals from each PMT are recorded with CAEN V1751 waveform digitizers with a sampling rate of 1 GHz and 10-bit resolution.
The gains of the PMTs are monitored weekly using a blue LED embedded in the inner surface of the detector. The scintillation yield response is traced with a 57 Co source [40] inserted along the central vertical axis of the detector every week or two. Through measurements with the 57 Co source at the center of the detector volume, the photoelectron (PE) yield was determined to be ∼15 PE/keV for 122 keV γ-rays. The nonlinear response of the scintillation yield for electron-mediated events in the detector was calibrated over the energy range from 5.9 keV to 2614 keV with 55 Fe, 241 Am, 109 Cd, 57 Co, 137 Cs, 60 Co, and 232 Th sources. Hereinafter, this calibrated energy is represented as keV ee where the subscript 3/17 stands for the electron-equivalent energy. The timing offsets for the PMT channels owing to the differences in their cable lengths and the electronic responses were also traced by the 57 Co calibration.
The LXe detector is located at the center of a cylindrical water Cherenkov detector, which is 11 m in height and 10 m in diameter. The outer detector is equipped with 72 20-inch Hamamatsu H3600 PMTs. This detector acts as an active veto counter for cosmic-ray muons as well as a passive shield against neutrons and γ-rays from the surrounding rock.
Data acquisition is triggered if at least four inner-detector PMTs record a signal within 200 ns or if at least eight outer-detector PMTs register a signal within 200 ns. A 50 MHz clock is used to measure the time difference between triggers. One-pulse-per-second (1PPS) signals from the global positioning system (GPS) are fed as triggers for precise time stamping. The GPS 1PPS triggers are also used to flash the LED for the PMT gain monitoring.

Expected signal and simulation
The process of 2νECEC on 124 Xe is If two K-shell electrons in the 124 Xe atom are captured simultaneously, a daughter atom of 124 Te is formed with two vacancies in the K-shell and this atom relaxes by emitting atomic X-rays and/or Auger electrons. Our Monte Carlo simulations of the atomic de-excitation signal are based on the atomic relaxation package in Geant4 [41]. On the assumption that the X-rays and Auger electrons emitted in the 2ν2K event are like those generated by two single K-shell vacancies, the signal simulation begins with two Te atoms with a single K-shell vacancy. In such a case, the total energy deposition is given by twice the K-shell binding energy of Te (2K b = 63.63 keV). On the other hand, the energy of the two electron holes in the K-shell of 124 Te is calculated to be 64.46 keV [42], which only varies by 0.8 keV.
Since the energy resolution of the 2ν2K signal peak is estimated to be 3.2 keV after all the detector responses mentioned below are accounted for, we judge that this difference is negligible in our analysis. The results actually do not change even if the peak position of the simulated signal is artificially shifted by this amount. According to the simulation, 77% of 2ν2K events emit two K-shell X-rays, 21% of events emit a single K-shell X-ray, and the remaining 1.6% of events emits no K-shell X-ray. These probabilities are consistent with those expected from the fluorescence yield for the K-shell of Te, ω K = 0.875 [43]. Auger electron cascades are also simulated. The energy deposition from the recoil of the daughter nucleus is ∼30 eV at most, which is negligible. Simulated de-excitation events are generated uniformly throughout the detector volume. The nonlinearity of scintillation yield is accounted for using the nonlinearity model from Doke et al. [44] with a further correction obtained from the γ-ray calibrations. The absolute energy scale of the simulation is adjusted at 122 keV. The time profile of the scintillation is also modeled based on the γ-ray calibrations [45]. Propagation of scintillation photons in LXe is also simulated. Optical parameters of the LXe such as absorption and scattering lengths for the scintillation are tuned by source calibration data at various positions. The group velocity of the scintillation light in the LXe follows from LXe's refractive index (∼11 cm/ns for 175-nm light [46]). Charge and timing responses of the PMTs are also modeled in the simulation based on the calibrations with the LED and the γ-ray sources. Finally, waveforms

Data set
The data used in the present analysis were collected between November 20, 2013 and July 20, 2016. The data set was divided into four periods depending on the detector conditions at that time as summarized in Table 2. Period 1 started two weeks after the introduction of LXe into the detector. At the beginning of the run, we observed neutron-activated peaks from 131m Xe and 129m Xe that were created when the LXe was stored outside the water shield. We also performed the 252 Cf calibration data collection twice in this period. Runs within 10 days after each calibration were excluded from the data set. We ended period 1 60 days after the second 252 Cf calibration because the 131m Xe and 129m Xe peaks caused by the 252 Cf disappeared. Period 2 then ran until the continuous gas circulation at a flow rate of ∼1.5 L/min with a getter purifier was introduced. Period 3 was ended so that the xenon could be purified by vaporizing xenon once to remove possible non-volatile impurities dissolved in LXe. During the purification process, LXe was extracted from the detector, and therefore, xenon was exposed to and activated by thermal neutrons outside the water shield. The purification process took 7 days and period 4 started immediately after completing the introduction of LXe into detector. We selected periods of operation under what we call normal data taking conditions with a stable temperature (172.6-173.0 K) and pressure (0.162-0.164 MPa absolute) of the LXe in the detector. After further removing periods of operation that include excessive PMT noise, unstable pedestal levels, or abnormal trigger rates, the total live time became 800.0 days.

Event reduction and classification
The event-reduction process comprises four steps: pre-selection, the fiducial volume selection, 214 Bi identification, and particle identification.

Pre-selection
Pre-selection requires that no outer-detector trigger is associated with an event, that the time elapsed since the previous inner-detector event (dT pre ) is at least 10 ms, and that the standard deviation of the inner-detector hit timing distribution in the event is less than 100 ns. The last two requirements remove events caused by after-pulses in the PMTs following bright events. The dT pre cut eliminates events by chance coincidence at a probability of 3.0% on average, which was estimated from the fraction of the GPS 1PPS events that are rejected 5/17 by this cut. The chance coincidence probability is counted as dead time, and the live time mentioned above is obtained after subtracting this dead time.

Fiducial volume selection
To select events that occurred within the fiducial volume, an event vertex is reconstructed based on a maximum-likelihood evaluation of the observed light distribution in the detector [39]. We select events whose reconstructed vertex has a radial distance of less than 30 cm from the center of the detector. The fiducial mass of natural xenon in that volume is 327 kg, containing 311 g of 124 Xe and 291 g of 126 Xe.

214 Bi identification
222 Rn emanates from the detector's surface and contaminates the LXe within the detector. Thus, its daughters, 214 Bi and 214 Pb, become one of the major sources of β-ray background. 214 Bi can be tagged using the 214 Bi-214 Po delayed coincidence (T 1/2 = 164 µs) and is used as a good control sample of pure β-ray events in the relevant energy range. To remove the 214 Bi events, events whose time difference from the subsequent event (dT post ) is less than 1 ms are rejected from the 2ν2K signal sample. This cut reduces the 214 Bi background by a factor of ∼70, while consequently discarding only 0.4% of all other events. The counterpart sample, i.e. events with 0.015 ms < dT post < 1 ms, is referred to as the 214 Bi sample and is used to constrain the 214 Bi and 214 Pb backgrounds.

Particle identification
The scintillation time profiles of LXe can be used for particle identification. We use them to eliminate both the α-ray and β-ray backgrounds from the 2ν2K signal sample.
After the fiducial volume selection, the largest source of background in the relevant energy range are β-rays coming from radioactive impurities within the LXe. 2ν2K or γ-ray events can be discriminated from these β-ray events by utilizing the energy dependence of the scintillation decay time for electron-induced events. The scintillation decay time increases from 28 ns to 48 ns as the kinetic energy of an electron increases from 3 keV to 1 MeV as summarized in Fig. 3 of [45]. In the case of the 2ν2K or γ-ray events, the X-ray or γray is converted into multiple low-energy electrons in the LXe; this shortens the effective scintillation decay time by a few ns from that of an event caused by a single electron with the same deposited energy. Especially, events caused by 2ν2K or γ-rays with energy close to twice of the K-shell binding energy are easily-distinguishable from the β-ray events by this effective scintillation decay time.
The particle-identification parameter βCL is formulated as follows. First, waveforms in each PMT are decomposed into single-PE pulses using the single-PE waveform template [45]. scintillation time profile. The variable βCL is defined as where n is the total number of single-PE pluses after truncating the first 20 ns, P = n−1 i=0 CL i , and CL i (i=0, 1, 2, . . . , n − 1) is the CL of each pulse timing under the assumption that the event is caused by a β-ray. The probability-density function of the pulse-timing distribution for a β-ray event including its energy and position dependences is modeled from measurements in the 214 Bi data sample over the energy range between 30 and 200 keV ee . This formula is in general used to combine p-values from a set of independent tests for a certain hypothesis [47]. Figure 2 shows distributions of the variable βCL for the 214 Bi sample in the energy range from 30 to 200 keV ee along with the 241 Am 59.5 keV γ-ray events. While the β-ray events in the 214 Bi sample are distributed between 0 and 1, the distribution of the 59.5 keV γ-ray events in the 241 Am sample peaks at βCL = 0. Events with βCL less than 0.05 are classified as the β-depleted sample, and the rest is referred to as the β-enriched sample. When selecting events with βCL less than 0.05, 42% of the 2ν2K signal events are selected, while only 6% of the β-ray events from the 214 Bi decay in this energy range are selected. Thus, the signalto-noise ratio is improved by a factor of 7 by this selection. The cut position is tuned based on the simulated data to maximize the sensitivity for the 2ν2K signal.
α-ray events often occur in the grooves of the inner surface of the detector, so that only some of their energy is detected. These events are sometimes incorrectly reconstructed within the fiducial volume [48]. Above 30 keV ee , α-ray events can be clearly separated from β-ray or γ-ray events using the scintillation decay time. Therefore, the waveforms from all PMTs are summed up to form a total waveform of the event after correcting for the relative gain and timing of each PMT. Then, the falling edge of that total waveform is fitted with an exponential function to obtain the decay time for each event. Events with fitted decay times of less than 30 ns are deemed α-ray events and are rejected. Figure 3 shows the energy spectra plotted after each event-reduction step for the observed data and the simulated 2ν2K sample. For the simulated 2ν2K sample, T 1/2 = 4.7 × 10 21 years is assumed. 8/17 6. Spectrum fitting 6.1. Chi-square definition To extract the 2ν2K signal from the observed data, the energy spectra for the β-depleted samples, β-enriched samples, and 214 Bi samples are simultaneously fitted to the expected signal and background spectra. The energy range from 30 to 200 keV ee is used for fitting. The chi-square value is defined as where n data ijk , and n exp ijk (p l ) are the observed and expected number of events in i-th subsample and j-th period and k-th energy bin, respectively. N sample = 3, N period = 4, N bin = 85, and N sys are the number of sub-samples, periods, energy bins, and constrained systematic parameters, respectively. p l and σ l are a scaling parameter for the nominal value and its relative error, respectively.

Expected background
We consider three types of backgrounds: radioactive isotopes (RIs) in the LXe, neutron activation of xenon, and external backgrounds.
For the internal RIs, the 222 Rn daughters ( 214 Bi and 214 Pb), 85 Kr, 39 Ar, 14 C, and 136 Xe are considered. The 214 Bi activity during each period is determined from the fitting to the 214 Bi sample and the 214 Pb activity in each period follows from this 214 Bi activity as both originate from 222 Rn. While 85 Kr decays by β-decay (Q β = 687 keV, T 1/2 = 10.8 years) predominantly into the ground state of 85 Rb, 0.434% of their decays go into the 514-keV excited state of 85 Rb followed by a nuclear relaxation γ-ray (T 1/2 = 1.014 µs). 85 Kr contamination in the detector is measured to be 0.26 ± 0.06 mBq by the coincidence of β-ray and γ-ray events. In this analysis, the 85 Kr activity in each period is fitted with this constraint. We have also found argon contamination in the xenon through measurements of the sampled xenon gas using gas chromatography-mass spectrometry (GC-MS). The argon is thought to have adsorbed to the detector material when we conducted a leakage test of the LXe chamber using argon gas in 2013. 39 Ar undergoes β-decay (Q β = 565 keV, T 1/2 = 269 years). By comparing the energy spectra for periods 2 and 3 of the data set, we found a reduction in event rate below ∼150 keV ee in the β-enriched sample. The difference in the energy spectra between two periods is consistent with the β-decay of 14 C (Q β = 156 keV, T 1/2 = 5730 years). Hence, we assume that impurities containing carbon were reduced by gas circulation through the getter although its chemical form is not known. Finally, natural xenon contains 136 Xe with an isotopic abundance of 8.9%, and 136 Xe undergoes 2νββ decay (Q ββ = 2.46 MeV, T 1/2 = 2.2 × 10 21 years [49,50]).
Although the LXe detector is shielded against environmental neutrons by water, some of the detector components such as the cable feed-through box, calibration system, and cryogenic system lie outside the water shield and are filled with xenon gas [39]. The volume of xenon gas outside the water shield is estimated to be 2.6 × 10 5 cm 3 at the standard temperature of 273.15 K and pressure of 10 5 Pa. This xenon is activated by thermal neutron capture and returned to the LXe in the detector. The resulting 13 RIs, 125 Xe, 125m Xe, 125 I, 9/17  considered. Their activities are calculated based on the isotopic abundance of xenon and the cross sections of thermal neutron capture. Among those isotopes, 125 I is the most considerable background in this analysis. The 125 I is produced from 125 Xe and 125m Xe created by thermal neutron capture on 124 Xe with a total cross section of 165±11 barn [51]. 125 I decays by 100% electron capture via an excited state of 125 Te into the ground state of 125 Te with a total energy deposition of 67.5 keV. The flux of thermal neutrons (E < 0.5 eV) in the Kamioka mine has been measured to be (0.8-1.4)×10 −5 /cm 2 /s [52,53]. In this analysis, the thermal neutron flux during each period is fitted under the constraint of these measurements. 125 I, 131m Xe, and 133 Xe are the main RIs relevant to this analysis and the entire energy range of the beta-depleted samples has the power to constrain the thermal neutron flux. Activations of xenon by neutrons emitted from the (α, n) reaction or spontaneous fission in the detector material is negligible. In addition, occasional neutron activations of the LXe appear in periods 1 and 4 due to the 252 Cf calibration and the purification works. The data taken just after the 252 Cf calibration, which were excluded from this analysis, showed clear event rate increases in the energy range between 30 and 200 keV due to 131m Xe and 133 Xe. To accommodate these backgrounds, we introduce additional quantities of 131m Xe and 133 Xe in the fitting. For the external backgrounds, a detailed evaluation of radioactive backgrounds from each detector material has been conducted previously [48]. In the present data set, a small contribution of γ-ray backgrounds from impurities in the PMTs is expected. 238 U, 232 Th, 60 Co, and 40 K are considered, and the uncertainties in their activities are accounted for in the fitting. 10/17

Systematic uncertainties
Systematic uncertainties in the background yields, exposure, event selections, and energy scales are considered in the fitting as listed in Table 3. The upper part of Table 3 summarizes the systematic parameters used to determine the activities of RI backgrounds in the spectrum fitting. The 85 Kr activity, thermal neutron flux, and γ-ray backgrounds from the PMTs are constrained by the external measurements as described in the previous section since this spectrum fitting does not have the sensitivity for an independent evaluation.
Isotopic composition of the LXe was measured with a mass spectrometer, and the result was consistent with that of natural xenon in air [37]. The uncertainties in the measurement, ±8.5% for 124 Xe and ±12% for 126 Xe, are treated as a systematic error. The uncertainties in the LXe density and the detector live time are negligible. The uncertainties in the event selections and energy scales are estimated from comparisons between data and simulated samples for the 241 Am (59.5 keV γ-ray) and 57 Co (122 keV γ-ray) calibration data at various positions within the fiducial volume. The radial position of the reconstructed vertex for the calibration data differs from that for the simulated result by ±4.5 mm near the fiducial volume boundary, which causes ±4.5% uncertainty in the fiducial LXe volume. From the difference in the βCL distribution between the calibration data and simulated samples, the uncertainty in acceptance of the βCL cut for the γ-ray events is found to be ±30%. In the same manner, the uncertainty in the rejection power of the βCL cut for β-ray events is evaluated to be ±8.0% from the comparison of the βCL distributions for the 214 Bi sample in the energy range from 30 to 200 keV ee . By comparing the peak position of the γ-ray calibration data and simulated samples at various source positions and in different periods, the uncertainty in energy scale for the γ-ray events is estimated to be ±2.0%. Since we observe a small difference in the peak position of the γ-ray calibration data between the β-depleted and β-enriched samples, the energy scales for the β-depleted and β-enriched samples are treated independently.  Fig. 5. The 222 Rn concentration in LXe increased by ∼50% after gas circulation was initiated at the beginning of period 3. It is surmised that 222 Rn emanating from detector materials in the xenon gas volume mixes into the LXe by the gas circulation. An increase of the 39 Ar concentration in period 3 is thought to occur in the same manner. On the other hand, the 14 C concentration decreased with gas circulation. Figure 6 shows  LXe outside the water shield and caused by the 252 Cf calibrations. Increases in the 133 Xe yield in periods 1 and 4 are not significant compared with the increases in the 131m Xe yield. The fitted energy scales for the β-enriched and the β-depleted samples vary within ±2%, which is consistent with the evaluation before fitting.

Results and discussion
Closeup figures of energy spectra between 30 and 100 keV ee for the β-depleted samples are shown in Fig. 7. The peak found at 67.5 keV ee is attributable to the 125 I decay. The event rate of the 125 I decay is constrained by the thermal neutron flux. Figure 8 shows the normalized profile likelihood L/L max as a function of the inverse of the 124 Xe 2ν2K half-life, where L max is the maximum value of the likelihood. No significant excess over the expected background is found in the signal region. We calculate the 90% CL limit from the relation ξlimit 0 where ξ = 1/T 2ν2K 1/2 . This leads to T 2ν2K 1/2 124 Xe > 1 ξ limit = 2.1 × 10 22 years.   Fig. 7 Closeup figures of energy spectra between 30 and 100 keV ee for the β-depleted samples. The observed spectra (points) are overlaid with the best-fit 2ν2K signal and background spectra (colored stacked histograms). Colored histograms are the 2ν2K signal (red filled), 125 I (green hatched), 133 Xe (blue hatched), 14 C (orange filled), 39 Ar (magenta filled), 85 Kr (blue filled), 214 Pb (cyan filled), 136 Xe 2νββ (brown filled), and external backgrounds (gray filled).
The fact that we do not observe significant excess above background allows us to give a constraint on 2ν2K on 126 Xe in the same manner: T 2ν2K 1/2 126 Xe > 1.9 × 10 22 years (8) at 90% CL. Figure 9 shows a comparison of the experimental 90% CL exclusion limits on 124 Xe 2ν2K half-life overlaid with the theoretical calculations [21,[28][29][30][31][32] for comparison. The present result gives a lower limit stronger by a factor 4.5 over our previous result, and gives the most stringent experimental constraint reported to date. For the theoretical predictions, the reported 2νECEC half-lives are converted to 2ν2K half-lives, divided by the branching ratio for the two electrons being captured from the K-shell, P 2K =0.767 [54]. The lower and upper edges of the bands correspond to g A = 1.26 and g A = 1, respectively.
Note that the predicted half-lives will be longer if quenching of g A is larger. These experimental results rule out a part of the relevant range of the reported half-life predictions, and future experiments with multi-ton LXe targets will have improved sensitivity to further explore this parameter space.

Conclusion
We have conducted an improved search for 2ν2K on 124 Xe and 126 Xe using 800.0 days of data from XMASS-I. For this search, a novel method to discriminate γ-ray/X-ray or 2ν2K signals from β-ray backgrounds using LXe scintillation time profiles was developed. With spectrum fitting in the energy range from 30 to 200 keV ee , no significant 2ν2K signal appeared over the expected background. Therefore, we set the most stringent lower limits on the half-lives for these processes at 2.1 × 10 22 years for 124 Xe and 1.9 × 10 22 years for 126 Xe at 90% CL.