Determination of Double Beta Decay Half-Life of 136Xe with the PandaX-4T Natural Xenon Detector

Precise measurement of two-neutrino double beta decay (DBD) half-life is an important step for the searches of Majorana neutrinos with neutrinoless double beta decay. We report the measurement of DBD half-life of 136Xe using the PandaX-4T dual-phase Time Projection Chamber (TPC) with 3.7-tonne natural xenon and the first 94.9-day physics data release. The background model in the fiducial volume is well constrained in situ by events in the outer active region. With a 136Xe exposure of 15.5 kg-year, we establish the half-life as 2.27 ± 0.03(stat.) ± 0.10(syst.) × 1021 years. This is the first DBD half-life measurement with natural xenon and demonstrates the physics capability of a large-scale liquid xenon TPC in the field of rare event searches.

background model in the fiducial volume is well constrained in situ by events in the outer active region.With a 136 Xe exposure of 15.5 kg-year, we establish the half-life as 2.27 ± 0.03(stat.)± 0.10(syst.)× 10 21 years.This is the first DBD half-life measurement with natural xenon and demonstrates the physics capability of a large-scale liquid xenon TPC in the field of rare event searches.
Two-neutrino double beta decay (DBD) is a second-order weak process with an extremely long half-life [1].Neutrinoless double beta decay (NLDBD), during which no neutrino is emitted, would be direct evidence of the Majorana nature of neutrinos and a clear violation of lepton number conservation [2,3].Therefore, the experimental search for NLDBD is an essential aspect of physics beyond the Standard Model (BSM) [4][5][6].Various experimental projects aim to search for NLDBD in different candidate isotopes.Among them, enriched 136 Xe has been used by EXO-200, KamLAND-ZEN, NEXT, etc [7][8][9][10].PandaX-II has conducted the first search for NLDBD of 136 Xe with a dark matter (DM) detector [11].XENON1T, another natural xenon-based TPC, has demonstrated excellent detector performance at the MeV range [12] and performed a search for NLDBD recently [13].The half-life of 136 Xe DBD has been measured [14][15][16] and the most precise result of 2.165 ± 0.016(stat.)± 0.059(syst.)× 10 21 years is given by EXO-200 [14].
Precise measurement of DBD half-life is an important step for NLDBD searches.For a finite detector resolution, DBD events contribute to the background of the NLDBD region of interest.
Furthermore, a measurement of DBD half-life can help determine the nuclear matrix elements (NME) of the decay process and may shine light on the NME of NLDBD [17][18][19].More recently, there have been growing interests in exploring BSM physics with the electron spectrum of DBD (e.g.[20][21][22][23]).
In this letter, we present a new measurement of 136 Xe DBD half-life with large statistics and a broad energy spectral range with the dual-phase xenon Time Projection Chamber (TPC) of PandaX-4T [24].This represents the first DBD half-life measurement with natural xenon as the target.PandaX-4T has been originally optimized for signals in the low energy keV range for DM searches [25].For the MeV energy range, we developed dedicated energy calibration procedures, reconstruction algorithms, and discrimination of single-and multiple-site events.We systematically select the innermost cleanest part of the detector as our fiducial volume (FV) while utilizing the outer region to robustly characterize the background distributions in situ.
PandaX-4T detector is a large dual-phase xenon TPC with approximately 3.7 tonnes of natural xenon in the active volume [24].The active volume is enclosed by the field cage on the side, and the gate and cathode electrodes on the top and bottom, respectively.The cross-section of the field cage In this analysis, we used the first data release of PandaX-4T, taken from Nov. 28, 2020 to Apr. 16, 2021.The total exposure is 94.9 days and divided into five different data sets, as detailed in Ref. [24].This analysis emphasizes the energy range from 440 to 2800 keV, where the detector response is calibrated with external high-energy gamma sources including 137 Cs, 60 Co, and 232 Th.The gamma sources are encapsulated in small, centimeter-scale stainless steel cylinders, and then inserted into dedicated calibration pipes near the detector.Due to the higher activities, 137 Cs and 60 Co are placed at 397 cm and 375 cm away from the central axis of the TPC respectively, while 232 Th with lower activity is placed at 84 cm away.All of them are approximately at the same height as the geometrical center of the TPC active volume.
For gamma events in the MeV range, energy loss via multiple Compton scattering dominates and can be identified via the timing profile of S2 signals.In the summed S2 waveforms of PMT arrays, multiple site (MS) events register different pulses as shown in Fig. 1 Top.The S2 pulses of MS events can overlap in time due to the diffusion effect and misidentified as one single site (SS) event.
To discriminate SS and MS events, the summed pulse of all PMTs for a S2 is firstly smoothed by a locally weighted regression algorithm [27].We then identify peak(s) on the smoothed pulse, and 0 0.2 0.4 0.6 0.  a SS event is defined where all peaks are identified in a time window of ±0.75 µs, corresponding to ∼2.2 mm in the vertical direction.Otherwise, the events are labeled as MS events.
PMTs may suffer from the saturation effect with high energy S2 signals.For events in the MeV energy range, a typical PMT may collect up to thousands of photoelectrons (PEs).We observe obvious waveform distortion when the total charge collected is above 900 PEs or so, as illustrated in Fig. 1 Bottom.PMT saturation effect is corrected by matching the rising slope of the saturated PMT waveforms to the non-saturated waveforms in the same event.We use the summed PMT waveforms with charges in the range of 50 to 900 PEs as templates of non-saturated waveforms.The waveform template is then scaled to match the rising slope of the saturated waveforms in a dynamic reference interval by minimizing the reduced χ 2 between the two waveforms.The scaled template offers an estimate of the true charge of saturated waveforms [12].We also validate the effectiveness of the desaturation with a bench measurement using PMTs illuminated by high-intensity photons with S2-like timing profiles.The desaturation protocol is applied to the top PMTs for all events with distorted waveforms.For events in the energy range of 1 to 3 MeV, the average correction factor is ∼3.0.The S2 waveforms for the bottom PMTs are not corrected for saturation, but a residual correction is applied to the reconstructed energy (see later).On the other hand, we do not observe apparent saturation for S1 signals in the energy region of interest.
Position reconstruction is significantly improved with desaturated PMT charges.At high energy, we use the photon acceptance function (PAF) method [28] to reconstruct the position.PAF describes the expected charge detected by a certain PMT as a function of the event's position and is determined iteratively from data.A likelihood function is constructed to describe the charge distributions among each PMTs for an event, and the horizontal position is determined by maximizing the likelihood.
In Energy calibration coefficients of SS events are inherited from the low energy analysis and applied Data quality selection cuts to remove unphysical events and select electron recoil events are developed based on calibration data, including that from 137 Cs, 60 Co, 232 Th, 83m Kr, and 220 Rn.
Non-electron recoil and alpha events can be removed with a cut on the charge of S1 and/or S2 as well as the ratio between the two.The relative size of S1 charges collected by the top and bottom PMT arrays and drift time can be cross-compared to reject accidental coincidence events and events originating from the gate electrode.All cuts were optimized and determined with the last 9.6 days of data after going through the high-energy specific data processing, while the rest of the data is blinded.The total cut efficiency is calculated to be (99.4± 0.4)%, while the uncertainty is estimated by the difference in different detector volumes.These cuts are applied to all physics data, and the efficiency is validated to be the same as that from 9.6 days of data.
Detailed simulation with Geant4-based Monte Carlo framework BambooMC [31] [24], the tail of its beta spectrum has a marginal impact on our result and has not been included in the fit.
DBD signal spectrum is also simulated with BambooMC.The energy of two electrons from DBD is generated with the Decay0 package [35] as input for our simulation.For DBD events with energy greater than 440 keV, the SS fraction is 97.4% with a fractional uncertainty of 1.7% according to our detector response simulation.
A cylindrical FV with a radius of 33.2 cm and a height of 66.3 cm in the geometrically center part of the detector is selected for the final fit, with the range in the Z direction pre-determined by the event rate distribution in the ROI.The FV is then determined from 220 Rn and 83m Kr calibration data.Both internal calibration sources are expected to be evenly distributed in the active volume.
The FV is defined where the proportionality between the event counts and geometrically calculated We constructed a simultaneous fit with the binned likelihood function defined as where N obs ij , N ij are the observed and expected event numbers in the j-th energy bin of the i-th region.N ij are modeled as Eq. 2 according to PDFs of the DBD spectrum S Xe ij and the background components B k ij given as (category, isotope) pairs and 214 Pb listed in Table 1.The PDFs are weighted by the number of counts n Xe and n k for signal and backgrounds respectively.For background, n k is fixed while the nuisance parameters η k denotes the fractional difference between the prior and fitted counts.In the likelihood function, a Gaussian penalty is added for each background except 214 Pb based on the prior fractional uncertainty σ k derived from Table 1.The prior 222 Rn level is measured in situ with the number of high-energy alpha particles emitted in the decay chain with small uncertainty.However, the number of 214 Pb may be smaller due to progenies attaching to electrodes and the inner surface of the detector [36].Therefore, the number of 214 Pb events is unconstrained in the fit.
The fitted spectra in Region 1 and the combined results of the other three regions are shown in Fig. 4. In the FV, 136 Xe DBD events are dominant.The total number of DBD events is 17468 ± 257.
A parallel fit is also performed with the RooFit package [37] giving consistent results.
The fitted results of all background sources are listed in Table 1 in comparison with the expected counts.The most noticeable background is from 214 Pb, the progeny of 222 Rn.The ratio of 59.6% between the fitted count of 214 Pb and the expected count of 222 Rn, which is estimated from its alpha events, represents the aforementioned depletion of 214 Pb.The large sensitive volume of PandaX-4T helps determine the external radioactive contaminations more accurately and robustly, as demonstrated by the smaller uncertainties on the fitted number of counts.Agreements within two sigmas between our fit results and radioactive assays are observed for most of the contributions.For comparison, the best fit results are used to calculate the expected spectrum for the regions outside of the FV, as outlined by the three dashed rectangles in Fig. 2. Within our DBD ROI, the largest difference between expected and measured rates is 2.3% and the agreement is within 1% when the three dashed regions are considered together.
We performed the fit by varying bin size from 1 keV to 40 keV and lower (upper) fit range from 440 (2600) to 600 (3000) keV.The impact of both changes is at the 1% level or smaller.Systematic uncertainties may also come from the mismatch between simulated and measured spectra in the  FV.Three independent effects are studied including uncertainties in energy resolution, energy scale, and the relative variation of weights among the four regions.To cleanly evaluate the systematic effect of energy reconstruction, the maximum deviation of both the reconstructed energy peak and energy resolution is included.PDFs from the simulation are generated with energy resolutions scaled from 0.8 to 1.2 times the measured resolutions.We also shifted the peak positions of PDFs by up to ±10 keV.The relative weight n k in Eq. 2 is calculated from the geometry or scaled by the measured 222 Rn rates in each region.Different impact on DBD half-life between the two sets of weights is treated as the regional weight systematics, as listed in Table 2.All the background PDFs are generated assuming secular equilibrium of the U/Th chains.If the early and late parts of the decay chains are allowed to float independently in our spectrum fit, we observe a 2.0% shift in our DBD result, mainly due to the thorium chain.The shift is treated as a systematic uncertainty from  In summary, we have measured the DBD half-life of 136 Xe with the 94.9 days of PandaX-4T data.Our analysis uses an FV at the very center of the detector with 649.7 ± 6.5 kg of xenon and a total number of 17468 ± 257 DBD events was observed in the energy range of 440 to 2800 keV.The final result 2.27 ± 0.03(stat.)± 0.10(syst.)× 10 21 years is one of the most precise measurements to date and the first result using a natural xenon dark matter detector with robust background control.Our DBD analysis is a major step forward for the NLDBD searches with a large-scale multiple-purpose liquid xenon TPC.PandaX-4T has continued data taking with the physics goal of dark matter direct detection, double beta decay searches, and measurement of neutrinos from astrophysical sources.We are also developing upgrade plans to enhance the detector response at the MeV level and possibly lower the background rate by replacing detector components.More and higher quality data will further improve our measurement of DBD half-life and understanding of the background relevant to NLDBD searches, and may be used to explore potential new physics [20][21][22][23].
is a polygon with 24 sides and its inscribed circle has a diameter of 1.185 m.The liquid xenon (LXe) inside the field cage is 1.185 m tall.The top and bottom PMT arrays have 169 and 199 three-inch Hamamatsu PMTs installed.All PMTs are read out with custom-designed bases for high voltage biasing and signal readout [26].PMT signals are recorded by commercial 14-bit, 250 MHz sampling rate CAEN digitizers.An energetic event in the active volume of the detector generates a prompt scintillation signal (S1) and a delayed electroluminescence signal (S2).Energy deposition in the active volume generates ionization electrons and detectable prompt scintillation light signals S1.The ionization electrons drift up to the top of the active volume where they are extracted into the gaseous region.The electric field in the gaseous region is much stronger, causing electrons to undergo an electroluminescence process to generate S2.The energy and position of each event can be calculated from the amplitude and temporal information of S1 and S2.The vertical (denoted as Z) position of an event can be determined from the drift velocity of electrons and the time delay between S1 and S2.The relative amount of photons seen by each PMT on the top PMT arrays is used to reconstruct the positions in the horizontal X-Y plane.

Figure 1 :
Figure 1: Top: Waveform of a sample MS event with one S1 peak and three S2 peaks identified by the MS algorithm.The separation is indicated by the blue dots; Bottom: An example of saturated S2 waveform (blue) with a charge of approximately 3.8×10 3 PEs.The desaturated S2 (red) from a waveform template is scaled by matching the rising slope in the reference interval (green), resulting in a corrected charge of ∼ 1.2 × 10 4 PEs.

Fig. 2
Top, the distributions of positions reconstructed from desaturated and uncorrected PMT charges, respectively, are shown as a function of R 2 , where R is the radial position.The wiggles in the uncorrected curve are significantly smoothed out by the desaturation process.The detector response to physical signals in different parts of the detector is non-uniform and the correction is done in situ with data.S2 signals suffer from the absorption of electrons on electronegative impurities in LXe while drifting up.The electron lifetime is calculated with the 164 keV γ-ray peak from 131m Xe run by run and used to correct S2 signals with an exponential function.Electron loss may be convoluted with a saturation effect for large S2 signals, especially for the top PMT array.Therefore, the electron lifetime correction is extracted based on S2 charges collected by the bottom PMT array (called S2 b ), and the corresponding correction is applied to S2 b in the physics data.The average electron lifetime by dataset ranges from 800.4 to 1288.2 µs, while the maximum drift time is from 800 to 841 µs [24].Gaseous 83m Kr calibration source is injected into the detector via the circulation loop [29] to map out the non-uniform response of the detector.To get a response map to S1 signals, 41.5 keV events emitted from 83m Kr are grouped into different voxels in the detector and the mean S1 charge in each voxel is calculated to represent the detector response at the corresponding positions.For the S2 signals, the response map is generated in the X-Y plane.The S1 and S2 signals for any event are corrected according to their position to achieve a uniform response throughout the active volume.

Figure 2 :
Figure 2: Top: Event position distribution in radius-squared (R 2 ) of the physics data in the energy range of [440, 2800] keV reconstructed from the PAF algorithm with PMT charges before (green) and after (red) desaturation.Bottom: Event distributions of the physics data in the energy range of [440, 2800] keV in Z vs R 2 coordinates, with the color bar showing counts in each bin.The FV is divided into 4 regions, projected as 4 rectangles with numbers accordingly.The outer regions outlined by the dashed rectangles are used for cross-validation of the results obtained from the FV.The top and bottom shaded regions close to the PMT arrays are excluded in the top panel.

Figure 4 :
Figure 4: Results of the simultaneous fit.For better visualization, Region 1 is shown individually (top) and the other three regions are shown together (bottom).Fitted spectral of background components of (category, isotope) pairs are grouped and plotted based on isotopes.In each figure, data and fit results are shown in the top panel and residuals in the bottom panel.The χ 2 /NDF (number of degrees of freedom) as an indicator of the goodness of fit is shown for each figure.
years.This result agrees with the previous measurements from KamLAND-ZEN, EXO-200, and NEXT collaborations[14][15][16], as shown in Fig.5.Our measurement with PandaX-4T, which is initially designed as a dark matter direct detection experiment, demonstrates the physics reach of large liquid xenon TPC on multiple fronts.Thanks to its emphasis on low energy, our analysis threshold is at least 260 keV lower compared to the previous measurements.With a large active volume and self-shielding, we can aggressively select the innermost and cleanest part of the detector for the measurement.The outer part of the active volume is used to cross-check the background and signal model in situ.Such advantages can also be exploited in other physics searches in the MeV range in the future.

Figure 5 :
Figure 5: Comparison with half-life measurements from enriched xenon experiments [14-16].The box and cap-tipped error bars represent statistical and systematic uncertainties respectively.
[30,32]blished with full detector geometry, including top and bottom PMT arrays with individual PMT and base, field cage, stainless steel vessels, supporting frame, and water shielding tanks.LXe responses to ER in the MeV range are modeled from standard NEST 2.0 construction[30,32], with light yield, charge yield, and recombination parameters extracted from calibration.Energy depositions simulated in Geant4 are converted into individual S2s, which are then smeared in time with drift-time-correlated Gaussian diffusion measured in the data.The pseudo-S2 waveforms are then piped through the SS/MS discrimination algorithm.The resulting SS and MS spectra are compared to the corresponding spectra in the calibration runs.No systematic deviation is observed as a function of energy, as top, bottom, and side.The top category includes the top flanges of the vessels and the top PMT assembly, which consists of the PMT array, readout circuits, cabling, and the mechanical supporting structure.The counterpart bottom PMT assembly and the bottom hemisphere of the vessels are grouped as the bottom category.The side category is composed of the field cage and cylindrical barrel of the vessels.Other detector components are found to have negligible background contributions and thus not included.The weighted sums of expected background counts in the ROI from four radioactive contaminations are listed in Table1.222Rn emanated from the inner surface of the detector and circulation pipes is the major internal contamination.214 Pb and 214 Bi, progenies of the 222 Rn, contribute mostly to the ROI of DBD.97% of the beta decays from 214 Bi can be rejected together with their subsequent alpha decay of 214 Po with a half-life of ∼ 163 µs [34].The remaining 214 Bi activity is less than 0.1 µBq/kg, which makes a negligible contribution to our ROI.Therefore, only the contribution from 214 Pb is considered and simulated with BambooMC.85 Kr beta decay also contributes to the lower end of the ROI with end-point energy of 687 keV.However, with an extremely low concentration level at 85 Kr/Xe ratio of 6.6 ± 4.2 × 10 −24 [33]within 440 keV to 2800 keV, the overall agreement is at 1.7% level (Fig.3).When taking the bin-by-bin deviation into account, the fit result (see later) shifts by 1.5%.Conservatively, we use 1.7% as a systematic uncertainty in the final DBD half-life calculation.The majority of the background events in our region of interest (ROI) from 440 to 2800 keV are from radioactive contamination of detector components and impurities in the xenon.Radioactive contamination of 60 Co, 40 K, and the U/Th chain from major detector components is assayed with dedicated setups such as high purity germanium detectors, with results summarized in Ref.[33].We generate expected background contributions in SS spectra according to simulation, assay results, and the discrimination algorithm.The major detector components are grouped into three categories, 232Figure 3: Top: SS (magenta) and MS+SS (cyan) spectra of 232 Th calibration data.Bottom: Comparison of SS fraction between MC (shaded green) and data (black), with an average difference of 1.7%.The uncertainties of the MC spectrum are shown in the green shaded band.denoted

Table 1 :
Expected and fitted contribution of background contaminations originating from the top, bottom, and side of the detector and LXe inside.All values are reported in the number of counts in the FV.5%, and the uncertainty of the FV cut is estimated as 1.0%.The FV is further divided into 4 regions, as shown in Fig.2Bottom and data spectra are reconstructed in each region.Region 1 is the innermost and cleanest region.Region 2, 3, and 4 are on outside of Region 1, where the external radioactive contaminations from the top, bottom, and side of the detector, respectively, have more impact.

Table 2 :
[39]ary of relative systematic uncertainties.equilibriumdecaychain.Moreover, a recent re-evaluation of 214 Pb beta decay spectrum pointed out that the spectrum calculated by Geant4 may not be precise[38]In our ROI, the new theoretical calculation changes the spectrum by approximately -6% to 30% from low energy to high energy, if we assume the same correction can be applied to 214 Pb decay to excited states as well as the ground state.Therefore, a sizable impact of 2.0% is observed in our fit result.Similarly, 136 Xe DBD shape may also influence the half-life measurement.An alternative shape assuming a single-state dominance hypothesis (SSD)[39]is used in the spectrum fit and the result differs from the baseline choice by 0.36%.The total relative systematic uncertainty is 4.5%, including all the systematics listed in Table2summed in quadrature.With a live time of 94.9 days, a selection and SS efficiency of 99.4% and 97.4%, respectively, and a DBD event fraction of 86.3% within our fit range, we obtain a final DBD half-life measurement of 2.27 ± 0.03(stat.)± 0.10(syst.)× 10 21 Aside from what has already been mentioned earlier, we have also evaluated the systematics related to the 136 Xe abundance.The fluctuation of gas phase pressure 2.064 ± 0.008 bar during the data taking leads to an LXe density of 2.8502±0.0036g/cm 3 , based on the saturated vapor properties of Xe.Therefore, the total xenon mass in the FV is 649.7 ± 6.5 kg, where the uncertainty comes from the FV cut and LXe density.We adopt the reported 136 Xe abundance of 8.86% in natural xenon and measure the xenon samples from the PandaX-4T TPC by a residual gas analyzer (RGA)[40]for cross-check and error estimation.With and without the relative ionization efficiencies of different xenon isotopes [41] taken into account, the measured isotopic abundance of 136 Xe is (9.03 ± 0.03)% and (8.79 ± 0.03)%, respectively.Conservatively, we use (8.86 ± 0.17)% in the half-life evaluation, resulting in a 136 Xe mass of 59.6 ± 1.3 kg in the FV.