Signal identification with Kalman Filter towards background-free neutrinoless double beta decay searches in gaseous detectors

Particle tracks and differential energy loss measured in high pressure gaseous detectors can be exploited for event identification in neutrinoless double beta decay~($0\nu \beta \beta$) searches. We develop a new method based on Kalman Filter in a Bayesian formalism (KFB) to reconstruct meandering tracks of MeV-scale electrons. With simulation data, we compare the signal and background discrimination power of the KFB method assuming different detector granularities and energy resolutions. Typical background from $^{232}$Th and $^{238}$U decay chains can be suppressed by another order of magnitude than that in published literatures, approaching the background-free regime. For the proposed PandaX-III experiment, the $0\nu \beta \beta$ search half-life sensitivity at the 90\% confidence level would reach $2.7 \times 10^{26}$~yr with 5-year live time, a factor of 2.7 improvement over the initial design target.


Introduction
Neutrinoless double beta decay (0νββ) is a hypothetical weak decay process that would confirm the Majorana nature of neutrinos and provide a direct evidence of lepton number violation [1,2]. Experimental search for 0νββ has been an active frontier in particle and nuclear physics [3]. The Standard-Model-allowed two neutrino double beta decay (2νββ) has been observed in more than ten isotopes. Searches for the 0νββ utilize those isotopes, including 136 Xe, 76 Ge, and 130 Te. The current lower half-life limits for 0νββ of the three isotopes are 1.07 × 10 26 yr, 1.8 × 10 26 yr, and 3.2 × 10 25 yr (90% confidence level, or C.L.), established by KamLAND-Zen, GERDA, and CUORE experiments respectively [4][5][6]. Most 0νββ experiments identify possible signals by event excess around the decay Q-value in the summed electron energy spectrum.
Xenon-based high pressure gaseous Time Projection Chambers (TPCs) can utilize topological features of event tracks and differential energy loss to identify possible 136 Xe 0νββ signals. Two emitted electrons carry the Q-value of 2458 keV and may travel O (20 cm) along meandering tracks in 10 bar xenon gas. Signal identification with tracks in a 136 Xe gaseous TPC has been exploited early on by the Gotthard experiment [7]. More recently, the NEXT experiment [8] demonstrated the effective track reconstruction with electroluminescence amplification signals in their prototype TPC [9]. Efficiencies of signal and background identification in the proposed PandaX-III gaseous TPC [10] have been studied with simulation data [11,12]. The aforementioned experiments utilize the prominent Bragg Blob (BB) feature of event tracks. Each end of a 0νββ track has a BB, in which rapid energy loss in unit volume happens because of increased differential energy loss (Bragg peak) and larger scattering angles right before an electron stops. However, only one BB are on the ends of the track to the background. Fig. 1 shows typical tracks of 0νββ signal and background events from simulation in gaseous xenon TPC, where the track of 0νββ has two BBs while that of background from the 238 U decay chain has one. Our study aims to reconstruct the meandering tracks with Kalman Filter in a Bayesian formalism (KFB) [13,14] and extract more signal-identifying features besides BB [7,9]. For example, differential energy loss along tracks and particle momentum by the end of the tracks can be calculated and used for the signal and background identification. When combined with parameters such as track energies, we demonstrate that almost all the background in the region of interest (ROI) of 0νββ search can be rejected and the search sensitivity further improved in a typical low-background gaseous TPC.
The paper is organized as follows: In Section 2, the generation of simulation data is presented. Five detector configurations are studied to compare the effect of the readout schemes and detector energy resolutions. In Section 3, we show the preprocessing and the method of track reconstruction with KFB in details. Then topological parameters are extracted for signal identification in Section 3.3. In Section 4, the signal and background discrimination power and the improvement for 0νββ search sensitivity of a proposed detector are presented. At last, we discuss the potential improvement and broader applications of KFB.

Simulation
To demonstrate the performance of the KFB approach, we construct a conceptual high pressure gaseous TPC with the Geant4 simulation framework [15]. The signal events are the 0νββ of 136 Xe. Background events from the decay chains of 232 Th and 238 U are considered in the simulation. To explore the impact of event discrimination power of different readout schemes and detector energy resolution, five detector configurations are considered.

Geometry and event generation
The detector geometry is similar to the PandaX-III TPC as outlined in Ref. [16]. The active volume (AV) is 1.6 m in diameter and 1.2 m high, which contains approximately 140 kg of xenon gas (with 90% 136 Xe) at 10 bar. The cathode and the readout plane are placed at the two bases of the cylindrical AV respectively. Outside of the AV, we construct an acrylic field cage, copper shielding liner, stainless steel vessel, lead shielding, and highdensity polyethylene shielding in sequence. The dimensions of each component are identical to Ref. [17]. The 0νββ signals are produced with the Decay0 package [18], which gives the energy distributions of two emitted electrons. Background events from the decay chains of 232 Th and 238 U in the detector components and shielding layers are considered.
The detector response of TPC, including electron diffusion, readout schemes, and energy resolution, is simulated in the REST framework [12]. While drifting to the readout plane, ionization electrons diffuse transversely and longitudinally which results in broadening of the tracks. The transverse (longitudinal) diffusion coefficient is set to be 1.0 (1.5)×10 −2 cm 1/2 , assuming xenon is mixed with a quencher gas such as Trimethylamine to have reduced diffusions. Detector response blurs event tracks spatially and energetically, which decreases the discrimination power between signal and background events. Hence, it is critical to reconstruct the tracks accurately for effective background suppression.

Detector configurations
We have performed our studies under five different configurations by varying the readout schemes and detector energy resolutions.
In the first two high granularity configurations, the detector's readout plane is fully instrumented with 1 mm × 1 mm pixels. An energy resolution of 3% (Full Width at Half Maximum, FWHM) at the Q-value is assumed for the first one and 6% for the second. We compare such detectors with configurations of degraded spatial granularity of 3 mm × 3 mm and energy resolution of 3% and 1% respectively in the third and fourth configurations.
The last configuration replicates the PandaX-III detector specifications. The readout plane is covered with 52 pieces of 20 cm × 20 cm readout modules, each of which is equipped with 3-mm-wide readout strips. Each strip reads out signals from 64 connected pixels in horizontal or vertical directions. The setup significantly reduces the number of readout channels needed, but does introduce ambiguity in track reconstruction. It is also worth noting that efficiency loss is considered because the TPC's AV is not 100% monitored by readout modules in this configuration. The energy resolution is assumed to be 3%.
Later we will refer the configurations as (1 mm, 3%), (1 mm, 6%), (3 mm, 3%), (3 mm, 1%), and (3 mm strip, 3%) respectively. Once reaching the readout plane, ionization electrons register charge signal hits (or hits for short) in pixels or strips. The amplitude and timing of hits contain all the information that can be collected in a physical detector. The output hit data from our simulation would mimic detector data.

Track reconstruction with KFB
Kalman Filter is widely used in particle physics experiment [19][20][21]. It is used as an optimal estimator for track reconstruction which combines information from physical model prediction and measurement data. Besides track fitting, it can also be used for optimum extrapolation, identification of pseudo-points, and error adjustment [22]. The KFB approach combines Kalman filtering with the Bayesian formalism and can estimate the noise covariance matrix of model prediction and measurement. Focusing on the signal identification in 0νββ search, we develop a new method based on KFB to reconstruct the meandering tracks of MeV-scale electrons. The hits of the simulation data are grouped into principal and subordinate tracks firstly in the preprocessing steps, then the principal track is reconstructed with KFB.

Preprocessing
The 0νββ signal and background events may generate more than one track in gaseous TPC. We focus on the principal track which carries the most deposited energy and contains the most information. The other shorter tracks (if exist) are named subordinate track(s). The simulation data are preprocessed in two phases before reconstruction with the KFB approach.
Firstly, we group hits of an event into a principal track and subordinate track(s) based on the distance among the hits. The grouping is implemented with a widely-used clustering algorithm called DBSCAN [23], which clusters discrete hits based on proximity. The cluster with the most deposited energy would be selected as the principle track. An example of well-separated principal and subordinate tracks of a background event is shown in Fig. 2a.
Secondly, the principal track is reconstructed roughly. Hits of the principal track are divided into segments with the Birch clustering algorithm [24], shown as clusters of colored dots in Fig. 2b. The charge-weighted centers of segments, named Birch clusters (BCs), are shown as red dot dash lines in Fig. 2b. They are sorted by a modified Ant Colony Optimization algorithm [25] with random starting points. The algorithm enumerates different connections of BCs and finds the one with the shortest total track length of all enumerations. The number of segments is optimized for the balance of sorting quality and computing loads.

Kalman filter
The Kalman filter is then introduced to further refine the reconstruction of the tracks [22,26], which keeps on recursion of the covariance to estimate the optimal position of each sampled hit. We use the classical Kalman filter which describes a linear dynamic system  Figure 2: Track reconstruction procedures of simulated 0νββ signal (left side) and background (right side) events, including (a) identification of the principal track, (b) rough reconstruction, and (c) reconstruction with KFB, respectively. The tracks are projected in the XZ plane, where X is one of the transverse directions and Z is the longitudinal direction along the field lines in a TPC. The MC-truth track from Geant4 (black dashed lines), simulated readout hits (orange points), Birch clustering hits (multi-colored points), BCs (red dot dash lines), and the reconstructed tracks by the ends (red lines) are labeled. We filter on both sides of the track, but one of the filter result is drawn for illustration. KFB stops at the BCs on each end. with given process and measurement noise [27]. Two equations, the prediction equation and measurement equation need to be determined firstly.
For the state prediction, we define a state vector s = (x, y, z, u x , u y , u z ) T at any given hit, including the three dimensional position and unit velocity vector. Considering the physical process of multiple scattering between two adjacent steps k − 1 and k, the prediction equation can be written as: where F is the propagation matrix of uniform rectilinear motion and ω k is the process noise. The step size is the same as readout granularity, noted as λ. ω k can be described as where G is the standard Gaussian distribution with a mean value 0, and the standard deviation θ x,y,z represents the projection of scattering angle θ rms space in three dimensions. Then we can rewrite Eq. 3.1 as Furthermore, for an electron with momentum p, Coulomb scattering distribution follows Gaussian approximation for the Molière's formula [28]: where v is the particle velocity and λ 0 is the step size in the unit of radiation length in the medium. Energy deposition is not considered to simplify our physics model. The measurement equation is defined as follows. Measurement data are anchored by BCs. In between adjacent BCs, the steps m k = (x m , y m , z m ) T k are linearly interpolated hit points spaced by the detector granularity λ. The measurement equation will be presented as where H is the the measurement matrix determined by observation process of the detector. The associated measurement noise is δ k = (G(0, σ), G(0, σ), G(0, σ)) T k , where measurement uncertainty σ includes contributions from detector spatial resolution and uncertainties introduced in the pre-processing. Eq. 3.4 is expanded in the form of a matrix Kalman filter fuses the physical model prediction and measurement to optimize the state vector s k . Based on the least mean-square estimation, a linear weight is determinated to minimize the covariance, defined C k , between s k and its true value. s k and C k are updated in every step k through iteration. The iterative formulae of kalman filter can be derived from Eq. 3.1 and Eq. 3.4 and more details can be found in [26].
Bayesian formalism is introduced in our model of Kalman filter to determine the values of process noise ω k and measurement noise δ k [13,14]. The covariance matrices of two vectors are denoted as Q and R, respectively. We specify two sets of values Q = {Q 1 , Q 2 , ..., Q n Q } and R = {R 1 , R 2 , ..., R n R } representing the samples of possible process noise and measurement noise. Every possible pair [Q i , R j ], where Q i ∈ Q and R j ∈ R, will be selected as the input parameters of Q k and R k for Eq. 3.6. With all the measurements up to step k, M k = {m 1 , m 2 , ..., m k }, the probability P (Q i , R j | M k ) is calculated according to At last, The values θ rms space and σ at the step k are determined by maximizing the above probability Subsequently, we estimate the optimal state vector and covariance matrices at the step k with [Q k , R k ] and are ready for filtering at the step k + 1.
The differences between signal and background tracks are expected to be the largest at the two ends, because of the features of single-electron track and double-electron track. Therefore for both two ends of track, an optimal ratio tuned of the principal track length (20%) from the linearly interpolated hit points is selected for KFB. In Fig. 2c, the dashed red rectangle (L-shaped polygon) marks one end of the background (signal) track used for filtering. The dashed black track represents the true trajectory of events in Geant4 simulation. One can see the red KFB tracks match the MC-truth reasonably well.

Track parameters
Once the event track is reconstructed, we extract six topological parameters for signal identification. When the principal track is determined, the total deposited energy of the principal track E p can be calculated. Both the 0νββ signals and the background events may have the multiple tracks, but background events from γ-rays are more fragmented and hence have smaller E p in general. We extract three parameters from the KFB tracks. With θ rms space determined in KFB, momenta at the ends of the reconstructed track can also be estimated from Eq. 3.3. We defineP as the larger of momenta at two ends. BB and Energy loss per unit travel length dE/dx are calculated more precisely along the reconstructed track. The length of dE/dx and the radius of BB are tuned as 21 mm and 12 mm for the best discrimination power.
We refer the smaller of dE/dx and BB values at two ends of a track as dE dx and E BB respectively. However, compared to E BB , dE dx is more representative of the Bragg peak of ionization energy loss. Besides, other two parameters could be calculated accurately together with the subordinate tracks. The parameter E Space defines energy in a unit volume  around the energy-weighted center of the event and represents the fragmentation of tracks. For pixel (strip) readout, the unit volume is a sphere (circle) with a radius of 57 mm, the size of which is optimized for signal and background discrimination. In addition, the number of tracks N T racks is also calculated and used to further improve event identification. Fig. 3 shows two examples of event tracks to illustrate the definition of E BB , E Space , and N T rakcs . Because of the fragmentation of background event tracks, N T racks and E Space are smaller than these of 0νββ signal in general. The correlation maps of the six parameters for signal and background events under the (1 mm, 3%) configuration are shown in Fig. 4. As expected, dE dx and E BB are highly correlated. The correlation of background events is about 70%, while for signal the value is 44%. Correlations among E p , E space , and N T racks are also observed, because they are all related to the degree of event track dispersion. However, E p , dE dx , andP show almost no correlation, which demonstrates the validity and the leading role of these three parameters for signal identification.
The distributions of E p , dE dx , andP for the (1 mm, 3%) configuration are shown in Fig. 5. In the E p distribution of background, double escape peaks, Compton backscattering peaks, and a full absorption peak at 2447 keV are visible and marked. The majority of 0νββ signals deposit most of the energy in principal tracks and the E p distribution is concentrated within the ROI. The most effective parameter for the signal and background discrimination is dE dx . The dE dx of signals is mostly in the range 15 to 30 keV/mm, but background less than 10 keV/mm. Background events with dE dx more than 10 keV/mm usually deposit a small amount of energy on the principal track, resulting a small E p . The distribution ofP for signal is mainly concentrated in less than 1 MeV/c region, while background is above 1.5 MeV/c. Because of the wide range of energies in E p , the background events also cover a wide range ofP . In addition, an excess of background events around 1 MeV/c is due to mis-identification of the principal track. When subordinate tracks are too close to the principal one, DBSCAN may group them together and result in smallerP values from KFB. Signal events with dE dx less than 15 keV/mm orP over 1 MeV/c are mainly the ones that two electrons share extremely lopsided partition of the Q-value and resemble a background event physically.
The distributions of E BB , E Space , and N T racks for the (1 mm, 3%) configuration are shown in Fig. 6. Compared with dE dx in Fig. 5, the distribution of E BB is similar but the separation between signal and background is not as distinctive. The E space distribution of 0νββ signals is concentrated within the ROI while these of background events are near zero. From the N T racks distribution one can see that about half of 0νββ signal has a single track since emitted electrons lose energy via continuous scattering. The majority of background events have three or four tracks because γ background events deposit energy mainly by Compton scattering at multiple sites.

Classification significance
Three parameters from principal track, E p , dE dx , andP are used as rectangular cuts and the results are shown in Table 1 for the five configurations. We define the discrimination significance Ξ = s √ b , where s is the signal efficiency and b the background efficiency after the cuts. In the (1 mm, 3%) configuration, the background rate of 232 Th is suppressed by more than three orders of magnitude while keeping s at about 30%, and the corresponding Ξ is 12.7. The discrimination for 238 U chain is less significant because the 2447 keV full absorption peak from 214 Bi overlaps with signals in E p . In the (1 mm, 6%) case, Ξ decreases evidently because the Compton backscattering peak of 208 Tl contaminates the wider ROI, similar to what 2447 keV peak does for all configurations. Compare to the high granularity case with the same energy resolution, Ξ decreases only slightly in the (3 mm, 3%) configuration because the 1-mm granularity is an over-kill for diffused tracks in a large TPC. Under the (3 mm, 1%) configuration, Ξ is up to 13.8  Figure 5: Distributions of the three parameters from the principal track though the KFB workflow, including E p , dE dx , andP , for the (1 mm, 3%) configuration. E p represents the total deposited energy of the principal track. In the E p distribution of background, double escape peaks, Compton backscattering peaks, and a full absorption peak at 2447 keV are visible and marked. dE dx represents the energy loss per unit travel length. The statistical length is 21 mm.P represents the momenta at the ends of the reconstructed track. The Y axis is normalized by the highest value in each plot.
better determination of E p . For the strip-readout scheme, the discrimination significance Ξ of 232 Th and 238 U is reduced to 7.3 and 4.6 respectively.
The Toolkit for MultiVariate Analysis in ROOT [29]. BDT classifies signal and background with complicated contours in a multi-dimensional parameter space established by training data. Compared with the rectangular cuts, the improvement of Ξ is smaller than 5% for all configurations with pixel readout, which demonstrates the orthogonality of the three parameters. For the strip readout configuration, the tracks are reconstructed in two two-  Table 1: Effects of rectangular cuts and BDT cuts on s , b , and Ξ in the five Configurations. Background events from 232 Th and 238 U are considered separately. We define the discrimination significance Ξ = s √ b , where s is the signal efficiency and b the background efficiency after the cuts. dimensional planes independently and two sets of dE dx andP are obtained. The BDT cuts get more effective and increase Ξ by approximately 20%.
When E Space , N T rack , and E BB are added as input for BDT as well, we observe additional 20% improvement in discrimination power approximately for all configurations. The best Ξ we achieve is 17.4 (10.3) for 232 Th ( 238 U) in the (3 mm, 1%) configuration. Table 1 shows the results of BDT cuts with all the input parameters.

0νββ search sensitivity study
To illustrate the effectiveness of event discrimination, we calculate the background levels and 0νββ search sensitivity of the PandaX-III experiment with and without topological cuts. 232 Th and 238 U sources are considered. Without any topological cuts, the total count of background level is 152 Count per Year (CPY) in the ROI, which is defined as 0.85 × FWHM (i.e., twice the standard deviation of a Gaussian peak) around the Q-value. The 232 Th chain contributes to 74% of the background events and the 238 U chain 26%. The majority of the background are from the acrylic field cage, the copper liner, and the stainless vessel. Table 2 shows the background CPY in different geometry parts.
The BDT is re-trained with mixed 232 Th and 238 U backgrounds and cut criteria readjusted to maximize the search sensitivity of 0νββ. Fig. 8 shows the signal and background spectra without cuts, with BB cut, and with BDT cuts respectively. We can see the  Table 2: Within ROI, the background rates from the 232 Th and 238 U decay chains in different geometry parts for the configuration of (3 mm strip, 3%). The majority of the background are from the acrylic field cage, the copper liner, and the stainless vessel.  background count after BDT cuts goes down by an order of magnitude, compared with that after BB cut. Fig. 7 shows the classification results of BDT cuts. The cut line for the best Ξ is shown as the green dash line. Within ROI, the optimized BDT cuts suppress background by 3.2 × 10 −3 , resulting a background of 0.49 CPY. The corresponding signal efficiency of topological cuts is 50% and the overall efficiency is 35%. Compared with the design target [10], the background rate is 10.4 times lower, thanks to the BDT cuts based on KFB tracks. The background rate is 3.1 times smaller than the one with the original feature cuts developed in Ref. [12] and meanwhile the signal efficiency is 1.5 times higher. Assuming null results, the exclusion limits are established following a method outlined in Ref. [33]  Counts / 5 keV Background wo cut wi BB cut wi BDT cuts Figure 8: Effects of BB cut and BDT cuts for the signal (left) and background (right). The background count after BDT cuts obviously goes down by about an order of magnitude, compared with that after BB cut. In the BDT cuts, Events with energy less than 2410 keV are rejected because of the E p influence.

Conclusion and discussion
In summary, we present a Kalman Filter based reconstruction workflow for meandering tracks of MeV-scale electrons traveling in a high pressure gas medium. With Monte Carlo simulated data, we apply the technique to 0νββ searches, and demonstrate a significant improvement in signal and background discrimination. Five different configurations of readout schemes and energy resolutions for a gaseous TPC are compared in terms of discrimination power. Impact of the spatial granularity of 1 mm and 3 mm is marginal, partially because track width due to electron diffusion is nominally larger than detector granularity. On the other hand, deterioration of energy resolution can negatively affect the identification of background events from 232 Th and lower the discrimination power. For the best configuration of 3 mm granularity and 1% resolution, cuts based on the BDT classifier can suppress 232 Th ( 238 U) background by a factor of 8.3 × 10 −4 (1.5 × 10 −3 ), while keeping 50% (40%) of the signals. In our simulated detector, a background rate of 0.11 CPY is achieved with a signal efficiency of 40%. Equivalently less than one background event is expected in 5 years, pushing the search towards background-free regime.
With a strip-readout scheme as the one currently used in PandaX-III, we can still suppress the background by a factor of more than 300 while keeping at least 50% of the signal with topological cuts from our improved track reconstruction. After the traditional calorimeter cuts and BDT cuts, the background rate is 0.49 CPY for PandaX-III, and the corresponding half-life sensitivity to 0νββ is 2.7 × 10 26 yr (90% CL) for an exposure of five years. Compared to the previously published topological cuts [12], our new approach improves the search sensitivity by a factor of 2.4.
Further improvement of the workflow can be expected by Extended Kalman filter (EKF) [22]. EKF may take into account of ionization energy loss together with Coulomb scattering processes and estimate the momenta of the particle more accurately. One can reconstruct the entire tracks with momentum information along the track, which would be a powerful background suppression tool for background-free searches of 136 Xe 0νββ.
Broader application of the reconstruction of meandering track with KFB is warranted.
We could estimate the measurement noise in KFB to study the diffusion effect of the tracks, which can help determine the position in drift direction for the PandaX-III experiment. The track reconstruction can also possibly locate the 0νββ vertex, which facilitates the tagging of barium ions in a gaseous TPC [34]. In addition, EKF can be used to determine the energy and angular distributions of the electrons emitted from 2νββ and 0νββ processes to help understand the decay mechanisms. We also envision possible application of the KFB track reconstruction in determining directions of meandering charged particles tracks in general. For example, the direction of MeV-GeV scale Compton or pair-produced electron tracks of γ-ray imaging telescopes [35,36] can be better reconstructed with our approach to improve the angular resolution.