Search for neutrinoless quadruple beta decay of $^{136}$Xe in XMASS-I

A search for the neutrinoless quadruple beta decay of $^{136}$Xe was conducted with the liquid-xenon detector XMASS-I using $\rm 327\; kg \times 800.0 \; days$ of the exposure. The pulse shape discrimination based on the scintillation decay time constant which distinguishes $\gamma$-rays including the signal and $\beta$-rays was used to enhance the search sensitivity. No significant signal excess was observed from the energy spectrum fitting with precise background evaluation, and we set a lower limit of the half life of 3.7 $\times$ 10$^{24}$ years at 90$\%$ confidence level. This is the first experimental constraint of the neutrinoless quadruple beta decay of $^{136}$Xe.


Introduction
In spite of a great success of the standard model (SM) in particle physics, the nature of the neutrino is not yet understood thoroughly.If the neutrino is a Majorana particle, processes which violate the lepton number (L) by two (∆L = 2) can take place.The neutrinoless double beta decay (0νββ) is one of the ∆L = 2 processes beyond the SM.The observation of the 0νββ would tell us the Majorana nature of the neutrino.It is also linked to the seesaw mechanism, which explains the extremely light neutrino mass, and baryon number asymmetry in the universe via leptogenesis [1].Although a number of 0νββ experiments using a variety of candidate nuclei have been performed, no evidence has been achieved so far [2,3,4,5,6].
On the other hand, Heeck and Rodejohann proposed that even if the neutrino is a Dirac particle, a decay which violates L by four (∆L = 4) can occur by adding three right-handed neutrinos in the SM [7].This ∆L = 4 process is naturally linked to the light Dirac mass terms of neutrinos [8], CP violation [9] and leptogenesis [10].The neutrinoless quadruple beta decay (0ν4β) (A, Z) is one of the ∆L = 4 processes.Here A and Z are the atomic mass number and atomic number, respectively.It is theoretically predicted that only three candidate nuclei, 150 Nd, 136 Xe, and 96 Zr, can undergo this process.The NEMO-3 experiment reported the result of a search for 150 Nd 0ν4β decay [11], but there is no experimental search with either 136 Xe or 96 Zr 0ν4β decay so far.The Q-value of the 136 Xe 0ν4β (Q 0ν4β ) is 79 keV [12] 9 .The XMASS-I accumulated more than 2 years' data at a low background (BG) rate of O(10 −4 ) counts/keV/day/kg in the energy region of interest [13,14,15] under a stable temperature and pressure of the liquid xenon (LXe) [16].These data are suitable for a 136 Xe 0ν4β decay search with a high sensitivity.In this paper, a first search for the 0ν4β decay in XMASS-I with a total exposure of 327 kg × 800.0 days is reported.

XMASS-I detector
The XMASS-I detector is located in the Kamioka mine under 1,000 m of rock, corresponding to 2,700 meter water equivalent.The inner detector (ID) contains 832 kg of LXe inside a pentakis-dodecahedron shape oxygen free high conductivity copper structure with an inscribed radius of about 43 cm.Scintillation light from the LXe is detected by 630 2-inch hexagonal R10789 photomultiplier tubes (PMTs) and 12 2-inch cylindrical R10789Mod PMTs with a total photocathode coverage of 62.4% of the detector's inner surface 10 .Signals from the 642 ID PMTs were recorded by CAEN V1751 waveform digitizers with a sampling rate of 1 GHz.The radioactivity of the R10789 PMT is summarized in Ref. [17].In order to reduce external γ-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 of the detector can be found in Ref. [18].
The individual ID PMT gains were continually monitored by blue LEDs embedded in the inner surface of the ID PMT holder.The LEDs are flashed by one-pulse-per-second signals from the global positioning system.Calibration data were taken every one or two weeks by inserting a 57 Co source along the detector's z-axis to monitor the stability of the optical parameters described in Ref. [13,14,16,19].The γ-ray and X-ray calibration sources of 57 Co, 241 Am, and 55 Fe were used for the measurement of the scintillation decay time constant of electron recoils in LXe [20], the evaluation of the scintillation efficiency of the detector, and energy calibration.Since the energy calibration was performed with γ-ray and X-ray sources and the visible energy for the same deposited energy depends on the particle, the electronequivalent energy unit keV ee is used in this analysis to represent the event energies.

Signal simulation
In a neutrinoless quadruple beta decay of 136 Xe, four β-rays with a total energy of 79 keV are emitted simultaneously [12].Since these β-rays are expected to deposit all their energy in the liquid xenon at each decay, a monochromatic peak is expected in the energy spectrum.The event rate R for a given half life of the decay (T 0ν4β ) is calculated as where M is the mass of the LXe in the fiducial volume (327 kg), N a is the Avogadro constant, A is the atomic mass of 136 Xe atom, and k is the natural abundance of 136 Xe (0.089).The momentums of the four electrons from each 0ν4β decay were simulated by DECAY0 event generator [21].In the DECAY0 event generator, the four-dimensional energy distribution for a single electron ρ 1 can be found as where e i , t i , and p i are the total energy, kinetic energy and momentum of the i-th electron, respectively.t 0 is the total energy available in the 4β process.F (t, Z) is the Fermi function [22].The scintillation yield from the four electrons and the detector response were simulated by the XMASS Monte Carlo (MC) simulation based on Geant4 [23].In the XMASS MC, the Doke model [24] with a further correction based on the total energy deposition of the γ-ray calibration data was used to take the energy-dependent scintillation photon yield into account [25].Figure 1 (left) shows the expected energy spectrum of the 0ν4β signal.The expected peak energy reconstructed from the observed light is seen at an energy 4% higher than 79 keV.This is because the energy of each electron is lower than 79 keV.According to the Doke model [24], these low energy electrons are expected to have high light yields.So the number of the ob-served total photoelectrons from four electrons with 79 keV total energy is larger than that of a single 79 keV electron.The light yield of the 0ν4β signal has uncertainty on both of the increasing and decreasing sides.It could decrease 3.0% at maximum due to the scintillation model difference between the Doke model and the NEST model at a zero-electric-field [26].It could increase 8.6% at maximum due to the enhancement of the recombination probability [27] since these four electrons are generated from one nucleus and their energies are deposited in a very small volume.The uncertainty of the light yield of the 0ν4β signal is treated as a constrained parameter in the spectrum fitting described in Sec. 5.

Data and event classification
A search for the 0ν4β decay was carried out with the data accumulated from November 2013 to July 2016.The data set was same as that for our two-neutrino double electron capture on 124 Xe and 126 Xe [13] analysis and the WIMP-129 Xe inelastic scattering search [14].During this period, the detector condition, specifically the temperature and pressure of LXe, was kept stable [16].The data set was divided into four subsets (subset 1-4) depending on the Xe gas circulation status and the BG rate due to the neutron activation.The data within about 10 days after the neutron calibrations was not used since the BG rate due to the activation was high.In the subset 1, the rate of the BG events from neutron-activated Xe isotopes was higher than other subsets because the data in this subset were taken only two weeks after the Xe installation and 252 Cf neutron calibrations were performed twice.The subset 2 started 60 days after the second 252 Cf calibration.Neutron-activated peaks from 131m Xe and 129m Xe due to the 252 Cf calibration disappeared and the rates of these peaks were stable in the subset 2. Gas circulation for the Xe purification was started at the beginning of the subset 3.During the circulation, gas Xe was extracted from the detector, passed through a hot getter, condensed into liquid and returned to the detector.The data in the subset 4 was taken after the removal of non-volatile impurities.In this removal process, the Xe was recovered from the detector to the reservoir tank located outside the detector.Then the evaporated xenon in the reservoir tank was returned to the detector through the hot getter.
The data went through several event selections.First, the event is triggered only by ID.Then the event should have both the time difference from the previous event more than 10 ms and the root mean square of the hit tim-ing in the event less than 100 ns.These selections are to remove the event caused by afterpulse following bright events.For the rejection of external γrays, we apply the fiducial volume selection using the position reconstruction of the events' vertices.In this reconstruction, measured light distributions on the PMTs were compared with the ones produced by the MC [18].The events whose vertices were reconstructed inside the fiducial volume within 30 cm radius, were selected.This fiducial volume contains 327 kg (29 kg) of LXe ( 136 Xe).The internal BG of 222 Rn daughter nuclei were quantified by counting 214 Bi events taking the coincidences of 214 Bi-214 Po chain-decays.To take the coincidence, the time to the next event (dT next ) selection was used.By this coincidence, 99.6% of the 214 Bi events can be tagged by selecting events with 0.015 ms < dT next < 1 ms since the half-life of 214 Po is 164 µs.The tagged events will be referred to as the 214 Bi samples.In this sample, 0.4% of non-214 Bi events were contaminated due to the accidental coincidence.
Further event selections were applied on non-tagged event to estimate the abundance of other BG and the signal.α-rays are typical BG events observed in the non-tagged events.A small fraction of scintillation light from α events would leak into the sensitive region via small gaps of the detector structure and would be detected by the PMTs.These events can be removed using the scintillation time constant obtained by the fitting of the summed-up waveforms with an exponential function.Those with scintillation time constants shorter than 30 ns were classified as α events.This scintillation time selection eliminated almost all α events in the energy range above 30 keV ee , while more than 99.9% of the signal events were retained.
The events which survived the α-event reduction were further classified into β-depleted and β-enriched samples using the particle identification selection.The particle identification uses the difference of the scintillation time profiles between β-rays and γ-rays.The scintillation time constant of LXe has two components; τ s = 2.2 ns, and τ t .The smaller energy β-rays have shorter τ t than the larger energy β-rays do [20].γ-rays have shorter τ t than β-rays with the same energy do since the γ-rays are converted into lower energy electrons via Compton and photoelectric effect.For example, a 79 keV β-ray has τ t of 38.5 ns.On the other hand, 79 keV γ-ray typically creates one 44 keV electron, one 25 keV electron, and a few keV electrons by photoelectric effect.44 keV and 25 keV β-rays have τ t of 34.4 ns and 31.6 ns, respectively.These values are evaluated from Table 2 in Ref. [20].This classification makes the evaluation of the abundance of β-rays and γ-rays BGs easier by improving the significance of the γ-ray events.Also the significance of the 0ν4β decay's signal is enhanced due to the shorter scintillation time of the signal since each electron has less energy than single β-ray events with the Q 0ν4β energy.For example, 20 keV (a quarter of the Q 0ν4β energy) β-ray has τ t of 30.3 ns.
We introduce a particle identification parameter βCL defined as where N is the number PMTs with pulses; t k is the timing of k-th pulse; E is the event energy, and CDF β (E, t) is the cumulative distribution function (CDF) for finding a pulse at time t in a β-event of energy E discussed in [13].
By definition, βCL distribution becomes uniform between 0 and 1 for βray events.On the other hand, a peak at 0 is expected for particles with shorter decay times than β-rays such as γ-ray and 0ν4β signal.α events also distribute around 0 in βCL distribution, although such events were already removed by the scintillation time selection.Thus, the ratio of the γ-ray and the signal events to β-ray events can be enhanced in the β-depleted sample.The probabilities that β-ray, γ-ray including the 0ν4β signal are classified as β-depleted sample are referred to as β-ray mis-identification probability (β mis-ID), and γ acceptance, respectively.We set the threshold of the βCL (βCL th ) at 0.05 for this analysis, following Ref.[13].The events with βCL less than 0.05 and greater than 0.05 are classified as β-depleted and β-enriched samples, respectively.MC results showed that 53±16% of 0ν4β signal events would remain after selecting events βCL less than 0.05, with 94±2% of βrays events from 214 Bi events are rejected by this selection, corresponding to a signal-to-noise ratio improvement by a factor 9. Here, the errors (±16%, ±2%) correspond to the uncertainty of the βCL selection referred to as the γ acceptance and the β mis-ID in Sec. 5, respectively.The signal MC and the data spectra after each treatment are shown in Figure 1.

Energy spectrum fitting
The events were classified into a β-depleted sample, a β-enriched sample, and a 214 Bi sample in the previous section.By fitting the energy spectra of these three samples simultaneously, we estimated the abundance of the signal and the BG.To estimate the activities of the BG, the energy range of Energy dependent as shown in Fig. 2 the fit was set from 30 to 200 keV ee .The energy bin's width is 2 keV ee .A χ 2 value is defined as where n MC ijk is the expected number of events including the BG MC and signal MC, and n data ijk is the number observed of events.The signal histogram was scaled by τ −1 0ν4β .Indices "i", "j", and "k" mean i-th sample (β-depleted, β-enriched, and 214 Bi), j-th subset, and k-th energy bin, respectively.Here, N sample = 3, N subset = 4, and N bin = 85.p const l (l = 1, 2, • • • , N sys = 42) and p free m are constrained parameters and free parameters, respectively.The systematic uncertainty for p const l is σ l .The systematic uncertainties are summarized in Table 1.
The BGs in this study are the radioactive isotopes (RIs) in the PMTs ( 40 K, 60 Co, 232 Th and 238 U), RIs uniformly distributed in the LXe ( 214 Pb, 214 Bi, 136 Xe 2νββ, 124 Xe 2ν double electron capture, 85 Kr, 39 Ar and 14 C) and RIs originating from thermal neutron captures ( 133 Xe, 131m Xe and 125 I).The constraints for RIs in the PMTs and 85 Kr were determined based on the results of the BG study in XMASS [19].The constraints for the abundance of 214 Pb and 214 Bi, 124 Xe 2ν double electron capture, 39 Ar, and 14 C were not given and their abundances were determined by the fitting.Here, 39 Ar is thought to be in argon gas which was used for a leakage test of the detector before starting data taking, and to have been absorbed in the detector material [13].To account for the 133 Xe and 131m Xe BG amount related to the special work of neutron calibration and purification work, we introduce the free parameters of the abundance of 133 Xe and 131m Xe in the fitting.The 2ν2β BG events of 136 Xe were constrained by the result of KamLAND-Zen experiment [28].Thermal neutron flux measurements [29,30] gave constraints of neutron-induced RI of 125 I. Since detector condition affects the amount of 85 Kr and BG from thermal neutron flux, we considered the subset dependence of these two BGs.
The uncertainty of the isotopic abundance of 136 Xe was evaluated from a modified VG5400/MS-III mass spectrometer measurement at the Geochem-ical Research Center, the University of Tokyo [31].The abundance was consistent with that of the natural xenon in the air.The uncertainty of ±1.3% in the measurement was treated as the systematic uncertainty of the abundance of 136 Xe isotope.The comparison of the MC and the data of γ-rays from 241 Am's 59.5 keV γ-rays and 57 Co's 122 keV γ-rays was used for the evaluation of the uncertainty of the fiducial volume (±4.5%), the energy scale for β-depleted and β-enriched sample (both ±2%), and γ acceptance (±30%) [13,14].Since energy scale varies depending on the detector condition, we introduced the energy scale parameter in each subset in Eq. 5.
The dead PMT effect (7% ± 14% for 30 < E < 35 keV ee and 19% ± 16% for 35 < E < 40 keV ee region), accounts for the mis-reconstructed events increase due to dead PMTs discussed in Refs.[15,19,32].The uncertainty of the β mis-ID was evaluated by the comparison of the data and MC of 214 Bi events selected by the 214 Bi-214 Po coincidence.Using the β-ray's continuous spectrum over the energy region of this analysis, the energy-dependent uncertainties were obtained.The energy region from 30 to 200 keV ee was divided into 17 bins.Probabilities of β mis-IDs for data and MC of 214 Bi delayed coincidence events were compared in each bin.To correct the difference of β mis-ID between the data and MC, β-ray BG MC histograms were scaled energy-dependently in the fitting using the β mis-ID ratio obtained from 214 Bi events between the data and MC shown in Fig. 2. The uncertainty of the light yield of the 0ν4β mentioned in Sec. 2 (-3.0%, +8.6%) is treated as a constrained parameter without the penalty χ 2 term in Eq. 5.

Results
The energy spectra of three samples (β-depleted, β-enriched, and 214 Bi samples) were fitted with the MC spectra of the signal and BG simultaneously so that the signal abundance and the BG abundances were determined at the same time.No significant signal excess over the expected BG was found and the null 0ν4β signal case gave the best-fit with χ 2 /ndf = 1096/997.Thus, the 90% confidence level (CL) lower limit on the T 0ν4β was calculated by the following equation.
The calculated 90% CL lower limit is 3.7 × 10 24 years.Measured βdepleted samples' spectra (black histograms) and best-fit BG spectra (color  histograms) are shown in Figure 3.It should be noted that 0ν4β spectra with the size of the 90% CL lower limit of the half life (red) are added on top of the best-fit BG spectra to illustrate the contribution to the BG ones.The LXe purification started at the beginning of subset 3 decreased the activity of 14 C (orange) as time proceeds.The cause of the increasing of 39 Ar (magenta) was presumably due to the emanation from the inner structure of the detector [13].This is the first experimental constraint on the half life of 0ν4β decay of 136 Xe and the longest half life limit of the 0ν4β decay.

Conclusion
A search for the neutrinoless quadruple beta decay of 136 Xe was conducted using 327 kg×800.0days of XMASS-I data.The particle identification based on the difference of the scintillation time constant between γ-rays and β-rays enhanced the sensitivity.The particle identification parameter (βCL) and timing-coincidence (dT next ) was used to create three separate subsets of the data (β-enriched, β-depleted, and 214 Bi sample).These spectra were fit with the BG MC + 0ν4β signal MC.In this analysis, no significant excess over the expected BG was found and a 90% CL lower limit of the half life of 3.7 × 10 24 years for the neutrinoless quadruple beta decay of 136 Xe was set.This is the first experimental constraint on the neutrinoless quadruple beta decay of 136 Xe.

Figure 1 :
Figure 1: Energy spectra of the signal MC (left) and measured data (right) at various selection stages.Histograms after the pre-selection (black), the fiducial volume selection, dT next selection for removing 214 Bi events, the scintillation time constant selection for removing α events (purple), and particle identification selection for selecting β-depleted sample (βCL < 0.05) (red) are shown.The right-hand figure depicts the entire energy region used for the analysis (30 − 200 keV ee ), while the left-hand figure shows a magnified view of the 40 − 120 keV ee region to make it easy to observe the 0ν4β signal.

Figure 2 :
Figure 2: Scaling factor of β-ray MC for the correction of β mis-ID events.Black points and their error bars represent the means and errors calculated from the difference between MC and the data.In the fitting process, β-ray MC spectra were scaled by factors "(mean)+ (p const l − 1)(error)" depending on the energy.The red points are the scale factor for the best fit and 1σ error bars.

Table 1 :
Summary of the systematic uncertainty used as p const l in the spectrum fitting.