Double-Weak Decays of $^{124}$Xe and $^{136}$Xe in the XENON1T and XENONnT Experiments

We present results on the search for double-electron capture ($2\nu\text{ECEC}$) of $^{124}$Xe and neutrinoless double-$\beta$ decay ($0\nu\beta\beta$) of $^{136}$Xe in XENON1T. We consider captures from the K- up to the N-shell in the $2\nu\text{ECEC}$ signal model and measure a total half-life of $T_{1/2}^{2\nu\text{ECEC}}=(1.1\pm0.2_\text{stat}\pm0.1_\text{sys})\times 10^{22}\;\text{yr}$ with a $0.87\;\text{kg}\times\text{yr}$ isotope exposure. The statistical significance of the signal is $7.0\,\sigma$. We use XENON1T data with $36.16\;\text{kg}\times\text{yr}$ of $^{136}$Xe exposure to search for $0\nu\beta\beta$. We find no evidence of a signal and set a lower limit on the half-life of $T_{1/2}^{0\nu\beta\beta}>1.2 \times 10^{24}\;\text{yr}\; \text{at}\; 90\,\%\;\text{CL}$. This is the best result from a dark matter detector without an enriched target to date. We also report projections on the sensitivity of XENONnT to $0\nu\beta\beta$. Assuming a $275\;\text{kg}\times\text{yr}$ $^{136}$Xe exposure, the expected sensitivity is $T_{1/2}^{0\nu\beta\beta}>2.1 \times 10^{25}\;\text{yr}\; \text{at}\; 90\,\%\;\text{CL}$, corresponding to an effective Majorana mass range of $\langle m_{\beta\beta} \rangle<(0.19 - 0.59)\;\text{eV/c}^2$.


I. INTRODUCTION
The XENON collaboration acquired science data with the XENON1T experiment at the INFN Laboratori Nazionali del Gran Sasso (LNGS) in Italy from November 2016 until December 2018. Its primary goal was the search for interactions between xenon nuclei and dark matter (DM) in the form of weakly interacting massive particles (WIMPs) [1,2]. In addition to these nuclear recoils, the detector was also sensitive to other rare processes that could be measured as energy depositions on atomic electrons in xenon, electronic recoils (ER). In particular, the collaboration reported the first direct observation of the two-neutrino double-electron capture (2νECEC) in 124 Xe with 4.4 σ significance [3]. The low background rate of the experiment and its good energy reconstruction and resolution up to the MeV region [4] also allow for the search for the neutrinoless doubleβ decay (0νββ) of 136 Xe. This potential will be extended with XENONnT, the latest experiment within the XENON program, owing to its approximately three times larger active xenon mass and six times smaller background rates [5].
The yet unobserved 0νββ is a nuclear transition predicted by extensions of the Standard Model (SM). Two neutrino double-β decay (2νββ) is allowed in the SM and has been observed in 136 Xe making it a candidate isotope to search for a 0νββ peak at the Q-value of Q ββ = (2457.83 ± 0.37) keV [6,7]. The currently best lower limit on the 136 Xe 0νββ half-life, T 0ν 1/2 is set by > 2.3 × 10 26 yr [8] at 90% confidence level (CL).
A detection of 0νββ or neutrinoless double-electron capture (0νECEC) would demonstrate the violation of total lepton number and prove the existence of a non-zero Majorana component of neutrino mass. Under the assumption of light Majorana neutrino exchange, the halflife is related to the effective Majorana mass, m ββ , by [9] m ββ 2 = m 2 e G 0ν |M 0ν | 2 T 0ν Here, G 0ν is the phase-space factor in units of yr −1 [10], M 0ν is the dimensionless nuclear matrix element (NME), and m e is the electron mass in eV/c 2 . Since m ββ can contain phase cancellations from the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix, it is sensitive to the neutrino mass hierarchy [11,12]. While the phase-space factor can be calculated with relative precision, theoretical uncertainties are associated with the choice of the NME. The central values of the most extreme 136 Xe NMEs presented in [13] range from M 0ν = 1.550-4.773 [14,15] and are considered when interpreting 0νββ decay limits in this work. This illustrates that M 0ν is a major source of uncertainty on m ββ . Although there is no direct correspondence between the neutrinoless and two-neutrino NMEs, the measured half-lives of two-neutrino decays such as 2νECEC can be used as a benchmark for different NME calculation approaches [13].
In this work, we perform 2νECEC and 0νββ peak searches in the measured ER energy spectrum of XENON1T and assess XENONnT's sensitivity to 0νββ of 136 Xe using simulated data. The paper is organized as follows. Sec. II A and Sec. II B give an overview of the XENON1T detector and the XENONnT detector, respectively. Sec. II C highlights the background components relevant for the 2νECEC and 0νββ decay searches and their constraints for background simulations. Sec. II D details the fitting method which is employed to derive results. Sec. III summarizes an updated search for 2νECEC in 124 Xe following an extension of the signal model to include captures from higher electron orbitals and using a larger exposure compared to the previous analysis [3]. Sec. IV reports on a 136 Xe 0νββ decay search in XENON1T. Sensitivity projections for the XENONnT experiment are discussed in Sec. V. A summary of the results and an outlook on 0νββ search with WIMP detectors are given in Sec. VI.

II. GENERAL ASPECTS OF THE ANALYSES
XENON1T [1] and XENONnT are designed as dualphase xenon time projection chambers (TPC). These cylindrical detectors are filled with liquid xenon (LXe) and have a thin xenon gas layer at the top. Several electrodes in the liquid and gas enable the application of electric fields. Photomultiplier tubes (PMTs) at the top and bottom read out the signals from particle interactions. A particle interaction in LXe produces excitation, ionization and heat. The total number of measurable quanta depends on the energy deposition of the incident particle and the interaction type, e.g., nuclear recoil or electronic recoil. Excitation occurs in the form of excited xenon dimers that decay to the ground state by emitting scintillation light at 175 nm [16]. The electron-ion pairs from the ionization process can recombine leading to further light emission. The resulting primary scintillation signal is registered by the PMTs and denoted S1. Full recombination is suppressed by an electric drift field that moves the electrons away from the interaction site and towards the liquid-gas interface. There they are accelerated into the gas gap by an extraction field and produce a secondary scintillation signal, S2, proportional to the number of extracted electrons. In these analyses, we consider events from particle interactions that have a summed PMT waveform containing at least one S1 and S2 pair. The time scale for such events is given by the maximum drift time of O(1 ms) in both detectors.
The 3D-position of an interaction is reconstructed using the distribution of the S2 light in the top PMT array (x, y) and the time delay between the S1 and S2 signals (z). Using the self-shielding of xenon and this position information allows for the definition of a fiducial volume (FV) with reduced background levels from external sources, i.e., located outside the detector and the detector materials themselves. In XENON1T data, S1 and S2 signal sizes were corrected accounting for the position-dependent measurement efficiencies. Moreover, the measured interaction positions were corrected for inhomogeneities of the drift field [17].
The total deposited energy of an event is characterized by the weighted sum of corrected S1 and corrected S2. The energy calibration in XENON1T was performed using monoenergetic peaks. The anti-correlation between the S1 and S2 signals at different energies allowed us to compute the photon detection efficiency g 1 and the charge amplification factor g 2 [4]. These parameters were used as weighting factors for the two observables S1 and S2 in the energy calibration, which we call combined energy scale (CES). The analyses in the following sections were carried out in the CES parameter space.
A. The XENON1T experiment The XENON1T [1] TPC had a height of 97 cm and a diameter of 96 cm. Two arrays of 127 and 121 Hamamatsu R11410-21 3-inch PMTs were arranged above and below the sensitive volume of the TPC, respectively. The active volume consisted of 2 t of LXe out of a total of 3.2 t in the detector. The TPC side walls were made of Polytetrafluoroethylene (PTFE) reflective panels to enhance the light collection efficiency. Two electrodes, a cathode placed at the bottom of the TPC and a gate ∼2.5 mm below the liquid-gas interface, produced a drift field of 81 V/cm. An anode placed 2.5 mm above the liquid-gas interface created an extraction field of 8.1 kV/cm. The cryostat was immersed at the center of a stainless-steel tank, filled with 700 t of ultra-pure demineralized water, used to shield environmental gammas and neutrons. The tank was instrumented with 84 PMTs to actively tag muons and muon-induced backgrounds through the detection of Cherenkov light.
The data used for the analyses presented here were acquired between February 2017 and September 2018 during the main science run of the experiment (SR1) and a second run (SR2) targeting research and development. Subsets of SR1 and SR2 were selected for the specific analyses and are described in more detail in Sec. III and IV.

B. The XENONnT experiment
XENONnT is the successor experiment to XENON1T. It was commissioned in the second half of 2020 and started operations shortly thereafter. It reuses several subsystems already developed for XENON1T, with additional radon removal, LXe purification, neutron veto and xenon gas storage systems. The TPC has an active region of 133 cm in diameter and 148 cm in height containing 5.9 t of LXe. The cryostat holds 8.4 t of LXe in total. Two hexagonal arrays contain 253 and 241 PMTs at the top and bottom, respectively. In order to avoid digitizer saturation from large S2s at MeV energies, the PMTs in the top array are read out with an amplification factor of ×0.5, in parallel to the ×10 amplification used for DM searches. This secondary readout was specifically installed for 0νββ searches. Muons are suppressed by means of the same tagging system developed for XENON1T. Additionally, a novel neutron veto (NV) system uses 120 PMTs inside an optically-separated volume around the cryostat to detect signals originating from the capture of radiogenic neutrons [5]. For the projections in section V, we assume the same detector operating conditions as in [5].

C. Electronic Recoil Background
The analyses presented in Sec. III, IV and V require modeling the individual background sources via Monte Carlo (MC) simulations. The simulated backgrounds include radioactive impurities in the detector components and the xenon target itself, as well as solar neutrinos [3,18]. The background composition in each analysis depends on the energy range and the chosen FV. In the following, we describe simulation aspects and backgrounds that are common to the 2νECEC and 0νββ analyses. More specific background contributions are discussed in the respective sections.
The background energy spectra were obtained with the XENON1T [19] and XENONnT [5] MC simulation frameworks, respectively. First, energy depositions from radioactive decays were simulated using the implementation of the detector geometry in Geant4 [20][21][22]. Next, individual energy depositions from the same Geant4 event were clustered based on their relative S2 sizes and their zseparation. This clustering mimics the finite resolution in the reconstruction of multiple nearby energy depositions in the detector. Events with a single energy cluster are denoted as single-site (SS) events while those with multiple clusters are multi-site (MS). Measured waveforms of SS and MS events would contain a single S1 signal but, while SS events have only one S2, MS events may contain multiple resolved S2s from energy depositions at different depths. The S2-based MS event resolution deteriorates for deeper events due to longitudinal diffusion of the electron cloud. Clustering distances for XENON1T range from 6.5 mm at the top to 11.5 mm at the bottom of the TPC. These were determined using simulated waveforms. Events with multiple S2s at the same depth can be identified using PMT hit pattern information. This information was not included in the simulations and accordingly not used for MC event selections.
After the clustering, SS events were selected from the MC simulation within the FV chosen for each analysis. The resulting energy spectra were convolved with a Gaussian representing the measured energy resolutions of the SR1 and SR2 data [4,18]. Monoenergetic peaks from target-intrinsic sources without significant Compton-scattering or Bremsstrahlung contributions were modeled as single Gaussian lines with the standard deviation given by the energy resolution.
The material backgrounds considered in XENON1T and XENONnT originate mainly from the 238 U and 232 Th decay chains as well as from 40 K and 60 Co. Their contributions were constrained using the radioassay results of the XENON1T [23] and XENONnT [24] detector materials. Due to its 5.27 yr half-life, the decay of 60 Co after the radioassay, is taken into account. We assume that production of 60 Co by activation underground is negligible. In the MC simulation of the uranium and thorium decay chains, we take possible decay chain disequilibrium into account. The early and late parts of the uranium chain were split at 226 Ra, and the 232 Th decay chain was split at 228 Th. For the XENON1T analysis, the full 232 Th chain was simulated. In order to account for disequilibrium, the partial chain starting at 228 Th was also simulated. In fits of simulated background spectra to measured data it could be added or subtracted from the full chain in the fit, depending on the observed disequilibrium. This different treatment with respect to the 238 U chain was caused by the internal processing of the decay chains and ensured that all expected γ-rays were present in the simulations. For XENONnT, both parts of the 232 Th chain were treated independently.
After applying the FV selection, no further spatial information was included in the background model. The energy distributions from the same isotope but originating from different materials are essentially identical. Therefore, the relative contributions from different materials to the background from each isotope were fixed in the analysis using the screening measurements. With this, a single scaling parameter for each isotope was required in each background model.
The 222 Rn emanation from materials of the detector to the LXe target induces an intrinsic background contribution. The two most relevant radon daughter isotopes identified for the analyses presented in this work are 214 Pb and 214 Bi. The latter predominantly undergoes β-decay to 214 Po [25]. Subsequently, 214 Po decays via α-emission with a half-life of 164 µs. The close timing coincidence of the two decays, with respect to the event time scale of O(1 ms), allows for their tagging and enables effective rejection of 214 BiPo events inside the active volume. A rejection efficiency of 99.8 % is assumed for XENONnT. For XENON1T, the 214 Bi background from radon emanation is discussed in section IV B.
The 2νββ decay of 136 Xe features a continuous energy spectrum with the endpoint at Q ββ . Theoretically calculated distributions of the energies and relative emission angles of the two electrons [10,26] were used as input to Geant4. The IBM-2 higher-state dominance (HSD) [10] model of the 2νββ process was used. The resulting energy spectrum was normalized according to the expected decay rate corresponding to a half-life of T 2νββ 1/2 = (2.165 ± 0.016 stat ± 0.056 sys ) × 10 21 yr [27] and the measured isotopic abundance η136 Xe = (8.49 ± 0.04 stat ± 0.13 sys ) × 10 −2 mol/mol in XENON1T. Considering these uncertainties, the background contribution over the whole energy range can be constrained with a relative uncertainty of 3 %. The difference between the HSD spectrum and an alternative single-state dominance (SSD) spectrum [28] was not considered as a source of systematic uncertainty due to the subdominant contribution of this background in the analyses. Due to a combination of sub-percent energy resolution at Q ββ [4] and low decay rate, the 2νββ contribution in the region of interest (ROI) is expected to be several orders of magnitude lower than the material background for both XENON1T and XENONnT.

D. Fit method and limit setting
The results presented below were derived using the standard procedure of a Poisson binned log-likelihood where nuisance parameters are profiled [29]. The likelihood reads N i is the measured number of events in each energy bin, λ i is the number of expected background events as a function of the nuisance parameters θ. The number of expected signal events in bin i is denoted with n s i and depends on the signal strength µ s and the nuisance parameters θ. Each constrained nuisance parameter θ j has an expected mean value µ j with a standard deviation σ j . Details regarding the set of nuisance parameters for each analysis are given in the specific sections. The binned likelihood for signal and background model fits to measured data in Secs. III and IV are normalized such that it can be interpreted as in a χ 2 fit [30], e.g., in terms of goodness of fit. Following the nomenclature in [30], this goodness of fit measure is labeled χ 2 λ . The test statistic for the 0νββ searches in Secs. IV and V is where quantities with a single hat denote the set of parameters which correspond to the unconditional maximum of the likelihood while quantities with two hats denote the set of parameters maximizing the conditional likelihood. Under certain conditions the test statistic q(µ s ) follows an asymptotic distribution which is given by a χ 2 distribution with one degree of freedom [29]. The distribution of q(µ s ) is estimated by toy-MC simulations to validate the assumption of asymptoticity. We report only the upper edge of the Feldman-Cousins confidence interval if an excess is smaller than 3 σ. Similar to [31] this imposes overcoverage for very small signals. The 3 σ significance threshold only serves as the transition point between reporting one-and two-sided intervals, and was decided prior to the analysis to ensure correct coverage. All derived limits correspond to a 90 % CL. The parameter of interest is the event rate A ββ of the double-weak processes 2νECEC or 0νββ. For a measured A ββ , the half-life is where M A = 0.131 kg /mol is the xenon molar mass, N A is Avogadro's constant, and η Xe is the isotopic abundance of the xenon isotope, 124 Xe or 136 Xe, under consideration. The SS efficiency SS is the fraction of signal events that are identified as SS events. With an electron mean free path in LXe smaller than 3 mm, the majority of the two electrons emitted in the 0νββ decay are detected as SS events. For reference, the spatial (z) resolution of XENON1T is ∼ 8 mm for S2 signals larger than 10 3 photo-electrons (PE) [17]. An efficiency loss occurs when Bremsstrahlung is emitted by one of the double-β electrons leading to an MS event. The SS efficiency for 0νββ was estimated with MC to be 90.3 % in XENON1T and 91.0 % in XENONnT. For the X-rays and Auger electrons emitted in the double-electron capture, the SS efficiency is 100 % due to the sub-millimeter mean free path of the quanta.

III. EXTENDED SEARCH FOR 124 XE TWO-NEUTRINO DOUBLE-ELECTRON CAPTURE IN XENON1T
The 2νECEC analysis presented here builds on the previous result [3] with an increased exposure and additionally considers double-electron capture contributions from higher atomic shells. As in single-electron capture, double-electron capture rates chiefly depend on the overlap between the electron and nuclear wave functions [35]. Since the electrons in the s-orbitals of the K-and Lshell (L1) feature the largest overlap, only they have been commonly considered in theoretical studies and in the interpretation of experimental data [5,18,32,36].  [3]. These were obtained using the Dirac solutions of the bound electron wave functions for a finite size nucleus in table V of [32] as in [18,33]. The improved double-capture fraction calculation in this work considers all shells up to the N5-shell with the squared amplitudes of the radial wave functions in [34]. The individual subshells are added together in the table rows. The energy ranges for the captures are from the orbital with the lowest to that with the highest binding energy. The captures with energies below the analysis energy threshold of 10  and the updated model considering shells from K to N5 (orange). The peak widths reflect the energy resolution of XENON1T. Uncertainties associated with the peak scaling, the peak positions and the peak widths are not shown for visibility. The energy region that was excluded due to the low-energy excess reported in [18] is marked in blue.
However, with xenon's 54 atomic electrons, the M-, Nand O-shells with s-, p-, d-and higher orbitals should also contribute to the total double-electron capture decay rate. Values for the squared amplitudes of the radial wave functions up to the N5-shell are tabulated in [34]. The corresponding relative capture fractions are given in Tab. I together with the results when considering Kand L1-captures only. The respective signal models are illustrated in Fig. 1. Taking the additional shells into account slightly decreased the relative fraction of KK-, KLand LL-captures with respect to all decays. No literature values for the O 1,2,3 wave functions were available. Their values were approximated by scaling the N 1,2,3 values with a factor 1/4 -approximately the scaling between the tabulated M-and N-shell wave functions. With this, O-captures were estimated to present O(0.1) % corrections to the other capture fractions. Accordingly, we do not include O-captures in the signal model. We note that this approach is still simplified compared to the calculation approach for single-electron capture outlined in [35], where the Q-value of the decay as well as the energies, parities, and angular momenta of the nuclear states involved in the decay are considered. However, more work is needed to extend this treatment for double-electron capture. The simplified treatment was considered as a systematic uncertainty in the analysis.

A. Reconstruction and cuts
The data used for this analysis comprises 226.9 livedays from SR1, subdivided into two partitions with nom- inal (SR1 a ) and increased (SR1 b ) background produced by neutron calibrations as in [18]. Moreover, 24.3 d from SR2 were added to the dataset. The total exposure of the combined dataset is 0.93 t × yr of which 0.68 t × yr in SR1 a overlap with the data used in [3]. With a measured abundance of η124 Xe = (9.94 ± 0.14 stat ± 0.15 sys ) × 10 −4 mol/mol the 124 Xe isotope exposure is 0.87 kg × yr. The energy calibration parameters from [4] were used for the SR1 data after an upgrade of the XENON1T data processor. They were derived anew for SR2 following the same method as laid out in [4]. Event positions were reconstructed with a neural network. As in [3], the SR1 a dataset was analyzed in a 1.5 t superellipsoid FV which was subdivided into an inner 1.0 t cylinder and an outer 0.5 t shell. The SR1 b and SR2 datasets were analyzed in cylindrical 1 t FVs [1] only. For SR2 no outer FV was used, since the respective cuts and energy calibration had been defined for the 1 t cylinder. For SR1 b the addition of an outer FV was not expected to yield a significant increase in sensitivity due to the larger background from neutron calibrations. The reconstruction-induced systematic uncertainties on the masses for all FVs were estimated using homogeneously distributed 83m Kr calibration data and are below 1 %. An overview of the four datasets with their respective live times and FV masses is given in Tab. II.
The data quality criteria (cuts) from [18] were applied to the SR1 data. Two cuts using PMT hit pattern information were adapted for SR2, taking into account changes in the PMT configuration between SR1 and SR2. A cut evaluating the difference in reconstructed event positions from two different algorithms was discarded since not all position reconstruction algorithms had been updated for SR2. As the cut targets singular outlier events, its omission has no impact. The energy ranges employed for the analysis are 10-200 keV for the central 1 t FV and 10-160 keV for the outer 0.5 t volume. The lower fit bounds were chosen in order to exclude the low-energy excess observed in [18]. Above 160 keV and 200 keV in the inner and outer FVs, respectively, Compton scat-tering lead to MS events. We were therefore unable to obtain a clean calibration sample to determine the acceptance of our SS cuts. Thus, we set these energies as our upper fit bounds.
In [3], misidentified events from the intrinsic calibration isotope 83m Kr led to the presence of a secondary lower-energy peak from this isotope (cf. Sec. III C). This obscured possible KL-, KM-and KN-capture peaks. A new cut from [18] reduces these events to the 10 −4 level. This enabled the addition of the double-electron capture peaks from higher shells to the signal model.

B. 125 I background
A key background in this search comes from 125 I, a daughter isotope of 125 Xe which is produced by neutron capture on 124 Xe in the detector as well as in the detectorexternal xenon purification loop [3,33]. The half-life of 125 Xe is 16.9 h [38] and its electron capture decay with subsequent γ-emission was observed outside of the energy region used for the 2νECEC search. The half-life of 125 I is 59.4 d and its three atomic deexcitation cascades from K-, L-and M-shell electron captures are merged with the 35.5 keV γ-ray from the daughter 125 Te [38,39] resulting in a peak energetically close to the 2νECEC signal peaks. The relative fractions for K-, L-, and M-shell electron captures as well as the merged peak energies are given in Tab. III.
The absolute 125 I background contribution in SR1 a and SR1 b was constrained by integrating the activity model from [3] over the new data selection. The integration yielded N 125 I SR1a = (10 ± 5) t −1 and N 125 I SR1 b = (100 ± 20) t −1 . The SR2 model was derived with the same method as in [3], which tracks the time evolution of the 125 Xe parent activity. The number of 125 I nuclei N125 I in the TPC is connected to the number of 125 Xe nuclei N125 Xe by the differential equation The decay constants λ125 Xe and λ125 I describe the production and decay of 125 I, respectively. The purification time constant τ pur accounts for continuous iodine removal in the xenon purification loop and is a fit parameter in the model. Combining decay and purification leads to an effective time constant τ eff . The replacement of the purification system's pumps with a new ultra-clean magnetically-coupled piston pump [40] during SR2 allowed increasing the purification flow from (50 ± 2) slpm to (79 ± 2) slpm. As the iodine removal was expected to be proportional to the purification flow, the ratio of purification time constants  [39]. The peak positions are obtained by adding the 35.5 keV γ-energy to the energies of the atomic relaxations [43,44]. Uncertainties on the X-ray transition energies are at the eV-level, far below the XENON1T energy resolution and not listed.

Decay
Capture With this, the 125 I model for SR2 was divided into a period before and after the pump installation. The pump replacement increased the 85 Kr background level in all subsequent data, most likely due to introduction of airborne krypton to the system during the operation. However, the pump replacement also reduced the 222 Rn background level due to its lower radon emanation [40]. A further reduction was achieved in the final SR2 data when the krypton distillation column was operated in a specialized radon distillation mode [41,42], but no significant change in the total background rate was observed due to the elevated 85 Kr background. Consequently, both of these periods were modeled by one background rate parameter. The respective background levels are discussed in the next section. The fit of the 125 I data is shown in Fig. 2

C. Additional background sources
Apart from 125 I, the target-intrinsic isotope 214 Pb from continuous 222 Rn emanation was the dominant background source in this analysis. For its direct β-decay to the 214 Bi ground-state we assume a branching ratio of (11.0 ± 1.0) % [25] in our simulations. An approximate activity concentration of 214 Pb was inferred from the αdecay rates of other radon daughters in the detector. For SR1 a and SR1 b , the corresponding 222 Rn activity concentration is (13.3 ± 0.5) µBq/kg [45]. Due to online radon distillation and the lower emanation of the new pump, the time-averaged 222 Rn activity concentration in SR2 was (10.1 ± 0.3) µBq/kg, (76 ± 4) % of the SR1 mean, leading to a reduction in 214 Pb background [45]. Since there is no direct measurement of the reduction, the scaling parameter for 214 Pb was left unconstrained in the fit. More details can be found in Sec IV B.
Anthropogenic 85 Kr was removed from the XENON1T target by cryogenic distillation before SR1 [41]. The natural krypton concentration in xenon was monitored over time by taking regular samples that were measured with rare gas mass spectrometry (RGMS) [46]. Considering that natural krypton contains 85 Kr at the 2×10 −11 level [31,41], the RGMS measurements were used to constrain the associated background. The mean concentra-tions in SR1 and SR2 were nat Kr/Xe = (0.7 ± 0.1) ppt and nat Kr/Xe = (1.2 ± 0.2) ppt in mass, respectively. The latter arises from nat Kr/Xe = (0.7 ± 0.2) ppt before the pump installation and nat Kr/Xe = (2.0 ± 0.4) ppt thereafter.
Several intrinsic backgrounds exhibited an explicit time-dependence. Neutron activation during calibrations led to backgrounds from the aforementioned 125 I and 125 Xe peaks, the metastable 131m Xe peak at 163.9 keV with a half-life of 11.84 d, and the merged γ+β spectrum of 133 Xe with a half-life of 5.25 d. Due to the proximity to neutron calibrations, the activation level in SR1 b was increased compared to SR1 a . The more frequent neutron calibrations in SR2 also led to a larger activation background contribution. Apart from 125 I, no constraints were placed on the rates of these backgrounds.
A 83m Kr peak at 41.5 keV was present in all datasets due to trace amounts of the parent isotope 83 Rb in the xenon recirculation system; its activity decreased with the 83 Rb half-life of 86.2 d. The decay of 83m Kr is a two step isomeric transition [34] with energy depositions of 31.2 keV and 9.4 keV from conversion electrons. The 157 ns half-life of the 9.4 keV state leads to the merging of the corresponding S1s and S2s in XENON1T which produces the 41.5 keV peak. If both energy depositions can be distinguished, the events are removed by cuts. In a fraction of 83m Kr events, the S1 from the 9.4 keV transition was wrongly classified as an S2 which was not included in the energy reconstruction leading to a secondary lower-energy 83m Kr peak. As outlined in Sec. III A, a dedicated cut was developed in order to address this reconstruction artifact.
The elastic scattering of solar neutrinos off atomic electrons constituted a subdominant background compared to the intrinsic and material background components and was implemented as in [19]. The backgrounds from detector materials and 136 Xe were implemented as outlined in Sec. II C.

D. Fit method and parameters
Following the methodology of Sec. II D, a binned likelihood was constructed for each of the four measured energy spectra in the SR1 a inner and outer FVs as well as the SR1 b and SR2 cylinders, and used for a simultaneous fit of the signal and background models to the measured data. The SR1 in a and SR1 out a datasets were fitted with a 1 keV binning while the SR1 in b and SR2 data were fitted in 2 keV bins due to the lower exposure. Different binnings of 0.5 keV, 1.0 keV, 1.5 keV and 2.0 keV were tested for all datasets, but did not significantly affect the results. Relative differences in the obtained 2νECEC half-lives were less than 3 % and small compared to the systematic uncertainties. The definition of the total likelihood for the signal and background model included 51 fit parameters. Of these, 29 parameters were constrained. Tables with constraints and best-fit values for all parameters can be found in the appendix. Twelve parameters were shared among all datasets. They included the 2νECEC signal and the background sources that were constant in time, such as detector construction materials, solar neutrinos and 136 Xe. The properties of the residual 83m Kr reconstruction artifact, its position (µ83m Kr,misID ), width (σ83m Kr,misID ) and relative frequency with respect to the main 41.5 keV peak (f83m Kr,misID ), were also shared. Its constraints were derived from 83m Kr calibration data.
The remaining 39 parameters were not shared among all datasets and can be grouped into four different categories. The first category includes homogeneously distributed backgrounds from 85 Kr and 214 Pb which could be averaged in time over SR1 and SR2.
The second category contains backgrounds that could not be averaged over the entire science run, but were sufficiently long-lived to distribute uniformly within the detector. The decay rates of the neutron-activated peaks of 125 I and 131m Xe, as well as those for 125 Xe and 133 Xe, were dependent on their temporal proximity to the neutron calibrations. Due to its 16.9 h half-life, 125 Xe was only present in SR1 b .
The third category contains the parameters describing the acceptance. These were constant over time in a given science runs, but differ in the inner and outer volumes of the detector. The SR1 acceptances for the inner and outer volumes were parameterized by linear functions of reconstructed energy that were fitted to 220 Rn calibration data in order to derive the parameter constraints. A constant parameterization was used for SR2.
The fourth and last category consists of parameters that were both time-and position-dependent, so they were fitted individually for each dataset. For the 83m Kr decay rate A83m Kr , the time dependence originated from the decay of a 83 Rb contamination with a half-life of 86.2 d. The spatial dependence is a feature of the event reconstruction. Since 83m Kr undergoes a two-step decay, events only appeared as a 41.5 keV peak if the S1s and S2s from the subsequent isomeric transitions were merged by the data processor. While the S2s were always identified as a single signal in the bulk of the detector due to their O(µs) width, the S1s could sometimes be distinguished depending on their separation in time and the decay position in the TPC. For the same 83m Kr activity in the inner and outer volumes of the TPC, this led to different 41.5 keV peak areas.
The remaining parameters implemented systematic uncertainties on the energy resolution and reconstruction in the fit. The energy resolutions σ E of the monoenergetic peaks with true energy µ E were parameterized as [4] with the constrained fit parameters a res and b res for each dataset. The SR1 parameters were constrained using the parametrization from [4]. For the SR2 resolution, fits of monoenergetic calibration lines were used. The simulated spectra were smeared prior to the fit with the same function.
In order to account for a possible bias in the energy reconstruction, the fitted energy E fit for each signal and background component could be shifted from the simulated energy E by adding a linear energy-dependent shift in the fit Here, a shift and b shift refer to fit parameters modeling this energy shift. The shifts were applied independently for all datasets in order to account for possible different behaviors in the FVs and for possible temporal drifts. This is obvious for SR2 where the energy reconstruction and resolution parameters were different from those determined for SR1. For SR1 these parameters were averaged over the entire science run and only determined in the inner detector volume. Since no independent calibration data is available over the whole duration of SR1 that could be used to formulate constraints, a shift and b shift were left unconstrained for all datasets.

E. Fit results
The best-fit combined signal and background model is shown in Fig. 3. The reduced χ 2 λ /n dof of the fit is 517/508 = 1.02. This considers 530 data points, 51 fit parameters and 29 pulls from constrained parameters with the sum of squared pulls being Σ = 7.6. The smallness of the pull contributions is attributed to the fact that a pull either does not originate from a statistical confidence interval or that the respective parameter is more strongly constrained by the auxiliarly constraint than by the science data. This is the case for the parameters describing the remaining part of the 83m Kr reconstruction artifact, as the total area of the peak is more than two orders of magnitude below the signal rate. Therefore, the smallness of Σ does not indicate a problem with the fit. The reduced χ 2 λ /n dof just from residuals and excluding pulls is 509/479 = 1.06.
The spectra in SR1 in a and SR1 out a are featureless except for the monoenergetic peaks from the 2νECEC signal, 83m Kr and neutron-activation. The spectrum in the outer volume has a larger slope due to the increased material background contribution. As expected, SR1 in b and SR2 exhibit larger neutron-activated peaks as well as the step from the merged β + γ signature of 133 Xe. The contribution from neutron activation in SR2 is lower than in SR1 b , due to the selection of datasets at least 50 d away from neutron calibrations, but higher than in SR1 a . The smallest 83m Kr peak is found in SR2 due to the largely decayed 83 Rb contamination.
The best-fit activity concentration of 214 Pb in SR1 is close to the XENON1T target value of 10 µBq/kg with 214 Pb SR1 = (9.3 ± 0.4) µBq/kg.
The significance of the 2νECEC signal was derived from the χ 2 λ profile of A 2νECEC as shown in Fig. 4. The   best-fit double-electron capture rate is The difference ∆χ 2 λ = 49.4 between the best-fit rate and a null result yields a significance of 7.0 σ for the presence of a double-electron capture signal. This marks the first significant detection (> 5 σ) of a two-neutrino doubleelectron capture in any isotope. Moreover, it is the first measurement of this process that leverages the signatures of higher-shell KL-, KM-, KN-and LL-captures.
The resulting 2νECEC half-life using Eq. (4) is The systematic uncertainty has four individual contributions given in Tab. IV: the cut acceptance, the exposure, the 124 Xe isotopic abundance, and the theoretical uncertainties on the relative fraction of KK-, KL and highershell captures. The first three contributions were calculated with Gaussian uncertainty propagation in Eq. (4) while the fourth was derived by comparing fit results using two different signal models. The uncertainty on the acceptance was obtained by calculating the average acceptance over the entire energy range for each dataset. The difference of the exposure-weighted sum of these average values from unity was taken as the systematic uncertainty. The uncertainty on the exposure was obtained from the FV uncertainties of each dataset that were multiplied with the corresponding live times and added in quadrature. The total uncertainty on the 124 Xe isotopic abundance was obtained by adding the statistical and systematic uncertainties in quadrature.
As stated in Tab. I, approximate values for the doubleelectron capture fractions from different shells were used. These were calculated from the overlap of the nuclear and electronic wave functions. In order to determine the systematic uncertainty arising from this approximation, the full analysis was also carried out with a simplified model including only captures from the K-and L1-shells. As the scaling of the 2νECEC model is predominantly determined by the KK-peak, the increased capture fraction in the simplified model lead to a half-life that is 6.3 % longer. The absolute difference is used as the systematic uncertainty.

F. Comparison with theory and other experiments
The new result can be compared with the previously measured 2νKK-half-life T 2νKK 1/2 = (1.8±0.5 stat ±0.1 sys )× 10 22 yr from [3]. We use the KK-capture fraction of 72.4 % to compute T 2νKK 1/2 = (1.5 ± 0.3 stat ± 0.1 sys ) × 10 22 yr for this work. The datasets partially overlap, so the statistical uncertainties are correlated. However, the analyses used different data processor versions, cuts and energy reconstructions. Consistency checks of both results were carried out using the 0.68 t × yr data contained in both analyses. It was found that the small difference between both results for T 2νKK 1/2 can be accounted for by the updated signal model, the improved energy reconstruction and the larger cut acceptance in this work, together with the independent systematic uncertainties as well as the 33 % larger exposure. XENON1T analysis [3], the agreement with the QRPA (2013) calculation [36] is improved. The value range from QRPA (2015) [49] is consistent with our new result at the 2 σ level. Both the ET and the NSM calculations are compatible with our new result [50]. While the central value of the first XENON1T result was less than 1 σ below the 90 % CL lower limit of XMASS, the new result is approximately 2 σ below the XMASS limit. Future xenon-based detectors with lower backgrounds and larger exposures will further probe 2νECEC to improve experimental constraints on NME calculations for proton-rich nuclides. The best-fit rate from this work would result in XENONnT detecting approximately 6000 double-electron capture events in its projected 20 t × yr total exposure. With a reduction in background by a factor of ∼ 6 [5], the half-life could be measured with a precision at the few-percent level and the relative capture fractions could be investigated. In this regard theoretical input on the relative capture fractions as well as the double-hole energies is needed. Moreover, with more exposure and less background the 2νECEC can be used as an ideal internal energy calibration source, and the remaining two-neutrino and hypothetical neutrinoless decays of 124 Xe [51] could become accessible.

IV. SEARCH FOR 136 XE NEUTRINOLESS DOUBLE-β DECAY IN XENON1T
In contrast to dedicated 0νββ experiments with xenon inventories enriched in 136 Xe [8,52], the isotopic composition of the XENON1T target was close to that of natural xenon with an abundance of 136 Xe as mentioned in Sec. II C. With a tonne-scale fiducial mass and two years of measurement time, an isotope exposure of 36.16 kg×yr was achieved, approaching exposures of dedicated experiments with enriched targets. The data used in this analysis are a subset of the SR1 dataset introduced in Sec. II A. Data periods when the neutron generator was in the water tank close to the cryostat were removed from the data selection to avoid an elevated high-energy γ-ray background level from thorium and radium decay chain isotopes in the neutron generator's materials. The total live time of the dataset is 202.7 days. The analysis was performed in the energy range between 1600 keV and 3200 keV in order to include multiple γ-peaks that helped to constrain material background components and covered the endpoint region of the 214 Bi β-spectrum. A blinding cut between 2300 keV and 2600 keV was applied to the dataset.

A. Event selection
The event selection criteria for signal-like interactions were developed on the blinded science data as calibration sources with energies close to Q ββ were not available in XENON1T. The applied cuts are based on those from [17] and were adapted to higher energies. Firstly, data quality criteria were applied to remove events in coincidence with muon veto triggers and data acquisition busy periods. Periods with light emission in the PMTs [53], causing abnormal data rates, were also removed. The 0νββ signal is expected to be an SS interaction, while events involving Compton scattering will typically be MS interactions. Thus, SS events were selected by rejecting events with a second S2 whose size, width, and PMT hit-pattern were compatible with the S1, such that the secondary S2 and the primary S1 would form a valid event. A multi-S1 cut, based on the size of the second largest S1 in an event, rejected interactions with multiple S1s originating from pile-up. One source of pile-up is 214 BiPo decay, discussed in Sec. II C, with two subsequent decays occurring in the same event.
Two different position reconstruction algorithms, a neural network using TensorFlow [54] and an algorithm using a fit of the S2 hit pattern on the top PMT array [17], were required to give compatible results. This removed events close to the edge of the TPC or in regions where non-functioning PMTs were located. We also required the reconstructed position to be compatible with the observed S2 hit-pattern of the top PMT array. Cuts based on the fraction of light detected by the top and bottom PMT arrays, for both the S1 and S2 signals, were ef-fective at removing events from energy depositions in the gaseous xenon layer. ERs were identified by the ratio of the S1 to the S2 signal size with 98 % efficiency. Finally, we applied a cut requiring the S2 width to be compatible with the expected diffusion of drifted electrons from the reconstructed depth (z) in the TPC.
The individual cut acceptances were determined with three different techniques. The exposure loss from data quality criteria was factored into the live time. Cuts whose acceptance was tested on controlled samples of data have fixed acceptances corresponding to the fraction of the parameter distribution passing the cut. The rest of the cut acceptances were determined iteratively by comparing the number of remaining events after a set of cuts with the number of events after applying the same cut set, except for the one under investigation. The combined cut acceptance was then determined by multiplication of the individual cut acceptances per energy bin. The result was interpolated with a quadratic spline weighted by acceptance uncertainties at each data point [55]. This provided a continuous acceptance parametrization over the full energy range of interest. Since it was not possible to differentiate between removed signal and background events, the iterative method provided only a lower limit on the signal acceptance. In the blinded region this was approximately flat and extrapolated to be > 88 % at Q ββ . An upper limit of 97.5 % was determined by considering only cuts with a fixed acceptance as outlined above. In the later fit of the signal and background models, discussed in Sec. IV C, the acceptance was allowed to float between the lower and upper limits.
An inner FV was selected based on a sensitivity figure of merit S vol in order to maximize the signal to noise ratio: where m is the target mass which scales linearly with the number of signal events and B is the number of expected background counts. Two control regions close to Q ββ in the science data were used for this study. They were defined as ±4 σ E intervals around the 214 Bi and 208 Tl peaks at 2204.1 keV [25] and 2614.5 keV [56], respectively, excluding the data in the blinded region. The resolution at these energies was σ E /E = 0.8%. The TPC's active volume was binned in a 9 × 9 grid in squared radius r 2 and depth z containing equal masses. The sensitivity figure of merit S vol was computed in each bin using the sum of events from both control regions. The resulting grid of sensitivity values was then smoothed to 100 contour levels of the same S vol , which was fitted with two semisuperellipsoid functions. The maximum allowed depth of the FV was −94 cm in order to avoid TPC regions with possible field distortion close to the cathode region. Finally, S vol was computed for the volume enclosed in the fitted contours. This resulted in an optimal FV containing (741 ± 9) kg, shown in Fig. 6  towards the bottom of the detector due to the presence of more material at the top of the TPC and less shielding from the xenon in gaseous phase.

B. Background model
The background model for the 0νββ search accounts for backgrounds from intrinsic and external sources. The background model was validated by a fit to the science dataset between 1600 keV and 3200 keV, excluding the blinded region.
The dominant background in the ROI, defined as the 2 σ region around Q ββ , was due to γ-rays from trace amounts of radioactive isotopes in detector components, as already discussed in Sec. II C. The main contributors were the late parts of the primordial 238 U and 232 Th decay chains as well as 60 Co. Most notable were the full absorption peak of the 2447.9 keV γ-ray of 214 Bi [25] and Compton scatters from the 2614.5 keV 208 Tl line [56]. Additionally, a peak at 2505.7 keV originating from two 60 Co γ-rays at 1173.2 keV and 1332.5 keV [57], detected in coincidence as an SS event, was expected. The early parts of the primordial decay chains did not contribute to the background in the ROI, but appeared in the 1600-3200 keV fit range.
The 2νββ decay of 136 Xe and the β-decay of 214 Bi, entering the FV by radon emanation, were considered as intrinsic background sources. Details regarding their contributions are discussed in Sec. II C. Since the 214 BiPo tagging efficiency was not known from external measurements, the intrinsic 214 Bi background component was not constrained in the fit. The measured spectrum is continuous up to the endpoint at (3270 ± 11) keV. Spectral features occur where the emitted β and subsequent γrays are merged into a single energy deposition. Decays occurring in the LXe shell outside of the TPC could not be tagged since the α from the 214 Po decay did not enter the active region of the detector. Only the γ-rays following the β-decay could be registered inside the active volume. Thus, they had to be treated separately from the TPC contribution. The background fit constraint of (10 ± 5) µBq/kg for the 214 Bi activity concentration was informed by α-decay measurements [45]. Neutron-activated 137 Xe and ERs induced by 8 B solar neutrinos are negligible compared to the material background in XENON1T and were not considered here.

C. Fit to the blinded data and sensitivity
The methodology for fitting and limit setting was introduced in Sec. II D. The set of nuisance parameters comprises scaling factors for all simulated backgrounds. Additionally, we considered a combined cut acceptance parameter , which was allowed to move between the lower bound of 88% and the upper bound of 97.5% in the ROI within 1 σ of its Gaussian constraint, as mentioned in Sec. IV A. The positions of high-energy peaks agree with the expected energy within ±0.5 % [4]. In order to correct for the remaining residual energy shift due to systematic uncertainties, we included two fit pa- rameters, ∆E slope and ∆E offset . The energies of the simulated events E MC were then allowed to move as where ∆E was parametrized as with constrained parameters for slope and offset. Fig. 7 shows the fit to the blinded data, which is well described with χ 2 λ /ndf = 311/259 = 1.20. In the highstatistics region below 2800 keV, the residuals are symmetric and centered around zero with a standard deviation of σ res = 1.05. In the low-statistics region above 2800 keV, the fit lies mostly above the measured data leading to negative residuals and an asymmetric distribution, as the fit function can only predict rates larger than or equal to zero. The parameter pulls are shown in Fig. 8. None of the pulls for the blinded fit exceeds 2 σ and the sum of the squared pulls is 4.6. The pull on 60 Co is close to zero since its double-γ peak is located in the blinded region. Due to the degeneracy with the 226 Ra spectrum and its small background contribution, the 214 Bi component in the LXe shell outside of the TPC is less sensitive to the data than to its constraint and is not pulled away from the expected value. The parameters for the uranium and thorium chains are within the expected range. No notable pulls on the systematic uncertainty parameters are observed. The acceptance parameter prefers a value close to the lower bound.
In order to compute the sensitivity, a 0νββ signal was added to the background model as a Gaussian peak. Its mean µ 0νββ = (2457.8 ± 0.4) keV is given by the Qvalue [6,7], the standard deviation σ 0νββ = (19.7 ± 0.3) keV is given by the energy resolution. The SS fraction SS = 90.3 % of signal events was determined with MC simulations. Initial momenta for 10 6 electron pairs were generated with DECAY0 [26], their tracks were propagated with Geant4 and clustered based on the zseparation of subsequent energy depositions.
The expected sensitivity for setting a lower limit on the 0νββ half-life T 0νββ 1/2 of 136 Xe was determined using toy-MC simulations as outlined in Sec. II D. We derive a median upper limit on the decay rate with A 0νββ < 144 t −1 yr −1 at 90 % CL. Using Eq. (4), the expected sensitivity on the blinded data is T 0νββ 1/2, expected > 1.7 × 10 24 yr at 90 % CL.
D. Post-unblinding changes and final results After unblinding the events in the 0νββ ROI, an unexpected excess of events was observed around 2550 keV, well above the Q-value. This excess increased over time and was localized at the edges of the active volume. This indicated that an external background source progressively leaked into the selected data. Our investigation pointed to a class of MS events that were not rejected by the previously defined cuts, but that were misidentified as SS events. These events had a secondary S2 signal which was smaller than and temporally close to the main S2, and likely caused by multiple Compton scatters of a single γ-ray. As the secondary S2 contained a part of the total deposited energy, the misidentified population was reconstructed at a lower energy with respect to the γ-peak. The effect was present for all peaks in the ROI, but only the 208 Tl peak with its rising edge in the blinded region was large and isolated enough to significantly affect the 0νββ search.
The time-dependence of the effect is assumed to originate from increased PMT afterpulsing rates over time: MS events were identified based on the peak area and the top PMT hit pattern of the second largest S2 signal that was found in an event. PMT afterpulses that occurred in coincidence with the S2 altered the hit pattern as well as the signal size. Although the original MS classification accounted for the growth of afterpulsing with time, a stricter cut was needed to remove pathological waveforms such as the one shown in Fig. 9. A post-unblinding cut was introduced, based on the peak area of the secondary S2 and its top PMT array hit pattern. The effect of the cut is shown in Fig. 10.
The acceptance of the new cut was determined by comparing the number of events contained in the 208 Tl peak and multiple 214 Bi peaks outside of the previously blinded region before and after the new cut. MS events had a smaller reconstructed energy in the main S2 signal, as part of the energy was deposited in the subsequent S2 peaks. Thus, the centers of the γ-lines -especially of 208 Tl as the highest energy line -provided pure samples of SS. We found an acceptance of (97 ± 2) % and updated the total cut acceptance accordingly. Fig. 11 shows the fit of the combined signal and background model to the unblinded data after the addition of the new MS cut with χ 2 λ /ndf = 392/318 = 1.23. Residuals below 2800 keV are symmetric around zero, while the model is mostly above the data points at higher energies. The parameter pulls are indicated by the black bars in Fig. 8, and the best-fit numbers of background events around the Q-value are given in Tab  Contrary to the expectation the rate of 60 Co events is pulled close to zero. The individual 60 Co peaks are present in the data outside of the fitting range and the best-fit components of the other backgrounds do not point to an overestimation of the acceptance for SS events from the 238 U and 232 Th decay chains. This suggests that the strong pull is a feature of the stricter multi-scatter rejection. With its double MeV-γ signature, 60 Co is different from the other background peaks, which do not feature secondary γ-rays of equally high energies. In order to detect the 2505.7 keV peak as an SS event, the γ-rays need to be emitted in the same direction and fully absorbed within a few millimeters in x-y-z. This makes an SS reconstruction of these events less likely than for the other background sources. In the SS vs. MS classification of MC events, only the z-separation of consecutive energy depositions was considered. However, the MS selection on data also uses the S2 hit pattern which is sensitive to the x-y separation of multiple scatters. With the stricter post-unblinding cut, most 60 Co double-γ events were identified as MS and removed from the energy spectrum. Since this was not modeled in the MC, the background model fit results in zero 60 Co events.
The tendency of the energy shift parameters is changed in the unblinded fit since the previously remaining MS events in the low-energy flanks of γ-peaks biased the energy reconstruction. In the original fit, a stronger shift of the spectrum towards lower energies was observed which is not present after the removal of the events in the flanks. The best-fit cut acceptance at the Q-value is 85 %. This is consistent with the acceptance-loss attributed to the new MS cut. V. Expected and best-fit background event counts in the 2 σE ROI around Q ββ . The best-fit numbers are given for the fit to the data before unblinding, and for the fit to the unblinded data with the post-unblinding cut on multi-site events. The 60 Co upper limit is given at 90 % CL. The expected 214 Bi events in the TPC and LXe shell are given for an assumed activity concentration of (10 ± 5) µBq/kg. The standard deviation and mean position of the signal peak exhibit pulls close to zero. Here, the fit is more sensitive to the constraint than to the data in absence of a signal. The best-fit 0νββ decay rate is (65 ± 87) t −1 yr −1 which translates to < 210 t −1 yr −1 at 90 % CL. From Eq. (4), the lower limit on the half-life is T 0νββ 1/2 > 1.2 × 10 24 yr at 90 % CL.
The resulting effective neutrino mass range, using NMEs from [14,15], is m ββ = (0.8-2.5) eV/c 2 . Dedicated xenon-based 0νββ experiments such as EXO-200 and KamLAND-Zen [8,52] have reported results that supersede our result by up to two orders of magnitude. In contrast to XENON1T, which is optimized for low background in the keV region, the dedicated detector designs are optimized to have low backgrounds at Q ββ . Together with 136 Xe enrichment, this leads to a more favorable signal-to-background ratio. The previously most stringent limit from DM direct detection experiments was set by PandaX-II with T 0νββ 1/2 > 2.3×10 23 yr at 90 % CL [58]. XENON1T improves on this result by an order of magnitude with lower background, larger exposure and four times better energy resolution. It illustrates the potential of current and future DM experiments such as LZ [59], XENONnT [5], DARWIN [60] and beyond [61] for 0νββdecay searches.

V. SENSITIVITY OF XENONnT TO 136 XE NEUTRINOLESS DOUBLE-β DECAY
With its larger projected exposure, XENONnT will improve upon the XENON1T sensitivity to the 0νββ process. This section discusses the sensitivity projections based on simulated XENONnT background data.

A. Background model
The material backgrounds considered here originate from 60 Co as well as the 238 U and 232 Th decay chains. The strategy to simulate the background events followed Sec. II C, but employed the XENONnT Geant4 geometry [5]. We also took the background contribution from the NV PMTs into account, which was at the same level as subdominant intrinsic backgrounds.
The intrinsic background sources considered for XENONnT include the daughter nuclei of 222 Rn isotopes emanated from the detector materials and the 2νββ decay of 136 Xe naturally present in the LXe, both discussed in Sec. IV B. The expected 222 Rn contamination level of the LXe target is further reduced with respect to XENON1T with a novel online Rn removal system. The 222 Rn activity concentration was assumed to be 1.0 µBq/kg.
Due to the addition of 5.7 t of xenon to the total inventory, we conservatively assumed a natural abundance of 136 Xe in xenon, η136 Xe = (8.9 ± 0.5) × 10 −2 mol/mol, using the difference with respect to the XENON1T measured abundance as a systematic uncertainty. The actual isotopic composition will be measured in the future.
With a half-life of 3.82 min and a Q-value at 4.17 MeV, far beyond Q ββ , the β-decay of 137 Xe is a relevant background source in XENONnT. It is produced through neutron capture on 136 Xe occurring either within the TPC itself or in the non-shielded regions -outside the water tank -of the experiment, especially in the LXe purification systems. The mean travel time of xenon through the LXe purification system is approximately 7 min, which is short enough for 137 Xe to be injected back into the detector before decaying. We estimated a total 137 Xe production rate of (6±5) t −1 yr −1 through a MC simulation. Muon-induced neutrons produced in the LXe are primarily responsible for the production of 137 Xe in the TPC, contributing 10% of the total rate. The thermal neutron flux induced by radiogenic decays in rock, concrete, and  materials is the dominant contribution to 137 Xe creation in the purification system, contributing 90% of the overall rate. The main source of uncertainty stems from different measurements of the thermal neutron flux at LNGS that are in tension with each other [62][63][64][65][66]. Neutrino-electron scattering is a potentially irreducible background source for the 0νββ decay search if the incident neutrino flux and energy are sufficiently high. While the contribution from atmospheric neutrinos, diffuse supernova neutrinos, or geoneutrinos can be excluded, since either their flux or their energy is too small, the contribution from 8 B solar neutrinos is relevant. We used the neutral current 8 B neutrino flux measurement from the Sudbury Neutrino Observatory with Φ = (5.25 ± 0.16 stat ± 0.12 sys ) × 10 6 cm −2 s −1 [67] and the 8 B neutrino spectral shape from [68] to derive the expected rate of neutrino-electron scatters in the detector, following the neutrino-electron elastic scattering cross-section calculation from [69]. The electron neutrino survival probability follows the large mixing angle solution of the Mikheyev-Smirnov-Wolfenstein effect [70][71][72] which takes into account matter effects in the Sun.

B. Analysis
To maximize the sensitivity to the 0νββ process, we employed a FV optimization analogous to Sec. IV A in the 2 σ E region around Q ββ . We found an optimal FV with a mass of 1088 kg. The simulated energy spectra for all backgrounds within the FV are shown in Fig. 12. The dominant background contribution around Q ββ arises from the detector materials. Each isotope considered in the MC simulation has a constraint term arising from the radioassay measurements, as indicated in Sec. V A, or from dedicated MC studies of the background. The fitting model was developed as described in Sec. IV C; the 0νββ-signal was modeled as a Gaussian peak with an area proportional to the decay rate, A 0νββ .
Differences to the XENON1T fitting model arise from additional backgrounds considered here and from omitting the nuisance parameters related to shifts of the peak positions. Tab. VI summarizes the expected event yields for an assumed live time of 1000 days in the 2 σ E region around Q ββ . The dominant background contribution arises from the energy deposition of 214 Bi βand γ-emission where the β was absorbed in the passive detector materials. The contribution of this decay product of the 226 Ra decay chain, is dominated by radioactive decays in the cryostat.
The physics reach of XENONnT was estimated with a profiled likelihood approach [29] similar to the XENON1T analysis with the difference that we used the asymptotic assumption of test statistics whose validity was verified, see Sec. II D. We performed a likelihood scan for a set of assumed live times of the experiment (between 10 days and 1000 days) and determined the intersection with the 90 % quantile of a χ 2 -distribution with one degree of freedom. This is the expected median lower limit that we report in Fig. 13. The median lower limits for the considered live times were interpolated with a square-root function.  [29] with its 1 σ statistical uncertainty. The projected sensitivity and the observed results from XENON1T, KamLAND-Zen [8] and EXO-200 [52] are shown as solid and dashed lines, respectively.

C. Results
For a live time of 1000 days, we obtain a median lower limit of T 0νββ 1/2 > 2.1 × 10 25 yr at 90 % CL. (16) This is below the expected sensitivities and the observed lower limits from KamLAND-Zen and EXO-200 [8,52] due to the large background contribution from detector materials and the low 136 Xe abundance. Other LXe TPCs such as LZ [59] or DARWIN [60] are expected to be more sensitive due to a lower background level and larger FV mass, respectively. From the derived lower limit of the half-life, we computed the effective Majorana neutrino mass m ββ with the relation reported in Eq. (1). We used the same assumptions regarding the NME and the phase-space factor as laid out in Sec. IV D. We summarize our findings in Fig. 14, where the green band indicates the range of the effective Majorana neutrino masses for our XENONnT half-life sensitivity. The masses range from 0.19 eV/c 2 to 0.59 eV/c 2 depending on the NME. While XENONnT is not yet competitive with dedicated experiments, this study shows that future xenon DM detectors can be competitive with optimized high-energy backgrounds and larger exposures.

VI. CONCLUSION AND OUTLOOK
In this paper, we reported on searches for double weak decays of 124 Xe and 136 Xe with XENON1T. The search for 2νECEC decay included a larger sample of data and an improved signal model compared to our previous result [3]. We detect 2νECEC in 124 Xe with The search for 0νββ of 136 Xe is compatible with the background-only hypothesis with an exclusion limit of T 0νββ 1/2 > 1.2 × 10 24 yr at 90 % CL. Due to a larger active mass and expected lower background rate, the XENONnT experiment will improve this result. With a live time of 1000 days, we expect a median lower limit of T 0νββ 1/2 > 2.1 × 10 25 yr at 90 % CL. While this is not competitive to dedicated searches, it demonstrates the feasibility of more sensitive searches in future xenon DM detectors.
The study of double-weak processes in LXe TPCs is not restricted to the two analyses discussed in this work and can be extended to a plethora of rare decays, such as the search for the 2νββ decay of 136 Xe to the 0 + 1 excited state of 136 Ba [76], the 2νββ and 0νββ decay of 134 Xe [77] or the neutrinoless second-order weak decays of 124 Xe [51]. Furthermore, a precise measurement of the 2νββ energy spectrum offers the possibility of experimentally testing the underlying nuclear models [78,79], but also to probe new physics beyond the SM [80][81][82]. The XENON project provides a broad science program ranging from DM searches to neutrino physics and properties of xenon, covering several orders of magnitude in energy.
gratefully acknowledge support from the National Science Foundation, Swiss National Science Foundation, German Ministry for Education and Research, Max Planck Gesellschaft, Deutsche Forschungsgemeinschaft, Helmholtz Association, Dutch Research Council (NWO), Weizmann Institute of Science, Israeli Science Foundation, Fundacao para a Ciencia e a Tecnologia, Région des Pays de la Loire, Knut and Alice Wallenberg Foundation, Kavli Foundation, JSPS Kakenhi in Japan, Tsinghua University Initiative Scientific Research Program and Istituto Nazionale di Fisica Nucleare. This project has received funding/support from the European Union's Horizon 2020 research and innovation programme under the Marie Sk lodowska-Curie grant agreement No 860881-HIDDeN. Data processing is performed using infrastructures from the Open Science Grid, the European Grid Initiative and the Dutch national e-infrastructure with the support of SURF Cooperative. We are grateful to Laboratori Nazionali del Gran Sasso for hosting and supporting the XENON project.  IX. Parameter constraints and best-fit parameters for the XENON1T 0νββ search. Best-fit parameters are given for the blinded and unblinded data. Unitless parameters state the relative change to the expected value. For the acceptance, the constrained and best-fit acceptances are given instead of the fit parameters. Due to the implementation of the acceptance, the Gaussian constraint on the acceptance scaling parameter yields an asymmetric acceptance range. Due to the internal processing of the thorium chain, the respective best-fit parameters cannot be translated directly to 228 Th event counts in Tab. V.