Searching for Neutrino-less Double Beta Decay of $^{136}$Xe with PandaX-II Liquid Xenon Detector
Kaixiang Ni Yihui Lai Abdusalam Abdukerim Wei Chen Xun Chen Yunhua Chen Xiangyi Cui Yingjie Fan Deqing Fang Changbo Fu Lisheng Geng Karl Giboni Franco Giuliani Linhui Gu Xuyuan Guo Ke Han Changda He Di Huang Yan Huang Yanlin Huang Zhou Huang Peng Ji Xiangdong Ji Yonglin Ju Kun Liang Huaxuan Liu Jianglai Liu Wenbo Ma Yugang Ma Yajun Mao Yue Meng Parinya Namwongsa Jinhua Ning Xuyang Ning Xiangxiang Ren Changsong Shang Lin Si Andi Tan Anqing Wang Hongwei Wang Meng Wang Qiuhong Wang Siguang Wang Xiuli Wang Zhou Wang Mengmeng Wu Shiyong Wu Jingkai Xia Mengjiao Xiao Pengwei Xie Binbin Yan Jijun Yang Yong Yang Chunxu Yu Jumin Yuan Dan Zhang Hongguang Zhang Tao Zhang Li Zhao Qibin Zheng Jifang Zhou Ning Zhou Xiaopeng Zhou

UTF8gbsn









Searching for Neutrino-less Double Beta Decay of ^136Xe with PandaX-II Liquid Xenon Detector
Supported by grants from the Ministry of Science and Technology of China
(No. 2016YFA0400301 and 2016YFA0400302), a Double Top-class grant from Shanghai
Jiao Tong University, grants from National Science Foundation of
China (Nos. 11435008, 11505112, 11525522, 11775142 and 11755001), grants from the Office of Science and Technology, Shanghai Municipal Government (Nos. 11DZ2260700, 16DZ2260200, and 18JC1410200), and the support from the Key Laboratory for Particle Physics, Astrophysics and Cosmology, Ministry of Education.
This work is supported in part by the Chinese Academy of Sciences Center for Excellence in Particle Physics (CCEPP) and Hongwen Foundation in Hong Kong.
Finally, we thank the following organizations for indispensable logistics and other supports: the CJPL administration and the Yalong River Hydropower Development Company Ltd.

    Volkan Dedeoglu and Raja Jurdak
CSIRO Data61, Australia

Email: {name.surname}@csiro.au
Guntur D. Putra, Ali Dorri and Salil S. Kanhere
UNSW Sydney, Australia
    December 30, 2023
====================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================






(PandaX-II Collaboration)



	 

	^1 INPAC and School of Physics and Astronomy, Shanghai Jiao Tong University, 
Shanghai Laboratory for Particle Physics and Cosmology, Shanghai 200240, China

	^2 Yalong River Hydropower Development Company, Ltd., 288 Shuanglin Road, Chengdu 610051, China

  ^3 School of Physics, Nankai University, Tianjin 300071, China

  ^4 Shanghai Institute of Applied Physics, Chinese Academy of Sciences, Shanghai 201800, China

  ^5 School of Physics & International Research Center for Nuclei and Particles in the Cosmos &
Beijing Key Laboratory of Advanced Nuclear Materials and Physics, Beihang University, Beijing 100191, China

  ^6School of Medical Instrument and Food Engineering, University of Shanghai for Science and Technology, Shanghai 200093, China

  ^7Tsung-Dao Lee Institute, Shanghai 200240, China

  ^8School of Mechanical Engineering, Shanghai Jiao Tong University, Shanghai 200240, China 

  ^9School of Physics, Peking University, Beijing 100871,China 

  ^10Department of Physics, University of Maryland, College Park, Maryland 20742, USA 

  ^11School of Physics and Key Laboratory of Particle Physics and Particle Irradiation (MOE), Shandong University, Jinan 250100, China

  ^12Center of High Energy Physics, Peking University, Beijing 100871, China







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 , or equivalently 22.2 of ^136Xe exposure.
At the region around ^136Xe 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 × 10^23 yr at the 90% confidence level, which corresponds to an effective Majorana neutrino mass m_ββ <  (1.3 - 3.5) eV.
This is the first NLDBD result reported from a dual-phase xenon experiment.





Neutrino-less double beta decay, dark matter, underground physics, liquid xenon detector



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



[0]† Corresponding author: ke.han@sjtu.edu.cn
[0]  Spokesperson: jianglai.liu@sjtu.edu.cn





2



§ 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 <cit.>.
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 <cit.>, while ^136Xe is usually the target of monolithic liquid detectors <cit.>. CANDLES-III <cit.>, GERDA <cit.>,  CUORE <cit.>, and KamLAND-Zen <cit.> 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 <cit.> 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 <cit.>.
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 <ref> will describe the PandaX-II detector and the data taking campaigns in year 2016-2018.
In Section <ref>, 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 <ref> and <ref>, 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.




§ PANDAX-II DETECTOR AND DATA TAKING CAMPAIGNS


Located at China Jin-Ping Underground Laboratory (CJPL) <cit.>, the PandaX-II experiment has been taking physics data since 2016.
The detector is a dual-phase xenon TPC <cit.>, as shown in Fig. <ref>.
Detailed description of the detector can be seen in <cit.>.
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. <ref>).
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.



   
    < g r a p h i c s >

		
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 <cit.>.






    
    < g r a p h i c s >

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





In this paper, three sets of physics data from the past 3 years are used for NLDBD analysis, as shown in Table <ref>.
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 <cit.>.
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. <ref>.
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.



        
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

    

    



      
    < g r a p h i c s >

     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 .
    The dash horizontal line indicates the total exposure at 242 , with a corresponding fiducial mass of 219 kg.
    The right Y axis is the electron lifetime while drifting in the field cage.

	







§ ENERGY RECONSTITUTION AND EVENT SELECTION AT MEV REGION


PandaX-II data analysis at MeV level follows the workflow of dark matter analysis closely <cit.>.
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.



 §.§ 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) <cit.>, 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±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±0.01 mm/µs due to the lowered drift field of 311.7 V/cm.
The derived drift velocities are consistent with Ref. <cit.>.



 §.§  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. <ref>.



  
    < g r a p h i c s >

  
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:

    E=W × (S1/PDE+S2/EEE×SEG).

Here W equals to 13.7 eV and  is the average work function for a particle to excite one electron-ion pair in xenon <cit.>.
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:

    S2=S2_Bottom× (T/B+1),

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 <cit.>.
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 × 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 × 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. <ref>.


  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. <ref>.
Peaks from 915 to 2615 keV are used to correct the residual non-linearity, as shown in Fig. <ref>.
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.




    
    < g r a p h i c s >

    	
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 × SEG.






    
    < g r a p h i c s >

		
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.






 §.§ 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 <cit.>.



  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. ±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. <ref> 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 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. <ref> 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.


    
    < g r a p h i c s >

		
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 <ref> and Fig. <ref>.
In Table <ref> 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±0.8%, averaged over the values of three runs.
By combining all the cuts, we calculate NLDBD detection efficiency to be 91.6 ± 0.8%.
Meanwhile, the number of background events has been suppressed by three orders of magnitude.




        
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           ≈100             

    Veto cut                        433482           -             

    Fiducial volume cut             32132            -             

    Pulse shape cut                 31900            99.3 ± 0.8 

    S1/S2 pattern cut             31837            99.8 ± 0.1 

    Position quality cut            31511            99.0 ± 0.1 

    Total                           -                91.6 ± 0.8 

    







    < g r a p h i c s >

    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.	





§ 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±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 <cit.>.
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. <cit.>.

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 <cit.>, and then the summed electron energy spectrum is added to the fit.
The 2NDBD event rate is fixed to 48 per in the fit region, which is calculated according to ^136Xe's 2NDBD half-life <cit.> and its abundance in natural xenon.

The fit is done with a software package called Bayesian Analysis Toolkit (BAT) <cit.>.
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. <ref>.
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±0.1%, 4.2±0.1%, and 4.3±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.




    < g r a p h i c s >

    	
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. <cit.>, 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 <cit.>.
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σ uncertainties.




§ 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 Γ_0 ν is in unit count per year per nucleus.
To add it into model, we convert the decay rate to total count C:

    C = Γ_0 ν× T × N_A × R × M ×ϵ/m ,

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, ϵ 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×4+1=13 (3 runs, 4 background components) fit parameters for the model.
The fitted ROI spectrum is shown in Fig. <ref>.
Fit results of all free parameters are listed in Table <ref>.
The fitted NLDBD rate is (-0.25±0.21) × 10 ^-23 yr^-1, with an equivalent event rate -0.93±0.79 per 




    < g r a p h i c s >

		
(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±0.21 × 10 ^-23 yr^-1 is obtained from our best fit model.
(Top) The normalized residuals of the fits.





        
Summary of best-fit parameters.
    Event rates are in the unit of counts per .
    
    Item                     Run 9               Run 10               Run 11           

    ^60Co                14.5 ± 1.0      13.3 ± 1.0       11.6 ± 0.6   

    ^222Rn chain         33.4 ± 6.3      31.9 ± 5.2       22.9 ± 3.3   

    ^232Th chain         63.6 ± 2.4      55.8 ± 2.0       67.6 ± 1.3   

    ^238U chain          36.6 ± 4.2      34.8 ± 3.4       41.0 ± 2.1   

    NLDBD                    3c-0.93 ± 0.79                  

    




We also estimate systematic uncertainties induced by different parameters, as listed in Table <ref>.
As in Ref. <cit.>, we split the uncertainty contribution into additive part (σ_add, which is independent of the decay rate) and scaling part (σ_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 Σ, we modify its value by its uncertainty in the background model to be Σ̃.
Hypothetical NLDBD decay rate Γ̃ is changed from 0 to 3×10^-23 yr^-1.
For each fake NLDBD rate Γ̃ we generate 1000 toy Monte-Carlo spectra according to the fitted background model but with Σ̃.
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. Σ), 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 σ_add and its slope's deviation from one is σ_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 <ref>.
Under the assumption that all the uncertainties are un-correlated, the total systematic uncertainty can be calculated as:

    σ_ syst (Γ_0 v)=√(∑_i(σ_add, i^2+σ_scaling, i^2Γ^2_0 v)).




        
Summary of the additive and scaling systematic uncertainties of different contributions.
    The unit of additive systematic contribution is (10^-23yr^-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 ± 0.8 (%)                                   0.83   

    Fiducial mass         219 ± 32 (kg)                                     14.7   

    Energy scale          96 ± 1 (%)                      0.045            -   

    Energy reso.          4.2 ± 0.1 (%)                   0.089            2.47   

    Fit bias              -                                    0.041            0.04  

    




Our final NLDBD rate is (-0.25 ± 0.21 (stat.) ± 0.11  (sys.)) × 10^-23 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. <ref>.
By integrating the positive part of the total likelihood curve to an area of 90% <cit.>, we obtain a limit on NLDBD rate Γ_0 ν< 0.29 × 10^-23 yr^-1 and NLDBD half-life T_1/2^0ν> 2.4 × 10^23 yr at the 90% confidence level (C. L.).
The half-life limit would be 2.9 × 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 <cit.> and the nuclear matrix elements calculated from models in <cit.>.
The upper limit of effective Majorana neutrino mass m_ββ is found to be in the range of 1.3 and 3.5 eV at the 90% C. L..




    < g r a p h i c s >

	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ν> 1.9 × 10^23 yr at the 90% C. L..



§ SUMMARY AND OUTLOOK


We search for possible neutrino-less double beta decay signals from ^136Xe with 244 '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ν>2.4 × 10^23 yr at the 90% C. L. and the corresponding upper limit for m_ββ 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 <cit.>.
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 <cit.>, comparable with other liquid xenon TPCs.
The PandaX collaboration is building the next generation PandaX-4T <cit.> 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 <cit.> 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_ββ phase space.





80mm0.1pt



2
34
    fxundefined [1]
     ifx#1
    
    fnum [1]
     #1firstoftwo
     secondoftwo
     
    fx [1]
     #1firstoftwo
     secondoftwo
     
    
    
    
    
    
    noop [0]secondoftwo
    
    ref[1]@startlink#1@href
    href[1]#1@endlink
    anitize@url [0]`
12`$12`&12`#12`1̂2`_12`%12
    startlink[1]
    endlink[0]
    
    rl [1]href #1
    
    
    
    
    
    
    
    
    
    
    
    
    @bib@innerbibempty
    
    [Olive et al.(2014)Olive
      et al.]PDG
      author author K. Olive et al. (collaboration Particle Data Group),      http://stacks.iop.org/1674-1137/38/i=9/a=090001 journal journal Chin. Phys. C volume
      38, pages 090001 (year 2014)NoStop
    [Agostini et al.(2018)Agostini et al.]Agostini:2018tnm
      author author M. Agostini et al. (collaboration GERDA), 10.1103/PhysRevLett.120.132503 journal journal Phys. Rev. Lett. volume 120, pages 132503 (year 2018), http://arxiv.org/abs/1803.11100 arXiv:1803.11100 [nucl-ex] NoStop
    
    [Aalseth et al.(2018)Aalseth
      et al.]Aalseth:2017btx
      author author C. E. Aalseth et al. (collaboration Majorana), 10.1103/PhysRevLett.120.132502 journal journal Phys. Rev. Lett. volume 120, pages 132502 (year 2018), http://arxiv.org/abs/1710.11608 arXiv:1710.11608 [nucl-ex] NoStop
    
    [Wang et al.(2017)Wang et al.]wang2017first
      author author L. Wang et al. (collaboration CDEX), https://link.springer.com/article/10.1007/s11433-017-9038-4 journal journal Sci. China-Phys. Mech. Astron. volume 60, pages 071011 (year
      2017)NoStop
    [Gando et al.(2016)Gando
      et al.]gando2016search
      author author A. Gando et al. (collaboration KamLAND-Zen), https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.117.082503
      journal journal Phys. Rev. Lett. volume 117, pages 082503 (year
      2016)NoStop
    [Albert et al.(2018)Albert
      et al.]Albert:2017owj
      author author J. B. Albert et al. (collaboration EXO), 10.1103/PhysRevLett.120.072701 journal journal
      Phys. Rev. Lett. volume 120, pages
      072701 (year 2018), http://arxiv.org/abs/1707.08707 arXiv:1707.08707 [hep-ex] NoStop
    
    [Umehara et al.(2006)Umehara
      et al.]umehara2006candles
      author author S. Umehara et al. (collaboration CANDLES), in https://iopscience.iop.org/article/10.1088/1742-6596/39/1/093/meta booktitle Journal of Physics: Conference Series,      Vol. volume 39 (organization IOP Publishing,      year 2006) p. pages 356NoStop
    [Alduino et al.(2018)Alduino
      et al.]alduino2018first
      author author C. Alduino et al. (collaboration CUORE), https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.120.132501
      journal journal Phys. Rev. Lett. volume 120, pages 132501 (year
      2018)NoStop
    [Tan et al.(2016a)Tan et al.]Tan:2016diz
      author author A. Tan et al. (collaboration PandaX), 10.1103/PhysRevD.93.122009 journal journal Phys.
      Rev. volume D93, pages 122009
      (year 2016a), http://arxiv.org/abs/1602.06563 arXiv:1602.06563 [hep-ex] NoStop
    
    [Tan et al.(2016b)Tan et al.]tan2016dark
      author author A. Tan et al. (collaboration PandaX), https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.117.121303
      journal journal Phys. Rev. Lett. volume 117, pages 121303 (year
      2016b)NoStop
    [Cui et al.(2017)Cui et al.]cui2017dark
      author author X. Cui et al. (collaboration PandaX), https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.119.181302
      journal journal Phys. Rev. Lett. volume 119, pages 181302 (year
      2017)NoStop
    [Wu et al.(2013)Wu et al.]wu2013measurement
      author author Y.-C. Wu et al., 10.1088/1674-1137/37/8/086001,
      10.1088/1674-1137/37/1/016001 journal journal
      Chin. Phys. volume C37, pages 086001
      (year 2013), http://arxiv.org/abs/1305.0899
      arXiv:1305.0899 [physics.ins-det] NoStop
    
    [Dolgoshein et al.(1970)Dolgoshein, Lebedenko, and Rodionov]DualPhaseTPC
      author author B. A. Dolgoshein, author V. N. Lebedenko,  and author B. U. Rodionov, @noop  journal journal JETP
      Lett. volume 11, pages 351 (year 1970)NoStop
    [Cao et al.(2014)Cao et al.]cao2014pandax
      author author X. Cao et al. (collaboration PandaX), https://link.springer.com/article/10.1007/s11433-014-5521-2 journal journal Sci. China-Phys. Mech. Astron. volume 57, pages 1476 (year
      2014)NoStop
    [Xiao et al.(2015)Xiao et al.]Xiao:2015psa
      author author X. Xiao et al. (collaboration PandaX), 10.1103/PhysRevD.92.052004 journal journal Phys.
      Rev. volume D92, pages 052004
      (year 2015), http://arxiv.org/abs/1505.00771
      arXiv:1505.00771 [hep-ex] NoStop
    
    [Yoshino et al.(1976)Yoshino, Sowada, and Schmidt]yoshino1976effect
      author author K. Yoshino, author U. Sowada,       and author W. Schmidt, https://journals.aps.org/pra/abstract/10.1103/PhysRevA.14.438 journal journal Phys. Rev. volume
      A14, pages 438 (year 1976)NoStop
    [Lenardo et al.(2015)Lenardo, Kazkaz, Manalaysay,
      Mock, Szydagis, and Tripathi]Lenardo2015NEST
      author author B. Lenardo, author K. Kazkaz,
      author A. Manalaysay, author J. Mock, author
      M. Szydagis,  and author
      M. Tripathi, 10.1109/TNS.2015.2481322 journal journal IEEE
      Trans. Nucl. Sci. volume 62, pages
      3387 (year 2015)NoStop
    [Agostinelli et al.(2003)Agostinelli et al.]agostinelli2003geant4
      author author S. Agostinelli et al. (collaboration GEANT4), https://www.sciencedirect.com/science/article/pii/S0168900203013688
      journal journal Nucl. Instrum. Meth. volume A506, pages 250 (year
      2003)NoStop
    [Wang et al.(2016)Wang,
      Chen, Fu, Ji, Liu, Mao, Wang, Wang,
      Xie, and Zhang]Wang2016counting
      author author X. Wang, author X. Chen, author C. Fu, author
      X. Ji, author X. Liu, author Y. Mao, author H. Wang, author S. Wang, author
      P. Xie,  and author
      T. Zhang, 10.1088/1748-0221/11/12/t12002 journal journal
      J. Instrum. volume 11, pages T12002
      (year 2016)NoStop
    [Li et al.(2017)Li et al.]li2017krypton
      author author S. Li
      et al. (collaboration PandaX), https://iopscience.iop.org/article/10.1088/1748-0221/12/02/T02002/pdf
      journal journal J. Instrum. volume 12, pages T02002 (year
      2017)NoStop
    [Ponkratenko et al.(2000)Ponkratenko, Tretyak, and Zdesenko]ponkratenko2000event
      author author O. Ponkratenko, author V. Tretyak,  and author Y. G. Zdesenko, https://link.springer.com/article/10.1134/1.855784
      journal journal Phys. Atom. Nuclei volume 63, pages 1282 (year
      2000)NoStop
    [Albert et al.(2014)Albert
      et al.]Albert:2013gpz
      author author J. B. Albert et al. (collaboration EXO-200), 10.1103/PhysRevC.89.015502 journal journal Phys. Rev. volume C89, pages
      015502 (year 2014), http://arxiv.org/abs/1306.6106
      arXiv:1306.6106 [nucl-ex] NoStop
    
    [Caldwell et al.(2009)Caldwell, Kollár, and Kröninger]caldwell2009bat
      author author A. Caldwell, author D. Kollár,  and author K. Kröninger, https://www.sciencedirect.com/science/article/pii/S0010465509002045
      journal journal Comput. Phys. Commun. volume 180, pages 2197 (year
      2009)NoStop
    [Alduino et al.(2016)Alduino
      et al.]alduino2016analysis
      author author C. Alduino et al. (collaboration CUORE), https://journals.aps.org/prc/abstract/10.1103/PhysRevC.93.045503 journal journal Phys. Rev. volume
      C93, pages 045503 (year 2016)NoStop
    [Casella and Berger(2002)]casella2002statistical
      author author G. Casella and author R. L. Berger, @noop  journal journal
      Pacific Grove, CA  (year 2002)NoStop
    [Kotila and Iachello(2012)]kotila2012phase
      author author J. Kotila and author F. Iachello, https://journals.aps.org/prc/abstract/10.1103/PhysRevC.85.034316 journal journal Phys. Rev. volume
      C85, pages 034316 (year 2012)NoStop
    [Barea et al.(2015)Barea,
      Kotila, and Iachello]barea20150
      author author J. Barea, author J. Kotila,       and author F. Iachello, https://journals.aps.org/prc/abstract/10.1103/PhysRevC.91.034304 journal journal Phys. Rev. volume
      C91, pages 034304 (year 2015)NoStop
    [Rodríguez and Martínez-Pinedo(2010)]rodriguez2010energy
      author author T. R. Rodríguez and author G. Martínez-Pinedo, https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.105.252503
      journal journal Phys. Rev. Lett. volume 105, pages 252503 (year
      2010)NoStop
    [Engel et al.(2014)Engel,
      Šimkovic, and Vogel]engel2014chiral
      author author J. Engel, author F. Šimkovic,  and author P. Vogel, https://journals.aps.org/prc/abstract/10.1103/PhysRevC.89.064308 journal journal Phys. Rev. volume
      C89, pages 064308 (year 2014)NoStop
    [Menendez et al.(2009)Menendez, Poves, Caurier, and      Nowacki]menendez2009disassembling
      author author J. Menendez, author A. Poves,
      author E. Caurier,  and      author F. Nowacki, https://www.sciencedirect.com/science/article/pii/S0375947408008233
      journal journal Nucl. Phys. volume A818, pages 139 (year
      2009)NoStop
    [Mustonen and Engel(2013)]mustonen2013large
      author author M. Mustonen and author J. Engel, https://journals.aps.org/prc/abstract/10.1103/PhysRevC.87.064302 journal journal Phys. Rev. volume
      C87, pages 064302 (year 2013)NoStop
    [Aprile et al.(2019)Aprile
      et al.]XENON:2019dti
      author author E. Aprile et al. (collaboration XENON), 10.1038/s41586-019-1124-4 journal journal Nature volume 568, pages
      532 (year 2019), http://arxiv.org/abs/1904.11002
      arXiv:1904.11002 [nucl-ex] NoStop
    
    [Zhang et al.(2019)Zhang
      et al.]zhang2019dark
      author author H. Zhang et al. (collaboration PandaX), https://link.springer.com/article/10.1007/s11433-018-9259-0 journal journal Sci. China-Phys. Mech. Astron. volume 62, pages 31011 (year
      2019)NoStop
    [Baudis et al.(2014)Baudis,
      Ferella, Kish, Manalaysay,
      Marrodan Undagoitia, and Schumann]Baudis:2013qla
      author author L. Baudis, author A. Ferella,
      author A. Kish, author
      A. Manalaysay, author
      T. Marrodan Undagoitia,  and author M. Schumann, 10.1088/1475-7516/2014/01/044 journal journal
      JCAP volume 1401, pages 044
      (year 2014), http://arxiv.org/abs/1309.7024
      arXiv:1309.7024 [physics.ins-det] NoStop
    
    



