Searching for Neutrino-less Double Beta Decay of $^{136}$Xe with PandaX-II Liquid Xenon Detector

TL;DR

This paper reports the first search for neutrinoless double beta decay of $^{136}$Xe using the PandaX-II liquid xenon detector, setting new limits on decay half-life and neutrino mass.

Contribution

It presents the first neutrinoless double beta decay search results from a dual-phase xenon experiment, with optimized event selection and analysis.

Findings

No evidence of neutrinoless double beta decay was observed.

Lower limit on decay half-life is 2.4 x 10^{23} years.

Effective Majorana neutrino mass is constrained to less than 1.3-3.5 eV.

Abstract

We report the Neutrino-less Double Beta Decay (NLDBD) search results from PandaX-II dual-phase liquid xenon time projection chamber. The total live time used in this analysis is 403.1 days from June 2016 to August 2018. With NLDBD-optimized event selection criteria, we obtain a fiducial mass of 219 kg of natural xenon. The accumulated xenon exposure is 242 kg$\cdot$yr, or equivalently 22.2 kg$\cdot$yr of $^{136}$Xe exposure. At the region around $^{136}$Xe decay Q-value of 2458 keV, the energy resolution of PandaX-II is 4.2%. We find no evidence of NLDBD in PandaX-II and establish a lower limit for decay half-life of 2.4 $ \times 10^{23} $ yr at the 90% confidence level, which corresponds to an effective Majorana neutrino mass $m_{\beta \beta} < (1.3 - 3.5)$ eV. This is the first NLDBD result reported from a dual-phase xenon experiment.

Equations

$$E=W\times (\frac{S1}{PDE}+\frac{S2}{\text{EEE}\times \text{SEG}}).$$

$$S2=S{2}_{\text{Bottom}}\times (T/B+1),$$

$$C=\frac{{\Gamma }_{0\nu }\times T\times N_A\times R\times M\times ϵ}{m},$$

$${\sigma }_{\text{ syst }}\left({\Gamma }_{0v}\right)=\sqrt{\sum _i\left({\sigma }_{\mathrm{add},i}^2+{\sigma }_{\mathrm{scaling},i}^2{\Gamma }_{0v}^2\right)}.$$

Full text

Kaixiang Ni(倪恺翔)1  Yihui Lai(赖奕辉)1  Abdusalam Abdukerim(阿布都沙拉木·阿布都克力木)1  Wei Chen(陈葳)1  Xun Chen(谌勋)1  Yunhua Chen(陈云华)2  Xiangyi Cui(崔祥仪)1  Yingjie Fan(樊英杰)3  Deqing Fang(方德清)4  Changbo Fu(符长波)1  Lisheng Geng(耿立升)5  Karl Giboni1  Franco Giuliani1  Linhui Gu(顾琳慧)1  Xuyuan Guo(郭绪元)2  Ke Han(韩柯)1,†  Changda He(何昶达)1  Di Huang(黄迪)1  Yan Huang(黄焱)2  Yanlin Huang(黄彦霖)6  Zhou Huang(黄周)1  Peng Ji(吉鹏)3  Xiangdong Ji(季向东)1,7  Yonglin Ju(巨永林)8  Kun Liang(梁昆)1  Huaxuan Liu(刘华萱)8  Jianglai Liu(刘江来)1,7,‡  Wenbo Ma(马文博)1  Yugang Ma(马余刚)4  Yajun Mao(冒亚军)9  Yue Meng(孟月)1  Parinya Namwongsa1  Jinhua Ning(宁金华)2  Xuyang Ning(宁旭阳)1  Xiangxiang Ren(任祥祥)1,8  Changsong Shang(商长松)2  Lin Si(司琳)1  Andi Tan(谈安迪)10  Anqing Wang(王安庆)11  Hongwei Wang(王宏伟)4  Meng Wang(王萌)11  Qiuhong Wang(王秋红)4  Siguang Wang(王思广)9  Xiuli Wang(王秀丽)8  Zhou Wang(王舟)10,12,1  Mengmeng Wu(武蒙蒙)3  Shiyong Wu(吴世勇)2  Jingkai Xia(夏经铠)1  Mengjiao Xiao(肖梦姣)10,12  Pengwei Xie(谢鹏伟)7  Binbin Yan(燕斌斌)11  Jijun Yang(杨继军)1  Yong Yang(杨勇)1  Chunxu Yu(喻纯旭)3  Jumin Yuan(袁鞠敏)11  Dan Zhang(张丹)10  Hongguang Zhang(张宏光)1  Tao Zhang(张涛)1  Li Zhao(赵力)1  Qibin Zheng(郑其斌)6  Jifang Zhou(周济芳)2  Ning Zhou(周宁)1  Xiaopeng Zhou(周小朋)5

(PandaX-II Collaboration)

3.40.-s, 95.35.+d, 14.60.pq

00footnotetext: $\dagger$ Corresponding author: ke.han@sjtu.edu.cn00footnotetext: $\ddagger$ Spokesperson: jianglai.liu@sjtu.edu.cn

1 Introduction

Neutrino-less double beta decay (NLDBD) is a hypothetical decay process, during which two neutrons in a nucleus decay to protons and emit two electrons but no neutrinos. This process is a direct evidence that neutrinos are their own anti-particles, i.e., Majorana Fermions. It also violates lepton number conservation since the emitted electrons are not accompanied by corresponding electron anti-neutrinos. NLDBD is clearly a physical process beyond the Standard Model (SM) of particle physics [1]. The SM-allowed double beta decay process with two electron anti-neutrinos emitted has been observed in a dozen or so isotopes, but NLDBD has not been convincingly observed yet.

There is a global effort to search for NLDBD with various candidate isotopes using different experimental approaches. For example, 76Ge is widely used in high purity germanium semiconductor detectors [2, 3, 4], while 136Xe is usually the target of monolithic liquid detectors [5, 6]. CANDLES-III [7], GERDA [2], CUORE [8], and KamLAND-Zen [5] experiments have recently established the most stringent half-life lower limit respectively for 40Ca, 76Ge, 130Te, and 136Xe decaying into the ground state of their daughters. Those limits are on the order of $10^{25}$ - $10^{26}$ years and show an extremely rare nature of NLDBD, if exists at all.

The PandaX-II detector [9] is a large dual-phase liquid xenon time projection chamber (TPC). The primary science goal is to directly detect possible collisions of dark matter candidate particle WIMP (Weakly Interacting Massive Particle) and xenon nuclei [9, 10, 11]. However, with a target mass of 580 kg of xenon and a natural isotopic abundance of 8.9%, it contains 51.6 kg of 136Xe and gives us the possibility to search for NLDBD with a reasonably large amount of isotope.

In this paper, we present the NLDBD search results with 403.1 days of data from the PandaX-II experiment. Section 2 will describe the PandaX-II detector and the data taking campaigns in year 2016-2018. In Section 3, we outline the data analysis procedure, with the energy reconstruction and event selections adjusted specifically for the high energy region of interest (ROI). In Section 4 and 5, we establish the background model from Monte-Carlo and fit the energy spectrum to identity possible NLDBD signal. We also discuss in details different contributions to systematic uncertainties. We summarize and give an outlook for NLDBD physics potentials with large liquid xenon detectors in the last section.

2 PandaX-II detector and data taking campaigns

Located at China Jin-Ping Underground Laboratory (CJPL) [12], the PandaX-II experiment has been taking physics data since 2016. The detector is a dual-phase xenon TPC [13], as shown in Fig. 2. Detailed description of the detector can be seen in [9]. The sensitive volume is defined by a cylinder of polytetrafluoroethylene (PTFE) reflectors of inner diameter 646 mm and two sets of photomultiplier tube (PMT) arrays at the top and bottom. The drift/extraction/amplification field in the TPC is exerted by the cathode, gate, and anode meshes, and is constrained by a set of copper shaping rings outside the reflector barrel. Particle going through the liquid generates scintillation photons and ionized electrons, the latter of which drift to the gas phase and generate electroluminescence photons. Light signals are used to reconstruct energy and position information of a particular event. Vetoing volume resides outside of the PTFE reflectors and is used to reject environment and detector background events originating from outside.

Multiple layers of passive shielding outside the TPC include copper, inner polyethylene, lead, and outer polyethylene (see Fig. 2). Bricks or panels of those materials are concatenated into octagonal shape. The gap between the innermost layer and the TPC vessel is flushed by nitrogen gas in order to suppress radon concentration inside. Ambient radiation such as gamma and neutron background, is attenuated effectively with the help of passive shielding.

\figcaption

Schematic drawing of the PandaX-II TPC. Key components such as top and bottom PMT arrays, PTFE reflectors, and electrodes are highlighted. The figure is adapted from [9].

\figcaption

Schematic drawing of the layered passive shielding system around the PandaX-II Detector. Thickness of different layers are labelled. The figure is adapted from [14].

In this paper, three sets of physics data from the past 3 years are used for NLDBD analysis, as shown in Table 2. The total accumulated exposure is 403.1 days. Run numbers in the table follow the convention used in dark matter direct detection physics analysis. Run 9 and 10 are the first two sets of physics runs for PandaX-II and have been previous published [10, 11]. Run 11 is a continuation of Run 10 with no change in run conditions. Run 11 data taking was undertaken smoothly from July 2017 to August 2018, interrupted only by a series of calibration runs, as shown in Fig. 2. The figure also shows that the drift electron lifetime, which is an indication for detector performance, was stable at above 400 µs most of the time. During this run, an additional 246.4 live-days of data was accumulated. This is the first time that we present an analysis with this new data set.

\tabcaption

PandaX-II data sets used for this NLDBD analysis.

Data set Begin End Live days

Run 9 Mar. 9, 2016 Jun. 30, 2016 79.6 d

Run 10 Apr. 22, 2017 Jul. 16, 2017 77.1 d

Run 11 Jul. 17, 2017 Aug. 16, 2018 246.4 d

\figcaption

Data taking history of Run 11. Other than 4 interruptions of calibration and detector performance study, we took data continuously from July 2017 to August 2018. The left Y axis is the exposure in the unit of kg$\cdot$yr. The dash horizontal line indicates the total exposure at 242 kg$\cdot$yr, with a corresponding fiducial mass of 219 kg. The right Y axis is the electron lifetime while drifting in the field cage.

3 Energy reconstitution and event selection at MeV region

PandaX-II data analysis at MeV level follows the workflow of dark matter analysis closely [15]. However, we developed new algorithms for energy reconstruction and new event selection criteria at MeV level specifically for NLDBD analysis. In this section, we summarize the analysis workflow and highlight the key differences.

3.1 PandaX-II data analysis workflow

PandaX-II data acquisition system (DAQ) records triggered waveforms of prompt scintillation light signal (commonly called $S1$) and electroluminescence light signal ($S2$). For each triggered signal, the readout window is 1 ms and the trigger point is placed in the middle of the signal window. Data are recorded by underground DAQ computers and immediately copied to a local computing farm aboveground at CJPL, where first level data preprocessing is performed in situ automatically.

The first step of data preprocessing chain process waveform of individual PMT. After the Single Photoelectron (SPE) correction, which we will describe later, we group these waveforms into individual pulses, known as hits, using a hit finding algorithm.

The second step is to group coincident hits of PMTs to construct light signals. Waveforms from individual pulses are summed up and the key parameters are calculated for the summed pulse. We then classify these signals into $S1$, $S2$ or noise using a decision tree based on the characterization of summed waveforms. A typical event may include one $S1$ signal and one or multiple $S2$ signals. Each $S2$ signal paired with the preceding $S1$ signal corresponds to one interaction vertex (usually called site) in the liquid xenon medium.

Finally we reconstruct the positions of sites in each event. The $Z$ direction follows electron drift direction and points upward while $X$ and $Y$ directions define the horizontal plane. The $X$ and $Y$ coordinates of each track are calculated from $S2$ signal. Multiple methods are used in this step, including center of gravity (COG), template matching (TM), and maximum likelihood finding using photon acceptance function (PAF) [10], the last of which gives best uniformity of event distribution of calibration source, and is used as standard position reconstruction algorithm in the following analysis. The $Z$ coordinate is calculated from electron drift velocity in the electric field and the drift time dt, defined as time difference between $S1$ and $S2$ signals. In Run 9 the drift velocity is $1.71\pm 0.01$ mm/µs with an electric field of 393.5 V/cm. In Run 10 and Run 11, the drift velocity dropped to $1.67\pm 0.01$ mm/µs due to the lowered drift field of 311.7 V/cm. The derived drift velocities are consistent with Ref. [16].

3.2 Energy Calibration and Reconstruction

The detector’s response at MeV level is calibrated using external 232Th source capsules. The source is deployed through the calibration loops inside the TPC outer vessel. A full absorption gamma peak from 208Tl at 2615 keV can be seen in both $S1$ and $S2$ signals. For MeV scale signals, a large $S2$ happens right below the top PMT array and the number of photons collected by a PMT may be beyond the linear response region. In rare cases when signal size exceeds the dynamic range of digitizers, the waveform is also clipped. Saturation of PMTs and digitizers deteriorates energy resolution. A group of saturated pulses with both PMT and digitizer saturation is shown in Fig. 3.2.

\figcaption

Examples of saturated waveform of PMTs. Different traces show waveforms from different top PMTs from a high energy event. The saturated PMT (black and yellow line) has asymmetric peak in its waveform, where the rising edge is steeper compared to the other ones. For PMT 10900, a typical PMT that experiences saturation problem, clipped peak is also observed as the digitizer saturates.

In our analysis, we used pulses from the bottom PMT array to reconstruct the energy of $S2$ signals to avoid the saturation problem. Saturation is not observed in the bottom PMT array, since $S1$ signal is usually much smaller than $S2$ and $S2$ signal is far away from the bottom array. We describe the whole process below and emphasis is given on last two steps, where pulse saturation has key impact.

SPE correction

Hit pulses are scaled by the individual SPE values of PMTs so that the area of pulse represents the number of photoelectrons (PE). This correction compensates possible PMT gain drift during detector operation.

Position uniformity correction

The number of PE each PMT collects for any specific event depends on the $X$, $Y$, and $Z$ position of the interaction vertex. The position dependence is mapped out with internal radioactive sources such as 83Kr or activated 127mXe. 127mXe emits 164 keV gamma rays and is better suited for position dependence correction at the high energy region. We divide the sensitive volume into small voxels to map out detector response at different position. In turn PMT response for any events is corrected with the mapping.

Secondary electrons from xenon ionization are recombined with or attached to the xenon medium in the drifting process. The electron loss is usually characterized by electron life time in the unit of µs. The measured $S2$ signal at liquid-gas interface is corrected with electron life time on a daily basis.

Energy reconstruction

To reconstruct the energy of an event, we follow the energy sharing formula used in dark matter analysis:

[TABLE]

Here W equals to 13.7 eV and is the average work function for a particle to excite one electron-ion pair in xenon [17]. $S1$ is the total detected PE number from both top and bottom PMTs, after event position correction. Since the saturated top PMT charge cannot be used, we replace $S2$ with bottom PMT PE number:

[TABLE]

where $T/B$ stands for ratio between the numbers of PE from the top and bottom PMT array. We expect that photon emission as well as collection is energy independent and set $T/B=2.3$, as derived in the dark matter analysis [10]. The other three parameters, photon detection efficiency (PDE), electron extraction efficiency (EEE), single electron gain (SEG), are derived by scanning the parameter space as follows. For each combination of PDE and EEE $\times$ SEG, a spectrum of 232Th calibration source is obtained and fitted. The combination of 7.4% and 22.75 PE/e yields the best energy resolution and is selected. The achieved best energy resolution is 4.2% at 2615 keV. The PDE and EEE $\times$ SEG are different from those used in dark matter analysis, which implies changes of detector response in this different energy region. The new set of values is verified with the anti-correlation of $S1$ and $S2$ bottom signals, as shown in Fig. 3.2.

Residual non-linearity correction

A spectrum with optimized energy resolution can be obtained with the above mentioned procedures. However, peaks could be identified slightly off the true energy of the characteristic gamma lines. Another correction is carried out to shift those peaks and correct the non-linearity of detector’s energy response. With a relatively large energy resolution, gamma lines on top of continuum background may not peak in the exact gamma energy. Monte carlo simulation is used to study such shifts. The detailed simulation process is discussed in Sec. 4. Peaks from 915 to 2615 keV are used to correct the residual non-linearity, as shown in Fig. 3.2. For each run, six peaks are used to fit a third-order polynomial function with zero intercept. Then energy of each event is corrected with these fit functions.

\figcaption

Anti-correlation between $S1$ signals (Top + Bottom PMTs) and $S2$ Bottom PMT signals in 232Th calibration data. The bright blob on the top-right is the 2615 keV gamma peak. The red line indicates the best anti-correlation between $S1$ and $S2$ signals, whose slope is calculated from the ratio of PDE and EEE $\times$ SEG.

\figcaption

Peak matching and residual nonlinearity of the three runs. The X axis is the detected peak position from 915 to 2615 keV before corrections. The Y axis is the expected peak position from Monte-Carlo simulation. Three-order polynomial functions are used to fit the data to determine the correction functions for corrections.

3.3 Event Selection Cuts

The original cuts from dark matter analysis are inherited but tuned to adapt to the high energy region in this study. We first apply the single-site cut, ROI energy cut, veto cut, and fiducial volume cut to reject a vast majority of background events. The signal quality cuts, including a few pulse shape cuts, $S1/S2$ band cut, top/bottom asymmetry cut, and position quality cuts, are used to further reject possible low quality signals. Here we briefly describe the cut criteria and more detailed definition can be found in [14, 15].

Single-site cut

We focus on single-site events from NLDBD and reject all events with multiple $S2$ signals, because in most cases the two electrons emitted from an NLDBD event deposit energy in one point-like site in liquid xenon. Monte Carlo shows that such selection accepts 93.4% of all the NLDBD events, and the inefficiency is dominated by the multi-site deposition due to the Bremsstrahlung. On the other hand, background events from gammas are predominately multi-site events, the cut reduces gamma background effectively. For example, in Run 10, 6.9 million out of 18.4 million of triggered events are tagged as single-site events.

ROI cut

In the NLDBD final spectrum fit, we use events with energy in the range of 2058 to 2858 keV, i.e. $\pm 400$ keV around the Q-value.

Veto cuts

The large majority of NLDBD events inside the field will not deposit any energy outside. On the contrary, penetrating particles such as muon may deposit energy in the veto region and is rejected with an offline coincidence analysis between veto signal and $S1$ signal.

Fiducial volume cut

Fiducial volume for NLDBD analysis is optimized according to the distribution of high energy events. Fig. Fiducial volume cut shows the spatial distribution of high energy events in the ROI. The X axis shows square of event radius. The Y axis is the negative drift distance, which is equivalent to the relative vertical coordinate in space. Number of events per bin remains relatively constant along the radius direction since xenon self-shielding is not as effective for high energy events. At the edge of the sensitive volume, an abnormal decrease of event density is observed, which is due to malfunctioning of PAF position reconstruction algorithm with saturated PMT pulses at the detector’s edge. To be conservative, we force $r^{2}<80000\,\textrm{mm}^{2}$. More background events concentrate at the top of the sensitive volume where the top PMT array, stainless steel flanges and pipes contribute significantly. Therefore, the fiducial volume cut in vertical direction is not symmetric. The vertical cut is within -233 and -533 mm, optimized by comparing fiducial mass and the number of background events under different ranges. The final fiducial region is shown in the red rectangle in Fig. Fiducial volume cut and the fiducial mass within is 219 kg. The fiducial volume cut is the strongest cut in the event selection, reducing large amount of background events by the self-shielding effect of liquid xenon.

We estimate the uncertainty of the fiducial mass with position reconstruction uncertainties. The uncertainty of position reconstruction in the radial direction is 21.5 mm and given by the mean difference between of PAF and TM method for events in the ROI. The uncertainty in the vertical direction is 3 mm, mainly from the uncertainty of drift velocity measurement. Overall, we obtained an uncertainty of 32 kg for our fiducial mass.

\figcaption

Spatial distribution of ROI event in the PandaX-II TPC. The Y axis is in the electron drift direction with the liquid and gas interface at 0 mm. Single scatter cut, quality cut, and veto cut are added. Fiducial volume used in the NLDBD analysis is marked as red rectangle.

Pulse shape cuts

Two pulse shape cuts are used to reject pulses with positive overshoot. Ripple ratio cut calculates the overshoot of $S1$ waveforms and rejects pulses with a ratio larger than 0.5%. Negative charge cut compares two methods of total PMT charge calculation. A discrepancy between them also indicates unwanted overshoot in corresponding waveforms. We reject events with such discrepancy larger than 1.5%.

$S1/S2$** band cut**

Noise signals of electrode discharge have pulses similar to $S2$ but with narrower width. Consequently, events with small “S3” are rejected, whose log${}_{10}(S2/S1)$ are smaller than 1.5.

Top/bottom asymmetry cut

The cut is defined as the difference between numbers of PE collected by top and bottom PMT arrays for $S1$ signals. Random coincidence events, usually with large asymmetry, are rejected.

Position quality cuts

This group of cuts reject events with poor position reconstruction quality. Events with unsatisafactory PAF fitting quality are rejected. We also require the difference between positions reconstructed by PAF and TM methods to be less than 40 mm.

Impacts of these successive cuts are shown in Table 3.3 and Fig. 3.3. In Table 3.3 we show the number of survived events and estimated signal efficiencies of the relevant cuts. Detector trigger efficiency at high energy region is estimated to be about 100% and neglected here. With a wide ROI, almost all the NLDBD events would fall in the region and the ROI cut efficiency is neglected. Veto cut efficiency is not needed since we only study single-site NLDBD events with full energy deposition. Fiducial volume cut is reflected in the final 136Xe mass calculation. Efficiencies of signal quality cuts are evaluated conservatively by comparing the number of events before and after applying corresponding cut. For example, the pulse shape cut efficiency is $99.3\pm 0.8$%, averaged over the values of three runs. By combining all the cuts, we calculate NLDBD detection efficiency to be $91.6\pm 0.8$%. Meanwhile, the number of background events has been suppressed by three orders of magnitude.

\tabcaption

Event selection tree and the efficiency for NLDBD signal. Signal efficiency stands for the ratio of NLDBD events that can pass through the step.

Events Signal efficiency (%)

Single scattering cut 38070871 93.4

ROI cut 771356 $\approx$100

Veto cut 433482

Fiducial volume cut 32132

Pulse shape cut 31900 99.3 $\pm$ 0.8

$S1/S2$ pattern cut 31837 99.8 $\pm$ 0.1

Position quality cut 31511 99.0 $\pm$ 0.1

Total

91.6 $\pm$ 0.8

\figcaption

Energy spectrum from 0 to 3 MeV after various cuts. Impacts of different cut can be seen in the plot. Backgrounds are rejected effectively by veto and fiducial volume(FV) cuts. The signal quality cut is tuned specifically for high energy events, and yields a lower efficiency for low energy events far away from the ROI.

4 Monte-Carlo simulation and background estimation

Gamma events from 238U, 232Th, 222Rn, and 60Co originating from detector and xenon medium contribute most to the background in the NLDBD ROI. Spectra from 222Rn and 60Co do not have any obvious peak in the ROI and can not be constrained well with fitting. The background rates are determined with a Q-value-blinded pre-fit with a range between 1 and 3 MeV, excluding data within the $2458\pm 100$ keV signal region. In the extended range, we include the 1173 and 1332 keV peaks from 60Co, as well as additional peak from 40K, 238U, and 232Th.

Background spectral shape in the ROI from various components is estimated with Monte-Carlo simulation based on the Geant4 toolkit [18]. For a given radioactive source originating from a specific detector component, we simulate a large number of events and obtain the energy deposition information in xenon. Detector response is also considered in the simulation. Hits are merged together when their distance in Z direction is smaller than 8.5 mm, which corresponds roughly to the width of a $S2$ defined by the gas gap, below which we could not identify a double scatter event. We also remove the events with multiple unmerged hits or with hits inside veto volume or outside fiducial volume. Background contribution from 222Rn dissolved in the xenon liquid is out of secular equilibrium and considered separately. The 214Bi in the 222Rn decay chain, a beta emitter with the associated 2447 keV gamma ray, is largely removed due to the delayed coincidence of 214Po (see later). Its contribution from the veto region (about 100 kg of liquid xenon) can also be neglected by the active veto of beta energy and the self-shielding effects of liquid xenon. For other isotopes, such as 238U, spectra from all the components are summed up with a weighting factor, which is derived from the relative radioactivity of each component, as reported in Ref. [19].

In addition, we take the two neutrino double beta decay (2NDBD) of 136Xe into consideration. Energy distribution of two emitted electrons from 2NDBD is simulated with the decay0 generator [21], and then the summed electron energy spectrum is added to the fit. The 2NDBD event rate is fixed to 48 per kg$\cdot$yr in the fit region, which is calculated according to 136Xe’s 2NDBD half-life [22] and its abundance in natural xenon.

The fit is done with a software package called Bayesian Analysis Toolkit (BAT) [23]. This software applies Bayes’ Theorem to optimize multiple variables when given a likelihood function. The fit is done for each run independently, and the summed fit result is shown in Fig. 4. The fitted isotope rates are propagated to the final NLDBD fit, used as inputs for Gaussian priors. Four parameters are added in the fit model to describe energy resolution at 1173 (60Co), 1447 (40K), 1764 (214Bi), and 2615 (208Tl) keV position in the spectrum. The energy resolution at any other energy is taken as the interpolation of resolutions of two nearby peaks. Fitted energy resolutions at 2615 keV peak are $4.8\pm 0.1\%$, $4.2\pm 0.1\%$, and $4.3\pm 0.1\%$ respectively. Energy resolution in Run 9 is slightly worse, since the uniformity correction map with the 164 keV gamma peak was obtained with a calibration period only towards the end of that run.

\figcaption

Data and best-fit curves between 1 to 3 MeV to evaluate the background contributions from radioactive sources. The spectrum is a summation of three runs. The shaded region around NLDBD Q-value is not used in the fitting process.

Background rate from radon in xenon medium is cross-checked with a so-called Bi-Po coincidence study. For example, 214Bi, which is after 222Rn in the 238U decay chain, beta decays to polonium nucleus, with a series of signature gamma rays emitted. 214Po quickly decays to 210Pb and a 7.7 MeV alpha particle with a half-life of 164.3 µs. Consequently, decay of 214Bi and 214Po is very likely to happen within a readout window width of 1 ms , which is easy to identify. As discussed in Ref. [20], we select the sequential Bi-Po decay events by finding two adjacent $S1$ signals with time difference consistent with the polonium decay. The concentration of 222Rn is measured to be (6.6 - 8.6) µBq/kg during the three runs, while 220Rn is an order of magnitude lower and ignored in the following analysis. The statistical uncertainty of radon concentration is less than 1% and neglected. The relative systematic uncertainty of 222Rn is about 50%, estimated from the discrepancy of 214Po single event rate and 214Bi-214Po coincidence event rate [11]. We see a reduction of 222Rn concentration from Run 9 to 11, possibly due to changes in the circulation and purification conditions. The results agree with our pre-fit results within $1\sigma$ uncertainties.

5 Double beta decay fit results

The background model of NLDBD ROI differs from the background pre-fit as described in the previous section in a few key aspects. In the ROI fit between 2058 to 2858 keV, energy resolutions are fixed to their values obtained in the pre-fit. We ignore the negligible 40K and 2NDBD contributions in this fit. Event rates from other isotopes are still to be fitted, but assigned with Gaussian priors with their sigma from the pre-fit results.

The fitted NLDBD rate $\Gamma_{0\nu}$ is in unit count per year per nucleus. To add it into model, we convert the decay rate to total count $C$:

[TABLE]

where $T$ is the detector live time, $N_{A}$ is the Avogadro constant, $R=8.9\%$ is the natural isotopic abundance of 136Xe, $M$ is the xenon mass in our fiducial volume, $\epsilon$ is the detection efficiency, and $m$ is the molar mass of xenon.

The best fit model combines data of three runs in one BAT fit. Background rate is independent in each of the runs, while the NLDBD decay rate is forced to be one parameter. In total we have $3\times 4+1=13$ (3 runs, 4 background components) fit parameters for the model. The fitted ROI spectrum is shown in Fig. 5. Fit results of all free parameters are listed in Table 5. The fitted NLDBD rate is $(-0.25\pm 0.21)\times 10^{-23}\,\mathrm{yr}^{-1}$, with an equivalent event rate $-0.93\pm 0.79$ per kg$\cdot$yr

\figcaption

(Bottom) Energy spectrum and best-fit model for the NLDBD ROI. Three runs are fitted simultaneously with shared NLDBD rates and the spectrum is the summation of all the data. The fits are dominated by 2204 and 2615 keV gamma peaks from 238U and 232Th decay chain respectively. An NLDBD rate of $-0.25\pm 0.21\times 10^{-23}\,\mathrm{yr}^{-1}$ is obtained from our best fit model. (Top) The normalized residuals of the fits.

\tabcaption

Summary of best-fit parameters. Event rates are in the unit of counts per kg$\cdot$yr.

Item Run 9 Run 10 Run 11

60Co 14.5 $\pm$ 1.0 13.3 $\pm$ 1.0 11.6 $\pm$ 0.6

222Rn chain 33.4 $\pm$ 6.3 31.9 $\pm$ 5.2 22.9 $\pm$ 3.3

232Th chain 63.6 $\pm$ 2.4 55.8 $\pm$ 2.0 67.6 $\pm$ 1.3

238U chain 36.6 $\pm$ 4.2 34.8 $\pm$ 3.4 41.0 $\pm$ 2.1

NLDBD $-0.93\pm 0.79$

We also estimate systematic uncertainties induced by different parameters, as listed in Table 5. As in Ref. [24], we split the uncertainty contribution into additive part ($\sigma_{\mathrm{add}}$, which is independent of the decay rate) and scaling part ($\sigma_{\mathrm{scaling}}$, which contributes to an uncertainty proportional to the true decay rate). Uncertainties from the signal detection efficiency and fiducial mass only have the scaling part, and are first added to the table.

We carry out a set of pseudo-experiments with toy Monte-Carlo spectrum. Toy spectra are generated according to the fitted background model and mixed with varied NLDBD events. For example, for energy resolution $\Sigma$, we modify its value by its uncertainty in the background model to be $\tilde{\Sigma}$. Hypothetical NLDBD decay rate $\tilde{\Gamma}$ is changed from 0 to $3\times 10^{-23}\,\mathrm{yr}^{-1}$. For each fake NLDBD rate $\tilde{\Gamma}$ we generate 1000 toy Monte-Carlo spectra according to the fitted background model but with $\tilde{\Sigma}$. Then the spectra are re-fitted by imposing the best fit background model with all systematic terms fixed at their best fit values (e.g. $\Sigma$), yielding an NLDBD event rate which is biased from the input. The fitted NLDBD rates are plotted against the input ones and the distribution is fitted with a linear function. Interception of the fit on the Y axis represents the $\sigma_{\mathrm{add}}$ and its slope’s deviation from one is $\sigma_{\mathrm{scaling}}$.

The energy scale at 2458 keV, estimated to be 96+-1% (see Fig. 6), also contributes to the uncertainty of the NLDBD rate. We modify its value by one sigma when producing the pseudo-experiment data, and repeat the constrained fits mentioned above. The uncertainty of the fitting process itself, known as fit bias, is derived similarly but with no change of the systematic sources.

Results of systematic uncertainty evaluation are listed in Table 5. Under the assumption that all the uncertainties are un-correlated, the total systematic uncertainty can be calculated as:

[TABLE]

\tabcaption

Summary of the additive and scaling systematic uncertainties of different contributions. The unit of additive systematic contribution is ($10^{-23}yr^{-1}$). The unit of scaling systematic contribution is (%). Energy resolution and energy scale vary in the three runs, therefore only approximate values are shown.

Fixed value Additive Scaling

Signal detection 91.6 $\pm$ 0.8 (%)

0.83

Fiducial mass 219 $\pm$ 32 (kg)

14.7

Energy scale 96 $\pm$ 1 (%) 0.045

Energy reso. 4.2 $\pm$ 0.1 (%) 0.089 2.47

Fit bias

0.041 0.04

Our final NLDBD rate is $(-0.25\pm 0.21\,(\mathrm{stat.})\pm 0.11\,(\mathrm{sys.}))\times 10^{-23}\,\mathrm{yr}^{-1}$. The result is consistent with zero and we find no evidence of NLDBD in the PandaX-II data. We can calculate a physical Bayesian limit from this best-fit NLDBD rate. The posterior probability distribution and the systematic uncertainties are combined, as shown in Fig. 5. By integrating the positive part of the total likelihood curve to an area of 90% [25], we obtain a limit on NLDBD rate $\Gamma_{0\nu}<0.29\times 10^{-23}$ yr*-1* and NLDBD half-life $T_{1/2}^{0\nu}>2.4\times 10^{23}$ yr at the 90% confidence level (C. L.). The half-life limit would be $2.9\times 10^{23}$ yr if we only consider the statistical uncertainty. Effective Majorana mass can be evaluated under the assumption of light Majorana exchange in the NLDBD process. We adopt the phase space factors from [26] and the nuclear matrix elements calculated from models in [27, 28, 29, 30, 31]. The upper limit of effective Majorana neutrino mass $m_{\beta\beta}$ is found to be in the range of 1.3 and 3.5 eV at the 90% C. L..

\figcaption

NLDBD rate probability distribution. The red scattered dots are from BAT fit result, while the red line is their gaussian fit.

The sensitivity of 136Xe NLDBD search in PandaX-II experiment can also be calculated. We carry out another set of pseudo-experiments, whose spectra are analyzed in the same way as our experimental data. The fitted probability distribution curves are associated with the systematic table and NLDBD half-life limits at the 90% C. L. are calculated. A total of 12000 MC spectra are analyzed and the median value of derived half-life limits is chosen as the sensitivity. The resulting NLDBD half-life sensitivity is $T_{1/2}^{0\nu}>1.9\times 10^{23}$ yr at the 90% C. L..

6 Summary and outlook

We search for possible neutrino-less double beta decay signals from 136Xe with 244 kg$\cdot$yr’s data of PandaX-II detector from March 2016 to August 2018. No evidence of NLDBD event is observed, and the half-life limit for 136Xe is given to be $T_{1/2}^{0\nu}>2.4\times 10^{23}$ yr at the 90% C. L. and the corresponding upper limit for $m_{\beta\beta}$ is the in range of 1.3 to 3.5 eV. This work is the first NLDBD result from dual-phase liquid xenon detector, as far as we are aware. We have developed new techniques of energy reconstruction and event selection, and studied systematically the detector performance in response of MeV-scale events at the PandaX-II detector. The background budget of PandaX-II at high energy region is also presented for the first time in this analysis, and will be used to guide the design and construction of future detectors.

Due to limited performance in energy resolution and background rate, our result is not as competitive as the dedicated NLDBD experiments, such as KamLAND-Zen [5]. However, it shows the multi-physics reach of dual-phase liquid xenon detector besides dark matter WIMP search. PandaX-II has been focusing on signals and background in the keV range while NLDBD region of interest is mostly above 2 MeV, where PMT saturation issue deteriorates resolution significantly. On the other hand, we recognize that the resolution is not intrinsically limited. The XENON1T experiment has shown excellent energy resolution [32], comparable with other liquid xenon TPCs. The PandaX collaboration is building the next generation PandaX-4T [33] detector and will use dual-dynode readout for PMTs, to avoid the PMT saturation issues altogether. PandaX-4T is expected to have a much improved energy resolution for double beta decay physics. PandaX-II also suffers from a high background rate in the NLDBD ROI since we can not take advantage of background suppression techniques developed for WIMP signals in the keV region. Techniques such as nuclear recoil and electron recoil discrimination are specific to WIMP searches. Effective self-shielding with detector medium at low energy become much less so at higher energy. In the proposed 30-ton detector such as DARWIN [34] and PandaX future generations, the fiducial volume in the core of the detector will be much more effectively shielded. At this future stage, dual-phase liquid xenon detector with natural xenon will be able to explore the inverted neutrino mass hierarchy region in the $m_{\beta\beta}$ phase space.