Suppression of Cosmic Muon Spallation Backgrounds in Liquid Scintillator Detectors Using Convolutional Neural Networks

TL;DR

This paper presents a convolutional neural network algorithm that effectively identifies $^{10}$C backgrounds in liquid scintillator detectors, significantly reducing cosmic muon spallation backgrounds for neutrinoless double-beta decay searches.

Contribution

The authors develop a CNN-based method that uses light emission correlations to identify $^{10}$C backgrounds, improving background rejection in liquid scintillator experiments.

Findings

Identifies 61.6% of $^{10}$C at 90% signal acceptance in a KamLAND-like detector.

Could identify 98.2% of $^{10}$C with perfect light collection.

Algorithm is independent of vertex and energy reconstruction, complementing existing methods.

Abstract

Cosmic muon spallation backgrounds are ubiquitous in low-background experiments. For liquid scintillator-based experiments searching for neutrinoless double-beta decay, the spallation product $^{10}$C is an important background in the region of interest between 2-3 MeV and determines the depth requirement for the experiment. We have developed an algorithm based on a convolutional neural network that uses the temporal and spatial correlations in light emissions to identify $^{10}$C background events. With a typical kiloton-scale detector configuration like the KamLAND detector, we find that the algorithm is capable of identifying 61.6% of the $^{10}$C at 90% signal acceptance. A detector with perfect light collection could identify 98.2% at 90% signal acceptance. The algorithm is independent of vertex and energy reconstruction, so it is complementary to current methods and can be expanded to other background sources.

Full text

††journal: NIM Section A

1 Introduction

Neutrinoless double-beta decay ($0\nu\beta\beta$) is a hypothetical decay process by which a nucleus ejects two electrons and no neutrinos, therefore violating lepton number by two units. The observation of this process would demonstrate that neutrinos are Majorana fermions, yield valuable insight into the mechanism behind neutrino mass generation, and support a theoretical framework for matter-antimatter asymmetry in the early universe. Experiments currently being planned will instrument up to a few tons of isotope aiming for $0\nu\beta\beta$ half-life sensitivities between $10^{27}-10^{28}$ years [1, 2, 3].

Large liquid scintillator detectors are attractive for these searches because they offer cost effective scaling to large volumes and effective background reduction through self-shielding, spatial and temporal coincidence analyses, and pulse shape discrimination. The KamLAND-Zen experiment has shown the effectiveness of this technique by setting the most stringent limit on the $0\nu\beta\beta$ half-life independent of isotope $T^{0\nu}_{1/2}>1.07\times 10^{26}$ for 136Xe [4].

Due to the size of these detectors, the typical 20% photocathode coverage, and $\gtrsim$1 ns timing resolution, all energy depositions are assumed to originate from a single point in space and in time. However, we know that different particle species have characteristic topologies which change in time. Even with current detector performance, these could be used for particle identification and background reduction in a variety of analyses. Recent advancements in photodetector technology, with the advent of Large Area Picosecond Photodetectors (LAPPDs) [5, 6], decrease the photon arrival time uncertainty to $\leq 100$ ps. This allows for particle identification by topology [7] and may also permit the reconstruction of particle direction through the separation directional Cherenkov light from the abundant isotropic scintillation light [8]. This has been demonstrated with muons [9, 10, 11] and recently with 90Sr $\beta$-decays [12].

One of the primary cost drivers of $0\nu\beta\beta$ experiments is the depth at which they must be located in order to minimize backgrounds due to cosmic muon spallation. In liquid scintillator detectors doped with $0\nu\beta\beta$ isotopes with endpoints below 3 MeV, such as 136Xe ($Q=2.458$ MeV), the critical long-lived light isotope is 10C ($Q=3.648$ MeV, $\tau_{1/2}=19.29$ s). In the current KamLAND-Zen result, the largest background is the two-neutrino double-beta decay ($2\nu\beta\beta$) of 136Xe, which is only reducible with improved energy resolution. In the absence of scintillator or photodetector upgrades, the largest reducible background is 10C [4]. Looking at future multi-ton scintillator experiments, the increased size and shallower depth make 10C the largest background for a $0\nu\beta\beta$ search with the JUNO experiment [13, 14]. The 10C background will also inform the choice of the depth and location of the proposed THEIA experiment [15].

The production of neutrons and light isotopes in muon spallation is an active area of study [16, 17, 18, 19]. 10C is a relatively common spallation product, predominantly made through $(\pi^{+},np)$ [16]. Results from the Borexino experiment suggest a three-fold-coincidence of muon, neutron, and 11C decay can be used to tag the 11C decay relative to the neutron capture vertex and muon track [20, 21]. This three-fold-coincidence, accompanied by a neutron in the final state, is also used to tag 10C  in the KamLAND-Zen analysis [4]. However, muon spallation is inherently chaotic and the copious number of neutrons produced through this process lead to frequent periods of detector instability. This is especially true for the highest energy events, like showering muons, which produce most of the light isotopes. The problem can be addressed with improved electronics, but due to the large dynamic range between muon events (volt-level signals) and neutron events (millivolt-level signals), difficulties remain. The problem is further amplified by theoretical uncertainty in the fraction of 10C  production with neutron final states, and the spatial distribution relative to the muon track. For these reasons, an independent method of identifying 10C is useful for both reducing the background and verifying our understanding of the spallation process.

In this work we demonstrate that an algorithm based on a ten-layer convolutional neural network (CNN) developed for machine vision applications can effectively separate 10C from $0\nu\beta\beta$ events in a kilo-ton scale liquid scintillator detector, like the current KamLAND detector, without relying on muon or neutron coincidences. Since the technique is spatially invariant, it goes a step farther and is also independent of vertex reconstruction. We then perform a series of studies to understand what information the CNN is using in its discrimination and how the performance of the algorithm changes as we improve the KamLAND-like detector’s performance.

The paper is organized as follows. Section 2 describes the topological differences between $0\nu\beta\beta$-decay and 10C events that allow them to be distinguished by the algorithm. Section 3 provides the details of the detector Monte Carlo (MC) simulation and Section 4 provides the details of the algorithm. These are followed by the results and conclusions in Sections 5 and 6.

2 Topology of $0\nu\beta\beta$-decay and 10C Events

In a liquid scintillator detector, charged particles deposit energy which excites organic molecules. These molecules subsequently de-excite by releasing photons which are detected by single photon detectors like photomultiplier tubes (PMTs). The intrinsic timing of these processes is on the order of ns. Neutrino detectors with large volumes (diameters $>$10 m) are typically instrumented with large PMTs with the capability to resolve photon arrival times on the order of 1-5 ns.

Given the similarity between the time scales of the scintillation process and PMT readout, energy deposits are assumed to be point-like when reconstructing energies and vertices of physics events. This neglects two effects. Gamma-rays at these energies scatter multiple times with a mean free path on the order 10 cm which leads to a smearing of the vertex in time and space. By comparison, electrons travel $\sim$1 mm at these energies, but they are above Cherenkov threshold. Therefore, the electrons produce some directional light which is not absorbed by the scintillation process. This information is encoded in the pattern of photons arriving at the PMTs and can be used to identify different categories of events.

In a large fraction of $0\nu\beta\beta$ events, the electrons exit the nucleus roughly back-to-back evenly dividing the available energy [22]. For 136Xe, this leads to two electrons each with a kinetic energy of roughly 1.23 MeV. The Cherenkov threshold for electrons in this liquid scintillator is 0.16 MeV. Examining the path of these electrons in Monte Carlo, we find that they travel $7.1\pm 0.9$  mm with a total distance from the origin of $5.6\pm 1.0$ mm in $26\pm$4 ps and drop below Cherenkov threshold after $24\pm 3$ ps. We also find that the final direction of the electron before it stops does not match the initial direction, however, the total scattering angle is small while Cherenkov light is emitted.

The decay of 10C is more complicated. 10C  is a $\beta^{+}$ decay that proceeds through one of two excited states as shown in Fig. 1. The event observed by the detector is a combination of a positron, the two subsequent 511 keV annihilation gammas, and either one gamma with energy of 718 keV (98.5%) or two gammas with energies of 718 keV and 1021 keV (1.5%). In addition, the 718 keV excited state is long-lived with a lifetime of 1 ns. This is significant on the time scale of events in liquid scintillator detectors.

The defining feature of the $0\nu\beta\beta$ is the two electrons above Cherenkov threshold and the ring-like features that algorithms, either traditional algorithms based on spherical harmonics or those based on machine learning, should pick out. Since the number of Cherenkov photons is low, the positron from 10C could mimic some of these features. However, the 10C vertex is significantly smeared in time and space by the three or more gammas in the final state. In the region of interest for 136Xe, the positron has an energy of $\sim$0.7 MeV. Examining the path of the positron in Monte Carlo, we find that they travel $3.6\pm 0.6$ mm for a total offset from the origin of $2.5\pm 0.6$ mm in $13\pm 2$ ps. It takes $11\pm 2$ ps for the positron to fall below Cherenkov threshold.

3 Detector Simulation

These studies are performed using a Geant4-based Monte Carlo (MC) with a simplified 6.5 m radius sphere of liquid scintillator. The properties of the liquid scintillator are chosen to match those of KamLAND. We record all photons that reach the outer surface and apply post processing to account for a variety of light collection configurations. This simulation is the same as those used in Ref. [7] and Ref. [8] where more details can be found on the modeling of the LS.

Two types of events are generated for this study: $0\nu\beta\beta$ decay of 136Xe and the $\beta^{+}$ decay of 10C background. The kinematics of $0\nu\beta\beta$ decay events are simulated using a custom MC event generator with momentum and angle-dependent phase space factors from [22]. 10C events are simulated using the default isotope decay generator in GEANT4 [24, 25]. This correctly accounts for the long-lived first excited state of 10B, but does not include the formation of positronium.

The energy spectrum of both event types is shown in Fig. 2 in terms of detected photo-electrons for a detector with perfect light collection. Since 10C event has a much broader spectrum comparing to 136Xe signal, an energy cut is placed at 2.2 MeV to 2.7 MeV prior to training in order to remove background events outside region of interest.

The simulation of these events are distributed either at the exact center of the detector, or uniformly within a 3-meter diameter spherical volume located at the center of the detector111Unless otherwise specified, the two types of event distributions will be referred as center and 3m Sphere in this text.. The latter matches the dimensions of the KamLAND-Zen mini-balloon, which contains the 136Xe-doped scintillator. Fig. 3 compares PE arrival times between $0\nu\beta\beta$-decay and 10C events. The time smearing has been included but the PMT quantum efficiency is added in later steps. In the events at the center, the excess early light from $0\nu\beta\beta$ and late light from 10C is evident. This pattern is repeated in the 3m sphere events but the shift in arrival time due to the vertex positions washes out some of the features. We note that the inclusion of positronium in the model would further increase the late light contribution in 10C decays, improving the results of this study.

3.1 Clock Latching Algorithm

The photon arrival time is an important parameter of the simulation, so we model the digitization of the signal using a clock latching algorithm. Before latching, the hit time of each photon is smeared by a Gaussian Probability Density Function with a 1 ns standard deviation in order to model the performance of these large PMTs. The sample time of the KamLAND electronics is 1.5 ns. When the first unsmeared photon reaches the 6.5 m boundary of the simulation, a clock will start ticking at an 1.5 ns interval. Any photon that arrives after the current clock tick but prior to next tick will be latched to this tick. A total of 30 ticks, 45 ns after the first photon, will be used to record the input event. Furthermore, 4 additional channels are added prior to zero time in order to take into account the backward time smearing of photon. A 2D image of the detector is then formed for each tick. We call each of these images a channel. Each event contains 34 channels, 33 from -6 to 45 ns in 1.5 ns increments and one larger bin for the remaining late photons. Fig. 4 contains four example channels for $0\nu\beta\beta$ and 10C. The distinctive double Cherenkov rings from $0\nu\beta\beta$ decay are visible in the first channel. The higher light levels from 10C are seen in the example channels corresponding to later times.

3.2 Gray Disk PMT Model

The photo-coverage of the detector is introduced by circumscribing a circular region around each PMT location, known as a gray disk. We used the KamLAND PMT locations for the 1,325 17-inch PMTs as a representative PMT layout[26]. The total gray disk area was adjusted to yield the desired photo-coverage. When a photon is produced in the MC, it is associated to the closest PMT. If this photon passes through the gray disk region of the associated PMT, one photon hit is recorded for that PMT position; otherwise, the photon is rejected. This method limits the total photocathode coverage to 40%, since at that point the PMTs start overlapping.

3.3 Quantum Efficiency

The PMT quantum efficiency (QE) is wavelength dependent and is another important simulation input. The QE dependence is accounted for during the event generation. In GEANT4, a QE bit associated with each photon is introduced to indicate whether or not it was recorded by the PMT. However, in this study, we need a varying QE to act as a pressure222In this context, the term pressure refers to the level of difficulty for a neural network to classify events. See also Sec. 4.2 parameter in order to demonstrate the neural network’s performance. A stepwise QE cut is introduced to accommodate this requirement. If the QE bit indicates that a photon was detected, then the photon always passes the QE cut. Otherwise, photons are randomly rejected based on the desired QE pressure. In the ideal situation the QE is 100% and all photons will be indiscriminately recorded.

For our baseline KamLAND-like detector model, the current model, we assume a QE of 23% and photocathode coverage of 19.6%. We define an example upgrade scenario where these parameters are roughly doubled to a QE of 36.2% and photocathode coverage of 42%. This is close to what is being proposed for the JUNO and upgraded KamLAND experiments.

4 Event Classification Algorithm

CNNs are a type of deep neural networks commonly used in the field of computer vision. To perform classification tasks using CNNs, a pixelized image is processed through several so-called layers, each containing a linear transformation step, followed by a non-linear activation. Among these layers, convolutional layers play the most important role. As the name implies they involve the application of a convolution operation between two functions. For continuous functions, it represents the Fourier transform of the products. Applying this procedure to a 2D discrete surface gives rise to the convolutional layer.

The convolutional filter, the kernel for the convolutional layer, is a fix-sized grid with specific values assigned to each block. The filter is scanned throughout the image body, and each image pixel is multiplied by a filter weight. Finally, the convolution operation is completed by summing the element-wise multiplications.

The convolutional layer is also capable of taking information from multiple channels. In this case, the convolution operation is conducted separately over each channel, and the output values are summed and fed into the next layer. For example, when a CNN is used to classify photographs, the input contains four channels: Red(R), Green(G), Blue(B), and Grey Scale. For this work, we have 34 time-based channels as described in Section 3.1.

The output of the convolutional layer contains features that are fed into the fully connected layer. This layer is where the high level decisions are made, ending with an image being classified into one of several categories with an assigned probability. Other layers in the network structure include pooling layers and dropout layers. A pooling layer reduces the image dimensionality by extracting only the maximum value from each pooling filter to represent the image. It significantly increases the processing speed with almost no sacrifice of classification accuracy. The dropout layer prevents overfitting by randomly disabling neurons in the hidden layer with a predefined dropout rate [27]. For extensive details on the functionality of CNNs, we refer the reader to Ref. [28].

4.1 Network Design

The CNN used in this work is implemented in Keras [29] with a Tensorflow backend [30]. The general outline of the network is shown in Fig. 5. It can be divided into a convolutional part and a fully connected part.

The convolutional part contains 5 sets of layers where each set includes a convolutional layer, a batch normalization layer, a pooling layer, and a dropout layer. For simplicity, one such set of layers is often simply referred to as a convolutional layer. Each time a pooling layer is introduced, the image reduces to 25% of it’s original size. Therefore, the convolution part is confined to five iterations.

The fully connected part also contains 5 sets of layers. Each set includes a fully connected layer, a batch normalization layer and a dropout layer. Again, this set of layers is often collectively referred to as a fully connected layer. With respect to the overall depth, the fully connected part is not limited by the image size or pooling layers. However, going too deep with the fully connected part will complicate the model and lead to overfitting. Dropout layers are inserted throughout the network in order to prevent this effect while decreasing the processing time [27].

Normalization of data plays an important role in classification tasks. Without normalization, the neural network will reveal some anomalous behavior, including non-convergance and high probability of misclassification. During a pre-training stage, each pixel is scaled to a value below unity. During the network design stage, several different normalization schemes were considered, including vector normalization, channel-wise standard scalar normalization, and batch normalization. For this work, batch normalization was chosen. This means that the normalization is performed for each incoming batch, transforming the input data such that the mean is zero with a standard deviation of one.

Hyperparameters refer to parameters in the network that are predefined before training, and stay constant throughout the training stage. Hyperparameters of neural networks define the structure of the network, and changing these parameters turns the neural network into a different model. Typical hyperparameters include, but are not limited to, the number of layers in the network, the number of nodes in each layers, the size of the filters, and the dropout rate.

A hyperparameter search can result in a significant improvement of performance [31]. For our model, we performed the search with hyperopt [32]. Three hyperparameters are tuned to achieve the best performance, including number of nodes, number of fully connected layers, and dropout rate. We selected a preliminary training data set of 20,000 3m sphere events with the current detector configuration. We then search across a continuous range for the dropout rate between 0 and 1 and several discrete values for the number of fully connected layers and nodes. A random search [31] of 50 attempts is executed to determine the best result. The 50 trials are evaluated and compared for the best background rejection capability. After tuning, the validation accuracy increases from 49.7% to 77.3%.

4.2 Training

The algorithm is trained and validated on MC data sets, generated according to Section 3. We study two different data sets: 70,000 centered events, and 50,000 3m sphere events. All events are stored in an 6 dimensional numpy array. The six dimensions are correspondingly: photocoverage pressure, QE pressure, event index, time channel, polar angle and azimuthal angle. The algorithm is trained and validated independently for each data set and each category contains equal amounts of $0\nu\beta\beta$ signal and 10C background events.

The data sets are separated into training and validation subsets with a 3:1 ratio. The network is trained on batches of 10 events over 30 training cycles. An RMSProp optimizer [30] is used to apply backward propagation optimization based on binary cross entropy [33]. Due to the large data volume, sparse matrix and batch generator technology is applied to reduce memory consumption. During training stage, a learning rate decay scheduler is incorporated to reduce systematic fluctuation. After training, the CNN is applied to the validation data set to study the out-of-sample performance. A bad validation rate in the presence of good training accuracy is indicative of overfitting and we do not see this effect.

While training the CNN on each event category, the photo-coverage and PMT QE were scanned over a wide range of values to better understand the CNN performance under different levels of classification difficulty. These levels of difficulty are referred to as pressure. The photo-coverage was allowed to vary from 20% to 40% and QE from 23% to 56%. A total of 99 CNN models were trained and evaluated for the pressure maps.

5 Results

The well-trained neural network outputs a single floating point number between 0 and 1 for every input event. This number is the sigmoid output, since it comes out of a sigmoid activation function. The sigmoid output serves as the probability or metric for event classification. If the sigmoid output of a given event is close to 1, it means the event is likely to be a signal event. The value of the sigmoid output used for the classification can be more or less stringent depending on the required signal purity versus signal acceptance. Fig. 6 shows the sigmoid output for the current configuration and a possible upgrade scenario.

After the CNN is applied to the validation data set, the value of the sigmoid output cut can be varied to generate a Receiver Operating Characteristic (ROC) curve. For our application, the result is simply the signal acceptance as a function of the background rejection. Fig. 7 shows the ROC curves for the current configuration and the upgrade scenario. We find that at 90% signal acceptance we can reject 61.6% of the 10C. For the scenario with increased coverage and QE this increases to 81.3%. The central events with the standard KamLAND configuration have a 97.7% rejection at 90% signal acceptance.

Fig. 8 shows the pressure maps which scan different QE and photo-coverage configurations for both the central and the 3m sphere events. As expected, the CNN performs best for centrally located events and for higher QE and photo-coverage. The results also indicate that increasing the total light collected, whether by increasing in QE or photocathode coverage, leads to improved performance.

Within the parameters of this study, we find that it is possible to reach $>$99.98% discrimination for central events. This indicates that higher isotope concentrations that lead to more centrally distributed $0\nu\beta\beta$ events are advantageous and this motivates future design studies. We also studied the algorithm with the 3m sphere events and perfect light collection and find 98.2% rejection at 90% signal acceptance.

We use the fluctuations observed Fig. 8 to understand the uncertainty in the algorithm. This is estimated by calculating at each point in the grid the standard deviation relative to the 4 adjacent neighbors and averaging this value over the pressure map. For central events, the uncertainty is 0.16%, while the 3m Sphere events give an 1.9%.

In order to understand what feature the CNN is using to discriminate, we produced a data set where the Cherenkov light was removed. The results for central events with the standard KamLAND detector configuration indicated a 4% decrease in the rejection. The 3m sphere events show a 2% increase in the rejection. This could indicate that Cherenkov light is interfering with the networks interpretation of the rise-time of the scintillation light, however this is also within our estimated uncertainty for the algorithm. As a whole, we find that the Cherenkov signal is not the dominant feature and that the scintillation light topology is driving the event separation. Within the scintillation light topology the photon timing information contains significant discriminating power. Removing spatial information and using a one-dimensional fully connected network with timing as an input we found background rejection drops from 61.6% to 55.0% at 90% signal acceptance.

6 Conclusion

Liquid scintillator detectors have been at the heart of many of the great discoveries in neutrino physics and have been a leading technology in the search for $0\nu\beta\beta$. However, the algorithms used to analyze their data have remained relatively unchanged for decades. In this work, we apply an algorithm from computer vision based on a CNN to extract the fundamentally different physics processes that take place in 10C background events and $0\nu\beta\beta$ signal events. With a standard detector configuration similar to the current KamLAND detector, we find we can reject 61.6% of the 10C background with 90% acceptance of the $0\nu\beta\beta$ signal. A detector with the same geometry and perfect light collection could achieve 98.2% rejection. We also find that the performance can be increased to better than 99.98% for centrally located events. The overall uncertainty of the algorithm is 1.9%.

These results are a basis for future studies combining machine learning techniques based on CNNs with liquid scintillator detectors. In short order, we intend to move to a spherically symmetric CNN [34] and a Bayesian classification that provides a posterior distribution for the classification. In future studies, this algorithm will be applied to other backgrounds with topologies distinct from $0\nu\beta\beta$ decay. These include 214Bi decays on KamLAND 3m balloon and elastic scattering of 8B solar neutrinos. Solar neutrinos are expected to be the dominant background in SNO+ [35]. We are also exploring algorithms which could move beyond simple classification to particle position and direction reconstruction. These studies are benefiting from an abundance of work being done for other applications both inside and outside of particle and nuclear physics and there are many new avenues to explore.

Acknowledgments

This material is based upon work supported by the National Science Foundation under Grant Numbers 1554875 and 1806440. We thank Taritree Wongjirad for contributions to earlier incarnations of this project and Jonathan Ouellet for the 10C decay diagram. We also thank Kazuhiro Terao for useful discussions. This work is done in support of the NuDot, THEIA and of course KamLAND experiments and we thank them for their input. The work at the University of Chicago is supported by U. S. Department of Energy, Office of Science, Offices of High Energy Physics and Nuclear Physics under contracts DE-SC0008172 and DE-SC0015367; the National Science Foundation under grant PHY-1066014; and the Physical Sciences Division of the University of Chicago. We thank Mikhail Hushchyn, Eugeny Toropov, and Ilija Vukotic for advice on application of machine learning techniques to image classification. This research was done, in part, using resources provided by the Open Science Grid [36, 37], which is supported by the National Science Foundation award 1148698, and the U.S. Department of Energy’s Office of Science. This research was also performed with the help of the Shared Computing Cluster at Boston University.