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

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.


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, 76 Ge is widely used in high purity germanium semiconductor detectors [2][3][4], while 136 Xe is usually the target of monolithic liquid detectors [5,6]. GERDA [2], CUORE [7], and KamLAND-Zen [5] experiments have recently established the most stringent half-life lower limit respectively for 76 Ge, 130 Te, and 136 Xe decaying into the ground state of their daughters. Those limits are on the order of -years and show an extremely rare nature of NLDBD, if exists at all.
The PandaX-II detector [8] is a large dual-phase liquid xenon time projection chamber (TPC). The primary scientific goal is to directly detect possible collisions of dark matter candidate particle WIMP (Weakly Interacting Massive Particle) and xenon nuclei [8][9][10]. However, with a target mass of 580 kg of xenon and a natural isotopic abundance of 8.9%, it contains 51.6 kg of 136 Xe 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 identify possible NLDBD signal. We also discuss in details different con-tributions 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) [11], the PandaX-II experiment has been taking physics data since 2016. The detector is a dual-phase xenon TPC [12], as shown in Fig. 1. Detailed description of the detector can be seen in [8]. 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. Veto- ing 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.
In this paper, three sets of physics data from the past 3 years are used for NLDBD analysis, as shown in Table 1. 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 previously published [9,10]. 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. 3. The figure also shows that the drift electron lifetime, which is an indication for detector per-formance, 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.

Energy reconstitution and event selection at MeV region
PandaX-II data analysis at MeV level follows the workflow of dark matter analysis closely [14]. 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 ) and electroluminescence light signal ( ). 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   Chinese Physics C Vol. 43, No. 11 (2019) 113001 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 , or noise using a decision tree based on the characterization of summed waveforms. A typical event may include one signal and one or multiple signals. Each signal paired with the preceding 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 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) [9], 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 and signals. In Run 9 the drift velocity is mm/μs with an electric field of 393.5 V/cm. In Run 10 and Run 11, the drift velocity dropped to mm/μs due to the lowered drift field of 311.7 V/cm. The derived drift velocities are consistent with Ref. [15].

Energy calibration and reconstruction
The detector's response at MeV level is calibrated using external 232 Th source capsules. The source is deployed through the calibration loops inside the TPC outer vessel. A full absorption gamma peak from 208 Tl at 2615 keV can be seen in both and signals. For MeV scale signals, a large 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. 4.
In our analysis, we used pulses from the bottom PMT array to reconstruct the energy of signals to avoid the saturation problem. Saturation is not observed in the bottom PMT array, since signal is usually much smaller than and 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 dependent correction
The number of PE that 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 83 Kr or activated 127m Xe. 127m Xe 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.

S 2
Secondary electrons from xenon ionization are recombined with or attached to the impurities in the drifting process. The electron loss is usually characterized by electron life time in the unit of μs. The measured 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: here W equals to 13.7 eV and is the average work function for a particle to excite one electron-ion pair or one UV photon in xenon [16].
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 with bottom PMT PE number: where stands for ratio between the numbers of PE from the top and bottom PMT array. We expect that photon collection is energy independent and set , For PMT 10900, a typical PMT that experiences saturation problem, clipped peak is also observed as the digitizer saturates.
Chinese Physics C Vol. 43, No. 11 (2019) 113001 as derived in the dark matter analysis [9]. 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 232 Th 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 shown as the red line in Fig. 5 overlaid with the vs. bottom signals.
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. 6. 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.

Event selection cuts
The original cuts from dark matter analysis are inher-S 1/S 2 ited 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, 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 [13,14].

S 2
Single-site cut We focus on single-site (SS) events from NLDBD and reject all events with multiple 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 (MS) 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. keV around the Q-value.
S 1

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 signal. Fiducial volume cut Fiducial volume for NLDBD analysis is optimized according to the distribution of high energy events. Fig. 7 shows the spatial distribution of Bottom PMT signals in 232 Th calibration data. The bright blob on the top-right is the 2615 keV gamma peak. The red line indicates the best anti-correlation between and signals, whose slope is calculated from the ratio of PDE and EEE × SEG. 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. Abnormal "stripes" of excess events are observed vertically, which is due to malfunctioning of PAF position reconstruction algorithm with saturated PMT pulses. This problem also leads to a decrease of event density at the edge of the sensitive volume. To be conservative, we force . 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. 7 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 selfshielding 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 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.

S 1
Pulse shape cuts Two pulse shape cuts are used to reject pulses with positive overshoot. Ripple ratio cut calculates the overshoot of 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 over-shoot in corresponding waveforms. We reject events with such discrepancy larger than 1.5%. band cut Noise signals of electrode discharge have pulses similar to but with narrower width. Consequently, events with small "S2" are rejected, whose log are smaller than 1.5.

S 1
Top/bottom asymmetry cut The cut is defined as the difference between numbers of PE collected by top and bottom PMT arrays for 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 unsatisfactory 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 2 and Fig. 8. In Table 2 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 136 Xe 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 %, averaged over the values of three runs. The uncertainty is derived from the standard deviation of cut efficiencies of three runs. By combining all the cuts, we calculate NLDBD detection efficiency to be %. Meanwhile, the number of background events has been suppressed by three orders of magnitude.

Monte-Carlo simulation and background estimation
(2458 ± 100) ±1 σ Gamma events from 238 U, 232 Th, 222 Rn, and 60 Co originating from detector and xenon medium contribute most to the background in the NLDBD ROI. Spectra from 222 Rn and 60 Co 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. The region keV (about ) is excluded in the pre-fit and NLDBD signal is not considered here. In the extended range, we include the 1173 and 1332 keV peaks from 60 Co, as well as additional peaks from 40 K, 238 U, and 232 Th.

S 2
Background spectral shape in the ROI from various components is estimated with Monte-Carlo simulation based on the Geant4 toolkit [17]. For a given radioactive source originating from a specific detector component, we first generate a large number of events and obtain their energy deposition information in xenon. Detector response is then considered. Hits are merged together when their distance in Z direction is smaller than 8.5 mm, which corresponds roughly to the width of a defined by the gas gap, below which we could not identify a double scatter event.
To compare simulation with calibration data, we simulate a set of 232 Th events originating from the calibration tube. Simulated data go through the same single scatter cut, veto cut and fiducial volume cut as data. The resulting SS/(MS+SS) ratio is 50.2% while the same ratio from calibration data is 55.2%. We assign a conservative ±9% as a systematic uncertainty to account for this difference. To verify the MC spectra, we compare the simulated 232 Th spectrum with calibration data. The mean difference of bin values between the two spectra is calculated to be 6%, which is regarded as the uncertainty of the Monte-Carlo model.
Background contribution from 222 Rn dissolved in the xenon liquid is out of secular equilibrium and considered separately. The 214 Bi in the 222 Rn decay chain, a beta emitter with the associated 2447 keV gamma ray, is largely removed due to the delayed coincidence of 214 Po (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 238 U, 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. [18].

·
In addition, we take the two neutrino double beta decay (2NDBD) of 136 Xe into consideration. Energy distribution of two emitted electrons from 2NDBD is simulated with the decay0 generator [19], and then the summed electron energy spectrum is added to the fit. The 2NDBD event rate is fixed to 48 per kg yr in the fit region, which is calculated according to 136 Xe's 2NDBD half-life [20] and its abundance in natural xenon. Cosmogenically produced 137 Xe [21] contributed about 4 orders of magnitude smaller than 2NDBD and neglected here, thanks mainly to the rock overburden at CJPL.
The fit is done with a software package called Bayesian Analysis Toolkit (BAT) [22]. 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. 9. The fitted isotope rates are propagated to the final NLD-BD fit, used as inputs for Gaussian priors. Four parameters are added in the fit model to describe energy resolution at 1173 ( 60 Co), 1447 ( 40 K), 1764 ( 214 Bi), and 2615 ( 208 Tl) 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 for run 9, 10, and 11. 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.
Background rate from radon in xenon medium is cross-checked with a so-called Bi-Po coincidence study. For example, 214 Bi, which is after 222 Rn in the 238 U decay chain, beta decays to polonium nucleus, with a series of signature gamma rays emitted. 214 Po quickly decays to 210 Pb and a 7.7 MeV alpha particle with a half-life of 164.3 μs. Consequently, decay of 214 Bi and 214 Po is very likely to happen within a readout window width of 1 ms, which is easy to identify. As discussed in Ref.
[23], we select the sequential Bi-Po decay events by finding two adjacent signals with time difference consistent with the polonium decay. The concentration of 222 Rn is measured to be (6.6 -8.6) μBq/kg during the three runs, while 220 Rn 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 222 Rn is about 50%, estimated from the discrepancy of 214 Po single event rate and 214 Bi-214 Po coincidence event rate [10]. We see a reduction of 222 Rn 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 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 40 K 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 is in unit count per year per nucleus. To add it into model, we convert the decay rate to total count C: where T is the detector live time, is the Avogadro constant, R = 8.9% is the natural isotopic abundance of 136 Xe, 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 runs, 4 background components) fit parameters for the model. The fitted ROI spectrum is shown in Fig. 10. Fit results of all free parameters are listed in Table 3. The fitted NLDBD (−0.25 ± 0.21) × 10 −23 yr −1 −0.93 ± 0.79 rate is , with an equivalent event rate per kg·yr. The fitted event rates of 232 Th and 238 U are correlated with 222 Rn and fluctuate in different runs. However, when forcing the Th/U event rates in the three runs to be the same, we observe effectively no change in the fitted NLDBD event rates. We also estimate systematic uncertainties induced by different parameters, as listed in Table 4. As in Ref.
[24], we split the uncertainty contribution into additive part ( , which is independent of the decay rate) and scaling part ( , which contributes to an uncertainty proportional to the true decay rate). Uncertainties from the signal detection efficiency, fiducial mass, and MC rate (SS/(SS+MS) ratio) only have the scaling part, and are firstly added to the table.

Σ ΣΓ
3 × 10 −23 yr −1 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 . 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 and its slope's deviation from one is . The energy scale uncertainties, estimated to be 1% at Q-value, also contribute to the uncertainty of the NLD-BD rate. We modify the energy scale by ±1% when producing the pseudo-experiment data, and repeat the constrained fits mentioned above. For Monte Carlo shape uncertainty, we modify the height of each characteristic peak either by +6% or -6% in the simulated spectrum, and repeat the steps 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 4. Under the assumption that all the uncertainties are un-correlated, the total systematic uncertainty can be calculated as: .
(−0.25 ± 0.21 (stat.)± (sys.)) × 10 −23 yr −1 Γ 0ν < 0.33 × 10 −23 yr −1 T 0ν 1/2 > 2.1 × 10 23 2.9 × 10 23 Our final NLDBD rate is 0.14 . 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. 11. By integrating the positive part of the total likelihood curve to an area of 90% [25], we obtain a limit on NLDBD rate and NLDBD half-life yr at the 90% confidence level (C. L.). The half-life limit would be 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 m ββ 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 is found to be in the range of 1.4 and 3.7 eV at the 90% C. L..

136
T 0ν 1/2 > 1.6 × 10 23 The sensitivity of Xe 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 yr at the 90% C. L..  136 Xe is given to be yr at the 90% C. L. and the corresponding upper limit for is the in range of 1.4 to 3.7 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 dualphase 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,33], comparable with other liquid xenon TPCs. The PandaX collaboration is building the next generation PandaX-4T [34] detector and will use dual-dyn-ode 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 [35] 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 provide an alternative approach besides the proposed KamLAND2-ZEN [36] and nEXO [37].