Search for two-neutrino double-beta decay of ¹³⁶Xe to the ${\bf 0^+_1} $ excited state of ¹³⁶Ba with the complete EXO-200 dataset

Double-beta ($ \beta\beta $) decay is a second-order weak transition in which two neutrons simultaneously decay into two protons. The observation of the decay mode without neutrino emission (0νββ) would demonstrate the Majorana nature of neutrinos, which is a fundamental question in particle physics. The current best half-life limits on 0νββ set by various experiments [1−5] for different isotopes are in the range of 10^25-26 yr. The $ \beta\beta $ mode with two accompanying neutrinos (2νββ) is a standard model process and has been observed in over ten isotopes [6] with half-life between 10¹⁹ and 10²⁴ yr. The 2νββ can also decay to the excited states of the daughter nucleus. These $ \beta\beta $ modes share the same initial and final nucleus as the 0νββ decay, and might have correlated nuclear matrix elements (NMEs) [7, 8]. Measuring the decays to the excited states offers additional experimental input to the calculation of 2νββ NMEs, which might help reduce the theoretical uncertainties for 0νββ NMEs.

The 2νββ decay to the excited states of the daughter nucleus are suppressed by orders of magnitude with respect to the decay to the ground state because of the reduced phase space from smaller Q values [9]. In addition, the decays to the $ 2^+ $ states are highly suppressed by angular momentum. Therefore, the decay to the first excited state of the $ 0^+ $ state, denoted as $ 0^+_1 $ hereafter, is the most viable. The first experimental limit derived for the decay to the excited state dates back to 1977 for ⁷⁶Ge [10]. The first positive signal was observed for ¹⁰⁰Mo decay to the $ 0^+_1 $ excited state of ¹⁰⁰Ru in 1995 [11], followed by confirmations from a series of other works [12, 13]. The ¹⁵⁰Nd decay to the $ 0^+_1 $ excited state of ¹⁵⁰Sm is the only other isotope where a positive signal has been observed [14]. Stringent limits have been set for other isotopes [15]. The CUORE experiment placed a limit of $ T^{2\nu}_{1/2}(0^+ \rightarrow 0^+_1) > 2.4 \times 10^{23} $ yr for ¹³⁰Te in 2019 [16]. More recently, the MAJORANA-DEMONSTRATOR placed a limit of $ T^{2\nu}_{1/2}(0^+ \rightarrow 0^+_1) > 7.5 \times 10^{23} $ yr for ⁷⁶Ge [17]. KamLAND-Zen and EXO-200 performed measurements for ¹³⁶Xe [9, 18], and the current best limit is set at $ T^{2\nu}_{1/2}(0^+ \rightarrow 0^+_1) > 8.3\times10^{23} $ yr.

The ¹³⁶Xe $ \beta\beta $ decay to the ground state of ¹³⁶Ba has a Q value of 2457.83 keV [19]. The $ \beta\beta $ decay to the $0_1^ + $ state of ¹³⁶Ba has a Q value of 878.8 keV, followed by the emission of two characteristic de-excitation γ's with energies of 760.5 keV and 818.5 keV from consecutive de-excitations from the $ 0^+_1 \rightarrow 2^+_1 \rightarrow 0^+ $ states of ¹³⁶Ba [19]. Two electrons accompanied by two γ's provide distinctive event signatures. The theoretically predicted half-life for 2νββ decay of ¹³⁶Xe to $0_1^ + $ excited state of ¹³⁶Ba ranges from 10²³ to 10²⁶ yr using different nuclear models [20, 21]. The sensitivities of current generation experiments are approaching this level, presenting a chance for discovery. In this paper, we report a new search for the decay to $0_1^ + $ state of ¹³⁶Ba using the complete EXO-200 dataset.

The EXO-200 detector is a single phase liquid xenon (LXe) time projection chamber (TPC). The experiment was located in the Waste Isolation Pilot Plant (WIPP) near Carlsbad, New Mexico (USA), with an overburden of $ {1624}^{+22}_{-21} $ meters of water equivalent [22]. EXO-200 operated between Sept. 2011 and Dec. 2018 in two phases. Phase I data taking was halted in Feb. 2014 due to underground incidents at the host facility unrelated to the experiment. Phase II operation began in May 2016 after underground access was restored. Among other things, the detector electronics were upgraded before Phase II to reduce coherent noise on the avalanche photodiodes (APDs) during the forced outage. A brief review of the detector is given in this section. More details on the experiment can be found elsewhere [23−25].

The EXO-200 TPC was housed in a cylindrical, thin-walled copper vessel with a cathode in the middle, resulting in two back-to-back drift regions with ~18 cm radius and ~20 cm drift length. Two wire planes crossing of 60° are placed at each anode end of the drift. The induction plane (V-wires) records the transit of charge, while the collection plane (U-wires) collects the charge. The scintillation photons produced simultaneously with the ionization were detected by arrays of large area avalanche photodiodes (LAAPDs) [26] behind the anode planes. The entire copper vessel enclosing the TPC was submerged in HFE-7000 cryofluid [27], enclosed by additional layers of passive shield made of 5.4 cm copper and 25 cm lead in all directions [23]. An active muon veto was composed of scintillator panels, covering four sides of the detector, providing >94% tagging efficiency for muons passing through the TPC [22].

Event reconstruction makes use of both the light and charge signals produced by particle interactions in the LXe. Signals from U-wires and V-wires are used to extract the position of the energy deposit in the plane normal to the electric field direction. The time difference between the registered light signal on the LAAPDs and the charge collection on U-wires is used to obtain the z position along the drift field, using the measured electron drift velocity [28]. This allows full three dimensional (3D) vertex reconstruction for individual energy deposits, defined as charge clusters in the analysis to distinguish them from the true energy deposits. Small charge deposits with energy below the V-wire but above the U-wire detection threshold do not have a reconstructed $ xy $ position. Events reconstructed with multiple charge clusters are referred to as "multi-site" (MS) while single charge cluster events are denoted as "single-site" (SS). The SS/MS identification can help distinguish $ \beta\beta $ events from γ's, as β-like events are primarily SS, while γ's tend to be MS. Only 12% (14%) of γ events are reconstructed as SS events near the $ {\rm{Q}}_{\beta\beta} $ of ¹³⁶Xe for the ²²⁸Th (²²⁶Ra) source [1].

The event energy is reconstructed by a linear combination of the light energy measured by the LAAPDs and the charge energy measured by the U-wires to fully exploit the anti-correlation between these two channels [29]. Such a linear combination cancels out anti-correlated fluctuations in light and charge signals, thus optimizes the energy resolution, which is especially important for the 0νββ search.

A comprehensive Monte Carlo (MC) detector simulation based on GEANT4 [30] was developed to model the responses to various signal and background interactions. This simulation models the entire chain from energy deposits produced by interactions in the LXe to waveforms created on the crossed-wire planes by the ionization propagating through the detector. The simulated waveforms are treated the same way in the reconstruction and analysis framework as the real data. The simulation is benchmarked by external γ calibration sources located ~10 cm away from the fiducial volume (FV) at set positions near the cathode and the anodes [24].

The dataset used in this analysis is the same as that used in [1], with a total livetime of 1183.1 days after run quality selections. All clusters of each event are required to lie within a FV, defined as a hexagon in the $ xy $ plane with an apothem of 162 mm and more than 10 mm away from the cylindrical polytetrafluoroethylene reflector, as well as the cathode and the V-wire planes. This FV contains $ 3.31\times10^{26} $ atoms of ¹³⁶Xe, with an equivalent mass of 74.7 kg. The systematic error of this value is included in the common normalization error, which will be discussed in Sec. V. After additional muon veto cuts to reduce cosmogenic events, the total ¹³⁶Xe exposure is 117.4 (116.7) kg·yr in Phase I (Phase II).

Compared to the analysis based on the Phase-I dataset in [9], a major improvement to the signal efficiency is achieved by relaxing the 3D cluster reconstruction requirement. Because the signal detection energy threshold on the V-wires (~200 keV) is higher than that on the U-wires (~90 keV), a large number of charge clusters with small energy deposits have no V-wire signals, resulting in incomplete $ xy $ positions for these clusters. In the previous search, all clusters were required to have fully reconstructed 3D positions. For this search, we only require the sum energy of all fully reconstructed clusters to be above 60% of the total event energy [1]. The efficiency for signal events increased by a factor of ~2, while the total background events only increased by a factor of ~1.4. In addition, the newly developed background discriminator, utilizing event topologies further improves the signal to background ratio. With a relaxed 3D vertex reconstruction cut, the signal efficiency increased from 24.5±3.1% (25.4±3.1%) to 58.0±2.9% (58.7±2.9%) for Phase I (Phase II). The remaining signal inefficiency is from the fiducial volume and energy cuts. Lowering the energy threshold has only a marginal gain in signal efficiency, which is out weighed by the large increase in background efficiency at those energies. As a result, there is little improvement to the sensitivity from lowering the energy threshold and the same 1000 keV threshold from the previous analysis is used.

A deep learning (DL) method based on a convolutional neural network (CNN) [31] is used to separate excited state signals from backgrounds in this work. Unlike the boosted decision tree (BDT) approach used in [9] where cluster information is combined into high level variables, the individual cluster information is directly provided as input to the CNN. This input contains more complete information about the event, which further improves background discrimination. The details will be discussed in Sec. IV.

The accurate determination of the energies for clusters is essential to identify the characteristic de-excitation γ's. The cluster energy is reconstructed from charge signals only, because the light signals are unresolvable among clusters in an event as they appears as one flash of light in the detector. Previous EXO-200 work already showed charge yield in liquid xenon is energy dependent in the MeV range [32]. A calibration curve is derived using SS events from ¹³⁷Cs, ⁶⁰Co, ²²⁶Ra and ²²⁸Th calibration sources. Using all available calibration data spanning the experiment lifetime, the calibrated cluster energy scale uncertainty is determined to be within 1% above 662 keV.

The fitting framework closely follows previous analyses [1, 33]. The signal search was performed with a negative log-likelihood (NLL) function to fit simultaneously the SS and MS events to two dimensional (2D) Probability Density Functions (PDFs) of event energy and a DL discriminator for different background and signal components. The relative fraction of SS and MS events for each component is set by MC predicted values and constrained by errors determined from data-MC difference in calibration source data (Sec. V). PDFs are constructed based on simulation validated by extensive calibration data. Residual discrepancies between data and simulation are taken into account as systematic uncertainties. Gaussian terms are added to the NLL to take into account various systematic errors, which will be discussed in Sec. V.

A profile likelihood scan was performed to derive the 90% confidence level (CL) limit at the negative logarithm of the likelihood ratio between a given number of counts and the best fit (ΔNLL) value of 1.35, under the assumption of Wilks’ theorem [34, 35] considering the large statistics of the dataset in the region of interest. A sensitivity study is performed to compare different analyses. The sensitivity is evaluated by generating a set of toy datasets based on the background model derived from a fit to the energy spectrum of the low background data without the excited state signal included in the PDFs. It represents our best $ a\; priori $ understanding of the background model. The 90% CL limit is derived for each individual toy dataset, and the median of the limit distribution from all toy datasets is defined as the sensitivity. Though the full EXO-200 dataset has been unblinded previously, we avoided tuning the analysis against the low background dataset used to search for excited state signals. The analysis strategy is chosen prior to the final fit based on these sensitivity studies which only rely on calibration data and simulations.

The excited state signals have distinct features due to the two accompanying γ′s of specific energies. The individual cluster energy and position variables for each event are used as inputs for background discrimination, as they are expected to maintain more information about the energy and topology of the β and γ related interactions. A DL based method was utilized to extract the correlation between both energy and position maintained by the two de-excitation γ′s from this more complicated input. The details are given in this section.

To train the background discriminator, the background model obtained from an energy-only fit to the low background data is used to represent the background compositions in the data. The relative fractions of different components of the background model are given by the best fit values. In total, 1.8 million MC events are used to train the network, composed of half signal events and half background events. During the training process, 80% of the sample was used for training and the remaining 20% for validation, Training and validation of the background discriminator is done in python using the Pytorch package [36].

The input information for the network are the energy and 3D position ($ x, y, z, E $) for every reconstructed cluster of the event, as shown in Fig. 1. As the number of charge clusters (i.e. event multiplicity) varies event by event, the dimension of the input matrix is set to 20×4, with 20 set to safely allow for the maximum number of clusters in any event. For events with multiplicity less than 20, the remaining rows are padded to zeros. The non-zero rows are arranged according to descending cluster energy order. The normalization of the input is done using the linear normalization formula $ v_{\rm{norm}}=v/(v_{\rm{max}}-v_{\rm{min}}) $, where v represents the $ (x, y, z, E) $ in each cluster. The maximum and minimum values are set by their corresponding physical limits in each variable so that the values of $v_{\rm norm}$ for different variables are comparable. The normalization was found to improve training stability and lead to more separated peaks in both the signal and background discriminator distributions, though there was negligible impact on the receiver operating characteristic (ROC) curve.

Figure 1.  (color online) Example of network input for a $0_1^ + $ excited state signal (left) and a 2νββ background event (right). The input image is fixed format of twenty rows by four columns. The four columns represent ($ x, y, z, E $) accordingly while each non-zero row represents a charge cluster. Incomplete $ xy $ information are set to placeholder values (–2.498 after normalization). The energy and spatial correlations among clusters is retained by this input with all available energy and vertex information for each cluster.

Clusters with no reconstructed $ xy $ positions are allowed in the analysis, with corresponding $ xy $ values set to an unphysical placeholder values (–999) before the normalization operation mentioned above. These placeholder values were found not to affect or bias the network performance. This is validated with two networks prepared by training on two different samples: one with only full 3D events, the other also including events with incomplete $ xy $ positions. The ROC area for the two networks on the same test sample of full 3D events are very similar. The ROC area of the latter network for events without full 3D positions is only slightly worse than full 3D events, presumably due to incomplete position information. The distributions of the ($ x, y, z, E $) for signal and background are shown in Fig. 2.

Figure 2.  (color online) Distributions of the input variables to the discriminator based on MC simulations in Phase II. The information of all clusters in the MS events is filled in the histogram. Useful information for background separation like number of clusters is not easily seen in 1D projections, but they are embedded in the image input.

CNN's are one of the most commonly used DL methods in high energy physics in recent years, with applications in particle identification and event reconstruction dealing with complicated input information in the format of images. EXO-200 has used this in the 0νββ analysis for background discrimination [1], as well as analyses for event reconstruction [37]. In this analysis, a CNN with a simple structure inspired from TextCNN [38] was developed to take the 20×4 array of cluster information. The architecture of the network is shown in Fig. 3. The network consists of a convolution layer, a max pooling layer and fully connected layers. The convolution part is composed of convolution kernels in six different sizes. The kernels have the same number of columns as the network input. The convolution is only done along the row dimension, each generating a one-dimensional array. A max pooling layer is applied, followed by two fully connected layers. The final output value after a sigmoid function gives the event a score indicating the type of event.

Figure 3.  (color online) The CNN architecture used for background discrimination. The convolution layer consists of six different kernel sizes. The kernel shape is n×4 with n ranges from 1 to 6. There are 100 convolution kernels for each size. A max pooling layer and two fully connected layers are followed.

The test sample has distributions consistent with the training sample. This implies no over-training of the network. The ROC curve, background rejection power vs signal efficiency, is shown in Fig. 4. The two phases give comparable background rejection.

Figure 4.  (color online) ROC curves of CNN discriminators in Phase I and Phase II. The Area Under Curve (AUC) is very similar between the two phases. The inset shows the distribution of the CNN discriminator for signal (red) and background (blue) events in Phase I as an example. The background rejection efficiency shown in this plot comes from information including event topology and spectrum shape. More details on spectrum shape effect can be seen in [1, 37, 39].

To understand the residual background contributions, the energy spectrum before and after a background cut on the discriminator variable is shown for illustration in Fig. 5. 2νββ decay to the ground state of ¹³⁶Ba dominates the low background data, but it can be rejected with very high efficiency using event topology information. γ′s from ²³⁸U, ²³²Th, ⁴⁰K and ⁶⁰Co, though much lower in rate, are harder to reject as they produce more clusters resembling $0_1^ + $ state signals. An obvious energy dependence on the discrimination was observed as well. Low rejection power is found for events around 1800 keV, where the broad spectrum of the excited state peaked. As a result, ²³⁸U-like backgrounds remain as a dominant background after background rejection. ⁶⁰Co is difficult to reject due to the two γ's, leading to higher multiplicity values, like for the signals.

Figure 5.  (color online) Background spectrum for MS events before (top) and after (bottom) a selection cut at a signal efficiency of ~30% in Phase I. The simulated spectrum for signal and background components are shown. The background rates are set by the background model described in Sec. III. For the 2νββ decay of ¹³⁶Xe to excited states, the number expected from the 90% CL sensitivity is plotted (will be discussed in Sec. VI). Other backgrounds consists of ⁴⁰K, ⁶⁰Co, ¹³⁵Xe, ¹³⁷Xe and background related to neutron capture.

Possible mismodeling of the spectral shape of the discriminator in Monte Carlo is a major source of systematic uncertainty. The agreement is studied by various calibration sources as shown in Fig. 6. The data is chosen to be binned with three equal efficiency bins for excited state signals. More bins are found to improve the background separation power, but with the risk of sacrificing data simulation agreement. Additionally, due to the powerful rejection ability, the background distribution decreases drastically towards a CNN value equal to 1. Therefore, fewer bins ensure enough statistics in the bin close to 1 for calibration data to constrain data simulation agreement. Data and MC agree within 15% for all source positions. The study of the binning method used only simulation and calibration data. The error introduced to the excited state signals from the residual disagreement will be discussed in Sec. V.

Figure 6.  (color online) Shape agreement of CNN between data (dotted line) and MC (Solid line) for MS events using calibration sources positioned near the cathode in Phase I (left) and Phase II (right). The choice of binning allows equal number of excited state signal events in each bin (middle and bottom panels). The shape agreement of uniform binning is also shown in the top panel as a comparison.

The impact of systematic uncertainties is accounted for by adding Gaussian constraints to the NLL. The list of systematic uncertainties are discussed below:

1. A common normalization error, caused by uncertainties in event reconstruction and selection efficiencies, is applied to all PDF components equally.

2. SS fraction, due to uncertainties in the relative fraction of SS events out of all the events. The mean SS fractions are derived from by simulation, and their error is determined by simulation and benchmarked with calibration data.

3. Uncertainty in the relative fraction of neutron capture related PDF components by dedicated simulations.

4. Uncertainty in the activity of radon in the LXe as determined in the standalone studies via measurement of ²¹⁴Bi-²¹⁴Po correlated decays.

5. Signal specific normalization error. An error associated only to signal events, allowing the signal to vary by the estimated error.

The first four errors were evaluated in previous EXO-200 analyses [1]. The common normalization errors are 3.1% (2.9%) for Phase I (Phase II), with the dominant contribution from the fiducial volume cut [1]. The SS fraction error is evaluated by the difference between data and MC simulation for various calibration sources at different positions. The errors are 5.8% (4.6%) for Phase I (Phase II) [1]. The relative capture fractions of cosmogenic neutrons is constrained with a 20% uncertainty [22]. The radon daughters-induced background in LXe is constrained by the measured rate of radon decays [24].

The estimate of the signal-specific normalization error follows the same methodology as in [40]. The main difference from the result in [9] is that a signal count dependent treatment is used in this analysis instead of a constant fractional uncertainty. This more accurately accounts for the error at small signal counts. This error varies as a function of signal counts, and consists of two main contributions: 1) shape error, caused by the level of the PDF shape agreement between data and MC; 2) the background model error due to not considering all the detailed locations of backgrounds originating from materials far from the TPC Vessel in the fit model. Instead, some representative positions are used to represent ²³⁸U, ²³²Th, ⁶⁰Co from these materials. The systematic errors caused by this approximation are estimated by replacing the PDF of the remote components at different locations in the fit. The shape errors are evaluated on an ensemble of toy datasets. Each toy dataset was generated from the MC PDFs weighted by the observed data/MC ratio based on the calibration data, but fitted with the original unweighted PDFs. The difference between the injected number of signals against the fitted number of signals is taken as the shape error. The background model errors are evaluated by comparing the difference of best fit signal counts by replacing a PDF component in the background model with its alternative one at a different position. The two contributions are added in quadrature in the end, with the shape error being the dominant one. The evaluated signal normalization errors ($\sigma_{\rm signal}$) at different injected signal numbers are found to be well described by $\sigma_{\rm signal}/N=a/N$, with N being the signal counts and a being the parameter used to quantify signal-specific normalization error. The errors are summarized in Table 1.

A possible energy scale difference between beta particles and gamma particles is considered. The energy scale of beta-like events is allowed to float freely with respect to γ events by multiplying a beta scale factor to all PDFs representing interactions of β-like events in the fit. The best fit value of beta scale is 1.0017±0.0017 (1.0008±0.0017) for Phase I (Phase II), suggesting a consistent energy scale within subpercent level above the 1000 keV analysis threshold.

The 90% CL sensitivity to the half-life of the excited state decay was evaluated to be $ 2.0\times10^{24} $ yr for Phase I. With the improved systematic uncertainty from the CNN discriminator, optimized selection cuts and slightly larger exposure, the Phase I sensitivity is improved by 15% from the BDT-based approach in [9], under the new treatment of signal dependent normalization error. The new analysis of Phase II data presented in this work has a slightly better sensitivity of $ 2.2\times10^{24} $ yr because of smaller systematic uncertainties. Considering the current sensitivities are dominated by statistical uncertainties, the combined sensitivity can be calculated by treating the systematic errors between the two phases as independent, which gives a combined sensitivity of $ 2.9\times10^{24} $ yr.

A final fit using energy and CNN as fitting dimensions was applied to the full EXO-200 dataset. We found no statistically significant signals in either phase (Fig. 7). A lower limit on the half-life is obtained to be $0.9 \times10^{24} $ yr and $1.4 \times10^{24} $ yr for Phase I and Phase II. The combination of the two phases gives a limit of $1.4 \times10^{24} $ yr. While there are large uncertainties from different nuclear models, the result in this work is in tension with the values predicted by QRPA as summarized in Table 2. The future nEXO experiment [41] is expected to greatly improve the search capabilities since it is expected to fully contain the de-excitation gammas and allow lower backgrounds, as well as have much more exposure [42].

Figure 7.  (color online) Best fit to the MS low background data for energy spectrum (left) and discriminator spectrum (right) in Phase I (top) and Phase II (bottom). The energy bins are 15 keV and 30 keV below and above 2800 keV, respectively. The best-fit residuals of the MS energy spectrum are shown for illustration, with only statistical uncertainty taken into account. The small deviations are taken into account in the spectral shape systematic errors.

		Reference	$ T^{2\nu}_{1/2} $ ($ 10^{23} $ yr)
Theory	QRPA	[20]	$ 0.14-13 $
	QRPA	[21]	$ 1.3-8.9 $
	IBM-2	[20]	$(1.5-3.6)\times 10^{2}$
	IBM-2	[43]	$ 2.5 \cdot10^2 $
	NSM	[20]	$(2.5-6.6)\times 10^3$
	EFT	[20]	$(0.62-16)\times 10^{2}$
Experiment	KamLAND-Zen	[18]	$>8.3 $
	EXO-200 (2016)	[9]	$>6.9 $
	EXO-200	This work	$>14 $

Table 2.  Theoretical and experimental results of $ ^{136}{\rm Xe} $ $ 2\nu\beta\beta $-decay half-life to the $ 0^+_1 $ state of $ ^{136}{\rm Ba} $.

The results of a refined search for 2νββ of ¹³⁶Xe to the $0_1^ + $ excited state of ¹³⁶Ba using EXO-200 are reported in this paper. No statistically significant evidence for this process is found and a limit on the half-life of $ T_{1/2}>1.4\times10^{24} $ yr at 90% CL is obtained. A CNN based discriminator was utilized in this analysis, which fully exploits the cluster information for background rejection while achieving good agreement between data and simulation, leading to an improvement by a factor of 1.7 relative to the current best constraint set by KamLAND-Zen. Future ton-scale experiments [4, 41] might have a chance to make a discovery.

We gratefully acknowledge the KARMEN collaboration for supplying the cosmic-ray veto detectors, and the WIPP for their hospitality.