[6,a]A. Simón, aaaCorresponding author. 11institutetext:  Department of
Physics and Astronomy, Iowa State University, Ames, Iowa, USA 22institutetext:
Department of Chemistry and Biochemistry, University of Texas at Arlington,
Arlington, Texas, USA 33institutetext:  Department of Physics, University of
Texas at Arlington, Arlington, Texas, USA 44institutetext:  Institute of
Nanostructures, Nanomodelling and Nanofabrication (i3N), Universidade de
Aveiro, Aveiro, Portugal 55institutetext:  Fermi National Accelerator
Laboratory, Batavia, Illinois, USA 66institutetext:  Unit of Nuclear
Engineering, Faculty of Engineering Sciences, Ben-Gurion University of the
Negev, Beer-Sheva, Israel 77institutetext:  Nuclear Research Center Negev,
Beer-Sheva, Israel 88institutetext:  Lawrence Berkeley National Laboratory,
Berkeley, California, USA 99institutetext:  Ikerbasque (Basque Foundation for
Science), Bilbao, Spain 1010institutetext:  Centro de Investigación en
Ciencias Básicas y Aplicadas, Universidad Antonio Nariño, Bogotá, Colombia
1111institutetext:  Department of Physics, Harvard University, Cambridge,
Massachusetts, USA 1212institutetext:  Laboratorio Subterráneo de Canfranc,
Canfranc-Estación, Spain 1313institutetext:  LIBPhys, Universidade de Coimbra,
Coimbra, Portugal 1414institutetext:  LIP, Departamento de Física,
Universidade de Coimbra, Coimbra, Portugal 1515institutetext:  Department of
Physics and Astronomy, Texas A&M University, College Station, Texas, USA
1616institutetext:  Materials Physics Center (CFM), CSIC & University of the
Basque Country (UPV/EHU), Manuel de Lardizabal 5, 20018 Donostia-San
Sebastián, Spain 1717institutetext:  Donostia International Physics Center
(DIPC), Donostia-San Sebastián, Spain 1818institutetext:  Faculty of
Chemistry, University of the Basque Country (UPV/EHU), Manuel de Lardizabal 3,
20018 Donostia-San Sebastián, Spain 1919institutetext:  Escola Politècnica
Superior, Universitat de Girona, Girona, Spain 2020institutetext:  Argonne
National Laboratory, Lemont, Illinois, USA 2121institutetext:  Departamento de
Física Teórica, Universidad Autónoma de Madrid, Madrid, Spain
2222institutetext:  Instituto de Física Corpuscular (IFIC), CSIC & Universitat
de València, Paterna, Spain 2323institutetext:  Pacific Northwest National
Laboratory, Richland, Washington, USA 2424institutetext:  Instituto Gallego de
Física de Altas Energías, Universidade de Santiago de Compostela, Santiago de
Compostela, Spain 2525institutetext:  Instituto de Instrumentación para Imagen
Molecular (I3M), CSIC & Univ. Politècnica de València, Valencia, Spain
2626institutetext:  Centro de Astropartículas y Física de Altas Energías
(CAPA), Universidad de Zaragoza, Zaragoza, Spain

# Boosting background suppression in the NEXT experiment through Richardson-
Lucy deconvolution

6,7 Y. Ifergan 6 A.B. Redwine 6,b R. Weiss-Babai,bbbNow in Soreq Nuclear
Research Center, Yavneh, Israel. 6,a L. Arazi 20 C. Adams 17 H. Almazán 25 V.
Álvarez 18 B. Aparicio 18 A.I. Aranburu 23 I.J. Arnquist 4 C.D.R Azevedo 20 K.
Bailey 25 F. Ballester 17 J.M. Benlloch-Rodríguez 14 F.I.G.M. Borges 3 N.
Byrnes 22 S. Cárcel 22 J.V. Carrión 26 S. Cebrián 23 E. Church 14 C.A.N. Conde
11 T. Contreras 9,18 F.P. Cossío 2 A.A. Denisenko 24 G. Díaz 22 J. Díaz 14 J.
Escada 25 R. Esteve 6,22 R. Felkai 13 L.M.P. Fernandes 9,17 P. Ferrario 4 A.L.
Ferreira 2 F. Foss 13 E.D.C. Freitas 9,27 Z. Freixa 17 J. Generowicz 8 A.
Goldschmidt 9,17,c J.J. Gómez-Cadenas,cccNEXT Co-spokesperson. 17 R. González
24 D. González-Díaz 11 S. Gosh 11 R. Guenette 10 R.M. Gutiérrez 11 J. Haefner
20 K. Hafidi 1 J. Hauptman 13 C.A.O. Henriques 24 J.A. Hernando Morata 16,17
P. Herrero 25 V. Herrero 11 J. Ho 3 B.J.P. Jones 24 M. Kekic 21 L. Labarga 3
A. Laing 5 P. Lebrun 22 N. López-March 10 M. Losada 13 R.D.P. Mano 11,22 J.
Martín-Albo 22 A. Martínez 17,22 M. Martínez-Vara 6 G. Martínez-Lema 3 A.D.
McDonald 20 Z.-E. Meziani 9,17 F. Monrabal 13 C.M.B. Monteiro 25 F.J. Mora 22
J. Muñoz Vidal 2 C. Newhouse 22 P. Novella 3,c D.R. Nygren 17 E. Oblak 17 M.
Odriozola-Gimeno 22,24 B. Palmeiro 5 A. Para 12 J. Pérez 22 M. Querol 22 J.
Renner 19 L. Ripoll 9,17 I. Rivilla 10 Y. Rodríguez García 25 J. Rodríguez 16
C. Rogero 3 L. Rogers 12,17 B. Romeo 22 C. Romo-Luque 14 F.P. Santos 13 J.M.F.
dos Santos 22 M. Sorel 11 C. Stanford 13 J.M.R. Teixeira 2 P. Thapa 25 J.F.
Toledo 17 J. Torrent 22 A. Usón 4 J.F.C.A. Veloso 2 T.T. Vuong 15 R. Webb 15,d
J.T. White,dddDeceased. 3 K. Woodruff 22 N. Yahlali<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Next-generation neutrinoless double beta decay experiments aim for half-life
sensitivities of $\sim 10^{27}$ yr, requiring suppressing backgrounds to $<1$
count/tonne/yr. For this, any extra background rejection handle, beyond
excellent energy resolution and the use of extremely radiopure materials, is
of utmost importance. The NEXT experiment exploits differences in the spatial
ionization patterns of double beta decay and single-electron events to
discriminate signal from background. While the former display two Bragg peak
dense ionization regions at the opposite ends of the track, the latter
typically have only one such feature. Thus, comparing the energies at the
track extremes provides an additional rejection tool. The unique combination
of the topology-based background discrimination and excellent energy
resolution (1% FWHM at the Q-value of the decay) is the distinguishing feature
of NEXT. Previous studies demonstrated a topological background rejection
factor of $\sim 5$ when reconstructing electron-positron pairs in the 208Tl
1.6 MeV double escape peak (with Compton events as background), recorded in
the NEXT-White demonstrator at the Laboratorio Subterráneo de Canfranc, with
72% signal efficiency. This was recently improved through the use of a deep
convolutional neural network to yield a background rejection factor of $\sim
10$ with 65% signal efficiency. Here, we present a new reconstruction method,
based on the Richardson-Lucy deconvolution algorithm, which allows reversing
the blurring induced by electron diffusion and electroluminescence light
production in the NEXT TPC. The new method yields highly refined 3D images of
reconstructed events, and, as a result, significantly improves the topological
background discrimination. When applied to real-data 1.6 MeV $e^{-}e^{+}$
pairs, it leads to a background rejection factor of 27 at 57% signal
efficiency.

##  1 Introduction

The search for neutrinoless double beta ($0\nu\beta\beta$) decay is the most
promising experimental path to determine whether the neutrino is a Majorana
fermion, with far-reaching implications in particle physics and cosmology
Avignone2008 ; Davidson2008 ; Blennow2010 ; Dell?Oro2016 ; Dolinski2019 .
Presently, several collaborations pursue different technologies for detecting
$0\nu\beta\beta$ decay with the leading experiments focusing on 76Ge GERDA2020
; Majorana2019 ; LEGEND2017 , 136Xe KamLAND-Zen2016 ; Loaded_scintillators2019
; EXO2018 ; nEXO2018 ; DARWIN2020 ; JJ2019 ; PandaX-III , 130Te CUORE2020 ;
SNO+2016 ; Loaded_scintillators2019 , and 100Mo CUPID2019 ; CUPID-Mo_results ;
AMORE2019 . The long half-life of $0\nu\beta\beta$ decay (above 1.8 $\times$
1026 yr in 76Ge GERDA2020  and 1.07 $\times$ 1026 yr in 136Xe KamLAND-Zen2016
) makes its detection extremely difficult, with only a few candidate
$0\nu\beta\beta$ events expected throughout the running life of an experiment,
calling for outstanding background suppression capabilities. The next
generation of $0\nu\beta\beta$ decay experiments will aim at half-life
sensitivities of $\sim 10^{27}$ yr, requiring, in turn, a background level
below $\sim 1$ counts/tonne/yr.

NEXT (Neutrino Experiment with a Xenon TPC) is a staged experimental program
aiming at the detection of $0\nu\beta\beta$ decay in 136Xe, using successive
generations of high-pressure gaseous xenon electroluminescent time projection
chambers (HPXe EL-TPC) Nygren:2009zz . The choice of gaseous rather than
liquid xenon is driven by two considerations: (1) the attainable energy
resolution at the Q-value of the decay, $Q_{\beta\beta}$ (which for 136Xe is
2458 keV), is a factor 3 better in gas than in liquid Bolotnikov:1997 ;
Renner:2019pfe ; EXO2018 ; XENON1T_Eres ; and (2) whereas in liquid xenon
events are point-like, the projected length of ionization tracks in high-
pressure Xe gas is $\sim 10$ cm, allowing for background discrimination based
on the track topology, as first pioneered by the Gotthard experiment
Gotthard1998 . The distinguishing topological feature of a two-electron double
beta decay event in high-pressure gas is the appearance of two blob-like
Bragg-peak energy depositions at the opposite ends of the track, in contrast
with single-electron background events which have only one such feature
(figure 1). Thus, by reconstructing the track in 3D and comparing the energy
contained in small regions at its extremities, one can effectively classify
the event as either signal or background.

The current stage of the NEXT program is the NEXT-White111Named after Prof.
James White, our late mentor and friend. demonstrator, whose TPC has an active
volume half a meter in diameter and length Monrabal:2018xlr . NEXT-White
(“NEW”) is a radiopure detector, operated underground under low-background
conditions at the Laboratorio Subterráneo de Canfranc (LSC), using xenon
enriched to $90\%$ 136Xe. Its purpose is to validate all aspects of the
technology on a large scale, including full characterization of the background
model and the technique’s background rejection power, and to demonstrate its
performance on two-neutrino double beta ($2\nu\beta\beta$) decay events. NEXT-
White, which has been running continuously since October 2017, will be
superseded in 2022 by the twice-larger NEXT-100 detector, which will deploy 97
kg of enriched xenon and demonstrate sensitivity to $0\nu\beta\beta$ decay
half-lives on the scale of $10^{26}$ yr JJ2019 . This will pave the way to a
tonne-scale experiment which will be sensitive to half-lives longer than
$10^{27}$ yr NEXT_tonne . Importantly, the NEXT Collaboration pursues in
parallel an extensive R&D program to develop the capability of detecting the
136Ba daughter resulting from 136Xe double beta decays inside a running TPC
using single molecule fluorescence imaging Jones2016 ; McDonnald_PRL ;
Byrnes_2019 ; Thapa2019 ; Rivilla2020 . If successful, these efforts could
boost the sensitivity of HPXe EL-TPCs to half-lives on the scale of $10^{28}$
yr.

NEXT excellent energy resolution – demonstrated to be $\sim 1\%$ full-width
half maximum (FWHM) at $Q_{\beta\beta}$ in NEXT-White Renner:2019pfe  – allows
to reduce background from natural radioactivity and $2\nu\beta\beta$ events in
the $Q_{\beta\beta}$ region-of-interest (ROI) by about 5 and 12 orders of
magnitude, respectively, making the latter completely negligible
NEXT100_sensitivity . NEXT also utilizes track multiplicity, i.e., the number
of distinct tracks in a given event, to reduce background from $\gamma$-ray
interactions or from single beta decays with energy in the $Q_{\beta\beta}$
ROI (from cosmogenically-produced 137Xe). This effectively removes multiple
Compton scatters, as well as fast-electron tracks with accompanying vertices
of bremsstrahlung and characteristic x-ray interactions, resulting in an
additional background rejection factor of $\sim 10$ while retaining $\sim
70\%$ of the signal NEXT100_sensitivity ; NEXT_tonne . Finally, using the
topology of the remaining single tracks, background is further suppressed by
measuring the energy deposition at the track ends as described above.

Figure 1: Simulated GEANT4 signal and background tracks near $Q_{\beta\beta}$
at 15 bar, from Ferrario:2015kta .

Tests of the topological background rejection in NEXT-White rely on the use of
a 228Th calibration source to produce 2615 keV $\gamma$-rays by the decay of
208Tl. Interactions of these $\gamma$-rays inside the gas lead to the
production of electron-positron pairs with a total kinetic energy of 1593 keV.
These leave a trace with two Bragg-peak energy depositions at the track
endpoints, mimicking double beta decay events. Previous topological analysis
of such events in NEXT-White (with Compton electrons serving as background)
yielded a background rejection factor of $\sim 5$, while retaining 72% of the
signal events Ferrario:2019kwg  for optimal performance. Extended by Monte
Carlo to $Q_{\beta\beta}$, this analysis provided a background rejection
factor 7.4 with similar signal efficiency. (Note that these values correspond
to the performance of the topological analysis alone, after the preceding
energy and single-track cuts.) Recently, these results were improved by using
a deep convolutional neural network yielding a background rejection factor of
$\sim 10$ with a signal acceptance of $\sim$65% Kekic:2020cne .

In this work we describe an improved methodology for track reconstruction. The
underlying idea is to enhance the sharpness of reconstructed tracks – degraded
by electron diffusion and spread of light produced in the electroluminescence
(EL) process – through the application of an image deblurring procedure,
namely the Richardson-Lucy (RL) deconvolution algorithm Richardson:72 ;
Lucy:1974yx . We begin by describing signal production and event
reconstruction in NEXT-White, with particular emphasis on previous work
regarding the topological analysis. This is followed by a discussion of the
blurring effects of diffusion and EL light production, and their
quantification through a spatially-dependent point spread function (PSF). We
then describe the implementation of the RL algorithm, which is outlined in
appendix A, within the experiment’s data processing chain, and validate it by
demonstrating the accurate reconstruction of point-like events (individual and
pairs) and muon tracks. The procedure is subsequently employed on 1.6 MeV
$e^{-}e^{+}$ pairs and gamma-induced background events in both Monte Carlo and
data recorded in NEXT-White, yielding a major improvement in topological
background rejection. We discuss the effect of the key parameters of the
deconvolution procedure on its performance, and present the optimized results.

##  2 NEXT-White: event reconstruction and prior work on topological analysis

The NEXT-White TPC Monrabal:2018xlr  consists of a cathode and gate grids
defining a 53 cm-long drift region, a transparent anode plate positioned 6 mm
behind the gate, and a field cage with an inner diameter of 45 cm. An array of
12 Hamamatsu R11410-10 3” photomultiplier tubes (PMTs), constituting the
energy plane, is located 13 cm behind the cathode. At the opposite end of the
TPC, 2 mm behind the anode plate, an array of 1792 SensL series-C 1 mm2
silicon photomultipliers (SiPMs) distributed at a pitch of 10 mm, serves as
the tracking plane. The entire tracking plane area is covered by
polytetrafluoroethylene (PTFE) to reflect light towards the PMTs, with holes
for the SiPMs. The interior wall of the field cage comprises a PTFE tube
coated with a thin wavelength-shifting layer of tetraphenyl butadiene (TPB) to
further optimize light collection. The anode is a 3 mm-thick fused silica
plate, coated on both faces with transparent resistive and conductive layers,
which are themselves coated by TPB. The detector is presently operated under
voltages which define a uniform drift field of 0.42 kV/cm between the cathode
and gate, and a nominal electroluminescence (EL) field of 12.8 kV/cm between
the gate and anode (the EL gap).

When an event occurs inside the sensitive volume, the associated primary
charged particles form Xe excimers and electron-ion pairs along their track.
De-excitation of the former produces prompt vacuum ultra violet (VUV)
scintillation light (“S1”) centered at 172 nm, lasting a few hundred ns. This
light, which is shifted to $\sim$430 nm by the TPB coating the inner surfaces
of the TPC and recorded by the PMTs with a sampling time of 25 ns, provides
the start time $t_{0}$ of the event. The drift field prevents electron-ion
recombination and drives the electrons at a uniform velocity of $0.91$
mm/$\mu$s (for $E_{drift}=0.42$ kV/cm) towards the EL gap. As they cross it
(in $\sim 1.6\;\mu$s), they produce EL light (“S2”), which is also centered at
172 nm and wavelength-shifted by TPB. For a nominal EL field of 12.8 kV/cm and
6 mm gap, this provides, at 10 bar, a light yield of $\sim 450$ VUV photons
per electron crossing the gap Freitas2010 . Depending on the track length, the
duration of the S2 signal varies from a few $\mu$s to a few hundred $\mu$s.
The S2 light is recorded by both the PMTs and SiPMs, with the latter
integrated in 1 $\mu$s time slices.

During operation, NEXT-White is continuously calibrated using 83mKr, which is
introduced into the gas system by an in-line 83Rb source. 83mKr spreads
uniformly throughout the TPC volume and produces point-like energy depositions
of 41.6 keV. These events serve for precise determination of the electron
lifetime and $xy$ S2 response of the detector, to generate 3D correction maps
which are then employed for accurate energy measurements of signal and
background events Martinez-Lema:2018ibw . Furthermore, as detailed below,
83mKr can also be used to obtain the point-spread function that describes both
electron diffusion and the optical response of the tracking plane, a vital
ingredient for image deconvolution.

When processing events with energy higher than 400 keV, SiPM waveforms are re-
binned into 2 $\mu$s, keeping only sensors recording more than 5
photoelectrons (PEs) within each slice. At this stage a second, optional,
higher charge threshold is applied. For each SiPM passing this threshold in
each slice, a 3D “hit” is generated, whose $xy$ coordinates are those of the
associated SiPM, the $z$ coordinate is the product of the electron drift
velocity and time difference between the slice and S1, and the magnitude is
the number of detected PEs. The energy measured by the PMTs in the same time
slice is divided among the reconstructed hits, proportionally to the charge of
their respective SiPMs. Afterwards the energy assigned to each hit is
multiplied by the correction factors derived from the 83mKr 3D maps. Figure 2
(left) shows an example of a 3D hit map of a reconstructed track with color
representing the number of PEs detected.

The resulting hits are grouped in cubic voxels with an event-dependent fixed
size. These are further grouped into connected tracks, using the breadth-first
search (BFS) algorithm Cormen2001_intro_algorithms ; Ferrario:2015kta . Voxels
are considered to be part of the same track if they have a common face, side
or corner. An event may contain more than one track (for example, in the case
of multiple Compton scatters). After building the tracks, the BFS algorithm
further identifies the end-point voxels of each track. These are defined as
the pair of voxels with the longest distance between them, where the distance
between any pair of voxels is defined as the shortest path along the track
that connects them. The SiPM hits contained in the end-point voxels are used
to define the voxel center-of-gravity (COG). Finally, two spheres of a fixed
radius are defined around the end-point voxels COGs, and the energy contained
in them is summed, designating the sphere containing more energy as blob1 and
the one carrying less energy as blob2. The final output for analysis consists
of a collection of tracks and their “blobs”.

To investigate the effectiveness of topological event classification in NEXT-
White, the experimental work described in Ferrario:2019kwg , to which we refer
below as the classical analysis, focused on $e^{-}e^{+}$ pairs produced in Xe
by the 2615 keV gamma of 208Tl, using an external source of 228Th (parent of
the decay chain containing 208Tl). The topological analysis of such pairs
provides an excellent “training arena” for that of $0\nu\beta\beta$ events.
Similar to the track structure of $0\nu\beta\beta$ decay, the pair electron
and positron tracks start from a common vertex, and both end with a Bragg peak
of dense ionization. Such pairs have a total kinetic energy of 1593 keV, and
their selection is based on the identification of the 208Tl double escape peak
in the energy spectrum. At 10 bar the combined continuous-slowing-down
approximation range of such pairs is $\sim 20$ cm, the same as that of
$0\nu\beta\beta$ events at 15 bar (the planned operation pressure of
NEXT-100); the typical projected length of the full event in both cases is
$\lesssim$10 cm. For the sake of the analysis, $e^{-}e^{+}$ pairs are
therefore considered as signal events, while single-electron events (from
Compton scatters) in the region of the double escape peak, are considered as
background.

For the classical analysis the SiPM charge cut was set at a high value of 30
PEs, SiPM hits were grouped in cubic voxels with an event-dependent fixed size
between 10 and 15 mm (figure 2, right), and the blobs were defined by 21 mm-
radius spheres. All choices are revisited in the present work. In the analysis
(classical and the new method presented here), events are required to be fully
contained within the fiducial volume (with all hits at least 2 cm away from
any of the drift volume surfaces) and consist of only one track.
Classification to signal or background is done by comparing the energy
contained in blob2 to a fixed threshold. If blob2 energy exceeds it the event
is classified as signal, and if not – as background.

To quantify the effectiveness of this topological cut on double escape peak
events, one defines the signal efficiency $\epsilon$ as the fraction of signal
events ($e^{-}e^{+}$ pairs) passing the cut, and the background acceptance $b$
as the fraction of background events that survive it. Both parameters depend
on the choice of the energy threshold for blob2. As discussed in
NEXT100_sensitivity , the sensitivity of the experiment to $0\nu\beta\beta$
decays is proportional to a topological figure of merit, defined as
$f.o.m.=\epsilon/\sqrt{b}$, where the signal efficiency and background
acceptance are calculated for events in the $Q_{\beta\beta}$ ROI. Considering
the similarity between 1.6 MeV $e^{-}e^{+}$ pairs at 10 bar and
$0\nu\beta\beta$ events at 15 bar, optimizing the $f.o.m.$ for $e^{-}e^{+}$
events – which are accessible experimentally at large numbers – can be
regarded as a proxy for optimizing the sensitivity of the experiment to
$0\nu\beta\beta$ decays. For a given configuration of reconstruction
parameters, the maximal $f.o.m.$ is obtained for an optimal choice of blob2
energy threshold. As shown in Ferrario:2019kwg  for the optimal figure of
merit, the classical analysis, when applied to experimental data, yields a
signal efficiency of 71.6% with a background acceptance of 20.6% (with small
statistical and systematic errors), i.e., a background rejection factor of
4.9. Similar outcomes were obtained in Monte Carlo (MC) simulations of double-
escape peak events, performed with the NEXUS Geant4-based NEXT simulation
framework Justo_thesis . When evaluating the performance of the classical
analysis at $Q_{\beta\beta}$ using MC events, the background acceptance was
further reduced down to 13.6% (background rejection factor of 7.4) while
keeping virtually the same signal efficiency Ferrario:2019kwg .

As mentioned earlier, a more recent work Kekic:2020cne  improved on these
results by using a deep convolutional neural network to classify 1.6 MeV
double-escape peak events to signal and background, starting from the
voxelized tracks, instead of using the BFS algorithm and blob-based analysis.
This method yielded a background acceptance level of 10% at 65% signal
efficiency. The reconstruction method described in the following sections is
closer in nature to the classical analysis, which therefore serves as the
reference for comparison. Future works will investigate the benefit of
combining image deconvolution and neural network-based analysis.

Figure 2: An example for a reconstructed event in NEXT-White. Left: SiPM hits,
with color representing detected photons; right: same event voxelized in
$15\times 15\times 15$ mm3 bins (from Ferrario:2019kwg ).

##  3 Image blurring effects in NEXT-White

The classical reconstruction using the SiPM hits and voxels is affected by two
blurring mechanisms which degrade the quality of the reconstructed track:
electron diffusion and the spread of EL light on the tracking plane.

The initial track structure is a thin trail of ionization, as shown in figure
1. As the ionization electrons drift towards the gate, elastic collisions with
Xe atoms lead to transverse and longitudinal diffusive spread of the charge
cloud around the track “backbone”. The root-mean-square (r.m.s.) diffusive
spread of each point-like element of the initial track is proportional to the
square root of the drift time, and therefore – since the drift velocity is
constant – to the square root of the distance between this element and the
gate (i.e., its $z$ coordinate). Under the operating conditions of NEXT-White,
this effect can be on the cm scale.

The second contribution to the overall blurring occurs as the electrons cross
the EL gap, where they emit light isotropically in $4\pi$. VUV photons emitted
towards the anode are absorbed in TPB and re-emitted (again isotropically) as
blue photons, which are subsequently transmitted through or reflected on the
interfaces of the multi-layered anode plate, with additional reflections from
the tracking plane PTFE cover. VUV EL photons emitted towards the gate and
absorbed on it may result in secondary photoelectron emission, creating an
additional discrete “halo” of diffuse light around the event. All of these
processes combine to further optically smear the image of the charge
distribution in the EL gap at a given time slice, with a similar relative
contribution as electron diffusion.

Figure 3: Simulated 136Xe $0\nu\beta\beta$ event. Left: initial track; Right -
same track after drifting 40 cm considering the longitudinal and transverse
diffusion coefficients to be 0.27 mm$/\sqrt{\textrm{cm}}$ and 1.07
mm$/\sqrt{\textrm{cm}}$ respectively. Color scale in the right plot represents
the number of ionization electrons at each $xy$ bin.

Both blurring effects can be characterized using 83mKr data. In fact,
transverse and longitudinal diffusion in NEXT-White were already studied using
83mKr events in Simon:2018vep . The r.m.s. transverse diffusion spread at 10
bar was found to be 1.07 mm${\times\sqrt{z(\textrm{cm})}}$ and the
longitudinal one – 0.27 mm${\times\sqrt{z(\textrm{cm})}}$. For the full drift
distance in NEXT-White, $z=53$ cm, the transverse and longitudinal FWHM spread
of an electron cloud starting from a point-like charge distribution are, in
this case, 18.3 mm and 4.6 mm, respectively. An example for the effect of
diffusion on a simulated $0\nu\beta\beta$ event for a drift distance of 40 cm
is shown in figure 3.

Both electron diffusion and the EL light spread can be quantified in terms of
point spread functions. The full diffusion PSF is three dimensional: a point-
like initial electron cloud transforms after diffusion to an oblate 3D
Gaussian (wider in the transverse plane than along the drift direction), where
both the transverse and longitudinal widths are proportional to $\sqrt{z}$.
This 3D PSF can be projected on the $xy$ plane to yield an effective 2D
transverse diffusion PSF, $F_{dif}^{2D}(x^{\prime},y^{\prime};z)$; (here
$x^{\prime}$ and $y^{\prime}$ are the $xy$ coordinates in a frame of reference
centered on the PSF axis). Similarly, integrating the total light hitting the
tracking plane for a point-like charge crossing the EL gap produces a 2D EL
PSF, $F_{EL}(x^{\prime},y^{\prime})$. Unlike the diffusion PSF, the EL PSF
does not depend on the drift distance $z$. Detailed analysis of 83mKr events
show that except for the TPC edges, both the diffusion and EL PSFs do not
depend, to leading order, on the absolute $xy$ position with respect to the
TPC axis, and both are axisymmetric.

Experimentally, the EL PSF can be determined from 83mKr events occurring
immediately in front of the gate, such that they do not suffer a diffusive
spread and can be considered point-like. The procedure, similar to that
described in Simon:2018vep , involves recording a large number of 83mKr events
over a small drift region (drift time $<25\;\mu$s). For each event, the SiPM
response is integrated over several $\mu$s, to include the full S2 signal. The
$xy$ location of the event is determined by calculating the center of gravity
(COG) of the SiPM hit map. The coordinates of all SiPMs participating in the
event are shifted to a reference frame whose origin coincides with the COG of
the event, and their charge is binned in $1\times 1$ mm2 pixels. The process,
repeated over a large number of events, such that in each step the SiPM charge
is added to the corresponding pixels, converges to the PSF. Figure 4 (top
left) shows the profile of the EL PSF constructed from 83mKr data. For
comparison, we also show the PSF extracted from a Monte Carlo simulation. The
MC PSF is narrower and the wings of the experimental PSF are somewhat larger.
This indicates that the present simulation does not provide a complete
description of the optical processes occurring in the EL gap and multi-layered
anode plate (in particular, it does not include, at this stage, a contribution
of secondary photoelectron emission from the gate, which may lead to single-
electron EL signals at some distance from the main event).

The effective 2D transverse diffusion PSF is given approximately by a Gaussian
whose standard deviation is: $\sigma_{t}=1.07\sqrt{z}$, where $z$ is in cm and
$\sigma_{t}$ in mm. This function is shown in the top right panel of figure 4,
for several values of $z$. The two blurring effects of diffusion and EL light
production can be combined in a single z-dependent PSF, which, to a good
approximation, is given by their convolution:
$F_{dif+EL}(z)=F_{dif}^{2D}(z)*F_{EL}$. This expression is only approximate,
because in reality longitudinal diffusion introduces “cross talk” of charge
between adjacent slices; however, since longitudinal diffusion is $\sim 4$
times smaller than transverse diffusion, and is, in all cases, smaller than
the width of the EL gap, we do not take into account this small effect. The
combined PSF can be determined experimentally similarly to the EL PSF, by
selecting 83mKr events over a range $[z,z+\Delta z]$, as was done in
Simon:2018vep  to determine the transverse and longitudinal diffusion
coefficients. In the present study, events were selected in 25 $\mu$s drift
time intervals. This step size was chosen to keep variations in transverse
diffusion within the interval sufficiently small. Thus, for the first 100
$\mu$s drift, the relative change within a 25 $\mu$s interval is $\sim$10%,
while at half the chamber, where the drift time is $\sim$225 $\mu$s, the
change falls below 5%. The bottom panel of figure 4 shows the combined
z-dependent PSF for the set of drift distances shown in the top right panel.
We show both the experimental PSF and the one obtained by convolving the Monte
Carlo EL PSF and the transverse diffusion Gaussian. As indicated by the
figure, convolution of the EL PSF with the diffusion PSF washes out most of
the differences between the two. Note that for the two different data sources,
MC and detector acquired data, double escape peak events are analyzed below
with their respective z-dependent PSFs.

Figure 4: Top left: EL PSF (from 83mKr data and Monte Carlo). Top right: The
2D transverse diffusion PSF (calculation based on the diffusion measurements
in Simon:2018vep ). Bottom left: The combined EL+diffusion PSF. Experimental
data (solid line) are from 83mKr events selected by their drift distance,
while MC (dashed lines) is a convolution of the EL PSF and the 2D transverse
diffusion PSF. Bottom right: The 2D PSF for experimental data for the
different drifts shown in the accompanying plots, along with the 2D EL PSF in
the first quadrant. The intensity value refers to the area-normalized value of
the distributions.

##  4 Track reconstruction in NEXT-White using Richardson-Lucy deconvolution

The Richardson-Lucy (RL) algorithm (also known as the Lucy-Richardson
algorithm) was developed independently by W. H. Richardson Richardson:72  and
L. B. Lucy Lucy:1974yx  in the early 1970’s in the context of observational
astronomy. Richardson focused his discussion solely on the recovery, by
deconvolution, of an underlying sharp image from an observed blurred one.
Lucy’s work, in contrast, considered a more general case, where one seeks to
recover an underlying frequency distribution from an observed one, with image
restoration as a particular application. The method has subsequently become a
main tool for image restoration in many scientific and engineering fields. The
algorithm is iterative, generating a sequence of improved approximations for
the underlying sharp image based on the observed blurred and noisy one, and
the (presumably known) point spread function. It can be applied on an
arbitrary number of spatial dimensions. Appendix A outlines the procedure for
two-dimensional images. Here we describe its implementation within the NEXT
event reconstruction scheme and a series of evaluations performed to cross-
check and validate the methodology.

###  4.1 RL implementation in NEXT-White

In NEXT-White, we apply RL deconvolution on individual SiPM time-sliced hit
maps (each integrated over a time interval $\delta t$, where, in the present
analysis $\delta t=2$ $\mu$s). Each slice is considered to be fully
independent from the others and longitudinal spread is not taken into account.
For a slice recorded at time $t$, we associate a physical slice of width
$\delta z=v_{d}\delta t$ of the original 3D track at the corresponding drift
distance $z=v_{d}\cdot(t-t_{0})$ (for $v_{d}=0.9$ mm/s, $\delta z=1.8$ mm). We
identify the physical slice, using the terminology of appendix A, with the
underlying sharp 2D image $W(x,y)$, and the corresponding SiPM hit map as a
sampled representation of the blurred image $\tilde{H}(x,y)$. The two images
are assumed to be related through the combined diffusion+EL PSF $F(x,y;z)$
corresponding to a drift distance $z$.

The implementation of RL deconvolution on each slice is done, for both
experimental data and MC events, in the following steps:

  1. 1.

Charge cut: SiPMs containing less charge than a predefined threshold $q_{cut}$
are removed from the slice. As discussed in appendix B, after exploring
several different choices for the charge threshold, we adopted a value
$q_{cut}=10$ PE (compared to 30 PE used in the classical analysis).

  2. 2.

Removal of isolated SiPMs: Single SiPM hits which have no adjacent non-zero
neighbors in the same slice are removed. This is done to avoid filling the
region between the main track and isolated SiPMs which fluctuate above the
charge threshold $q_{cut}$ by non-physical data in the subsequent
interpolation step.

  3. 3.

2D interpolation: We define a rectangular region surrounding the SiPMs which
have survived steps (1) and (2) with 10 mm margins. To estimate the full
pattern of photon hit points in this region, we apply bicubic 2D interpolation
on the “cleaned up” SiPM hit map over a $1\times 1$ mm2 grid (the SiPMs cover
only 1$\%$ of the plane). Note that no significant differences were observed
in the final outcomes of the analysis (signal efficiency and background
acceptance for double escape peak events) when replacing bicubic by linear
interpolation (see a brief discussion in appendix B).

  4. 4.

RL deconvolution: For each slice, we use the corresponding z-dependent
combined EL+diffusion PSF for the deconvolution process, following equations
(16)-(18) in appendix A, to find successive estimations $W^{(r)}(x,y)$ for
$W(x,y)$ in $N_{iter}$ iterations. Data and MC PSFs are used for data and MC
events, respectively. The process maintains the overall charge of each slice
constant in all iterations. It was implemented using the Richardson-Lucy
function from Python’s scikit-image library scikit-image .

  5. 5.

Cleaning cut: Once the iterative process is completed, a cleaning cut with an
adjustable threshold $\epsilon_{cut}$ is applied to the image intensity given
by the iterative process. This is done to remove non-physical backgrounds and
reconstruction leftovers, and sharpen the track edges for the topological
analysis. For the double escape peak analysis the cut was set at 0.008 a.u.
Details on the optimization can be found in appendix B. No cut was performed
when applying the method to Kr events as the reconstruction leftovers where
not found to have an impact on the performance.

  6. 6.

Energy allocation: Finally, based on the integrated S2 signal recorded by the
PMTs over the entire event duration, and using the 83mKr-based lifetime and S2
correction maps Renner:2019pfe , we find the total energy of each recorded
slice and divide it among all of the $1\times 1$ mm2 pixels of the deconvolved
image, proportionally to their interpolated charge.

The first three steps aim to generate a reasonable estimate for the actual
photon hit pattern on the tracking plane. They reflect a pragmatic approach to
bridge the empty spaces (and hence lack of information) between the SiPMs, and
to avoid distorting the image by distant effects, such as reflections from the
various TPC surfaces, or distant EL light emission by photoelectrons ejected
from the gate mesh. The interpolation step is justified as the smearing
effects of both electron diffusion and EL light spread produce gradual changes
in light intensity on the tracking plane. Rather than claim for absolute
mathematical rigor in this approximation, we provide a series of simple
demonstrations to support its practical value. These include the
reconstruction of individual Kr events, adjacent pairs of Kr decays, and
straight muon tracks.

###  4.2 Validation tests

As a first test, RL deconvolution was applied to individual Kr events from
both MC and detector data. A typical example, from data, is shown in figure 5,
with the SiPM sensor response for $q_{cut}=10$ PE on the first column, bicubic
interpolation on the second, and deconvolved images after 75 RL iterations on
the third. When comparing the deconvolved images to the MC true information,
the r.m.s. error in the reconstructed COG of all Kr-events was $\sim 0.5-1.0$
mm in both $x$ and $y$. For 75 RL iterations, the FWHM of the reconstructed
Gaussian-like charge distribution was $\sim 5$ mm for both MC and detector
data, independent of the drift distance, for $z>100$ mm. Although adding
iterations was found to reduce the FWHM further, the effect was quite marginal
(e.g., applying 150 iterations reduced the FWHM to $\sim 4$ mm).

Figure 5: Example of a reconstructed 83mKr event from NEXT-White data. The
event is centered at (0,0) for convenience. Left: raw sensor response, with a
charge cut of 10 PE. Center: bicubic interpolation. Right: deconvolved image
after 75 RL iterations.

Figure 6: Reconstructed pair of nearby 83mKr events from detector data,
synthesized by overlaying SiPM response maps for two nearby events. Left
column - raw (synthesized) SiPM data; center - bicubic interpolation; right -
deconvolved image after 75 iterations. Red dots represent the COGs of
individual events.

.

While the accuracy in COG reconstruction using RL deconvolution for point-like
Kr events was practically the same as obtained from the raw SiPM response, the
strength of the method lies in the ability to resolve nearby structures. To
demonstrate this, we artificially synthesized pairs of Kr events from detector
data, by overlaying SiPM sensor response maps for two individual events with
the same drift distance. Image synthesis was done by shifting the SiPM
response map of one event by an integer number of 10 mm steps in $x$ and $y$,
to bring its COG close to that of the other one. An example is shown in figure
6, where the COGs are 19.5 mm apart (for a drift distance of $\sim 360$ mm).
The dots represent the individual Kr COGs. As in figure 5, the left, center
and right columns show the raw (overlaid) sensor response maps, interpolated
images and deconvolved ones (with 75 iterations).

Figure 7: Two muon events in NEXT-White data after RL deconvolution, shown in
three Cartesian projections. The top event is of a “clean” muon, while the
bottom one also contains a delta electron.

To demonstrate the method over long tracks, we selected, from detector data, a
set of crossing muon events. A visual inspection was carried over an extensive
dataset with no observable deviations from the expected straight line tracks,
with occasional delta electrons branching out from the main track. Two
examples of muon tracks of NEXT-White data (Run V) are shown in figure 7.

Sensor response

2 iterations

10 iterations

75 iterations

Figure 8: Effect of successive RL iterations (2, 10 and 75) in three selected
slices. The event is an electron-positron pair (from data).

Figure 8 demonstrates the iterative refinement in the track sharpness obtained
by the RL process. The event consists of a 1.6 MeV electron-positron pair,
acquired in NEXT-White, which forms a U-shaped track, where the two ends are
located roughly in the same $xy$ plane. The first row shows the raw data
(sensor response) before deconvolution, binned in $10\times 10\times 1.8$ mm
voxels. The full track, projected on the $yz$ plane, is shown on the top left.
The dashed lines mark three particular slices where the effect of applying
successively more iterations is shown for each one individually. The
corresponding slice images are shown in columns 2-4. Rows 2-4 show the same
event after 2, 10 and 75 iterations, with the latter after application of the
cleaning cut. The disconnected artifact appearing in slice 1 after 75
iterations (around $x=70$ mm, $y=50$ mm) does not affect the analysis of the
main track, as explained in the following section.

##  5 Analysis of 1.6 MeV e-e+ events

###  5.1 Methodology

The detector data used for this analysis were taken during August 2019 under
the same conditions as in Ferrario:2019kwg , but with a much longer electron
lifetime of $\sim 9$ ms. The MC simulation dataset was exactly the same.
Detector data consisted of events generated by a 228Th source placed in an
external calibration port above the center of the drift region, and of
internal 83mKr events.

Similarly to Ferrario:2019kwg , event selection for the topological analysis
around the 208Tl double escape peak was done by applying the following
filters: (1) event energy in the range $1.4-1.8$ MeV; (2) full containment of
the event in the fiducial volume, with all hits at least 2 cm away from all
borders of the drift region; (3) the event comprises only a single track, when
binned into 15 mm voxels. The latter filter, to which we refer below as the
“gross single-track cut” (given the large voxels size), was employed on the
SiPM data before deconvolution rather than at the end of the process, to avoid
artificial track multiplicity which may be created during the RL iterations
and as a result of the final cleaning cut.

Events passing the above filters were processed according to the RL steps
listed in section 4. The set of deconvolved 2D images (which, in fact,
consisted of $1\times 1\times 1.8$ mm3 voxels with 1.8 mm along $z$) were re-
binned into larger voxels of adjustable size. As discussed in appendix B, 5-mm
voxels provided the optimal results. The BFS algorithm was then used – as in
the classical analysis – to combine the event voxels into one or more tracks.
At this stage, events may have contained more than one track either because
they consisted of multiple nearby physical tracks which were unresolved by the
gross single-track cut, or due to artificial breaking of a single track into
smaller segments during the RL process. Concretely, it was estimated, based on
the true information of simulated data, that 29% of the events passing the
gross single-track cut were actually multi-track events. At the same time, for
the optimal parameter configuration described below, 27% of the deconvolved
simulated events exhibited more reconstructed tracks than the number of
segments in the true information. These additional “satellites” were of low
intensity, with essentially no effect on the blob-based analysis.

Once the voxels were grouped into tracks, the BFS algorithm was used to find
the track ends. Those of the longest track of the event served as the centers
for two spherical blobs of a variable radius. As before, we defined blob1 as
the more energetic blob and blob2 as the less energetic one. We set a
threshold $E_{0,blob2}$ on blob2 energy such that events above it were
considered signal and below it – background. For the $i$-th value of blob2
energy threshold $E_{0,blob2}^{(i)}$ we defined the signal efficiency
$\epsilon_{i}$, background acceptance $b_{i}$ and figure of merit $f.o.m_{i}$
as:

$\epsilon_{i}=\frac{\mathrm{number\;of\;signal\;events\;with\;E_{blob2}>E_{0,blob2}^{(i)}}}{\mathrm{total\;number\;of\;signal\;events\;with\;no\;cut}}$
(1)

$b_{i}=\frac{\mathrm{number\;of\;background\;events\;with\;E_{blob2}>E_{0,blob2}^{(i)}}}{\mathrm{total\;number\;of\;background\;events\;with\;no\;cut}}$
(2)

$f.o.m_{i}=\frac{\epsilon_{i}}{\sqrt{b_{i}}}$ (3)

For MC, the true nature of the event is known: any event containing a positron
is signal, and any event without one is background. For each threshold
$E_{0,blob2}^{(i)}$, the signal efficiency, background acceptance and figure
of merit are found directly from equations (1)-(3) using the true information.
For experimental detector data, however, the nature of the event is unknown a
priori and one must resort to a different approach, which involves fitting the
ROI data around the double escape peak with an expression that describes both
signal and background, before and after the application of the cut
$E_{blob2}>E_{0,blob2}^{(i)}$ Ferrario:2019kwg .

Figure 9: Energy spectrum and event population in the double escape peak
region for MC events. In blue, events with a positron, in red – the rest of
the interactions (predominantly Compton electrons). Left: the spectrum before
applying any topological cut. Right: after a blob2 energy cut of 340 keV.

Figure 9 shows the MC energy spectrum around the double escape peak,
highlighting separately the contribution of $e^{-}e^{+}$ pairs (in blue) and
of events with no positron (in red). The figure shows the spectrum before and
after the application of the topological cut. The background spectrum in this
region can be well described by a decreasing exponential, as was done in
Ferrario:2019kwg :

$f_{bkg}(E)=A_{1}\exp(-A_{2}E)$ (4)

The signal spectrum consists of a Gaussian centered at 1593 keV, with flat
shallow wings for energies below and above the peak, as shown in figure 10.
Below the peak, signal events comprise $e^{-}e^{+}$ pairs created by photons
of energy below 2615 keV (which predominantly result from Compton scatters of
2615 keV gammas prior to pair production), as well as $e^{-}e^{+}$ pairs which
lose some energy by bremsstrahlung, where the emitted photon does not interact
in the sensitive volume. Above the peak there is a smaller population of
$e^{-}e^{+}$ events where one of the 511 keV gammas created when the positron
annihilates interacts close to the main track and is unresolved by the gross
single-track cut. The signal spectrum can therefore be approximated as:

$f_{sig}(E)=B_{1}\left(\frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{\left(E-\mu\right)^{2}}{2\sigma^{2}}\right)+C_{1}\mathrm{erfc}\left(\frac{E-\mu}{\sqrt{2}\sigma}\right)+C_{2}\right)$
(5)

where a complementary error function (erfc) with the same standard deviation
as the Gaussian is used to describe events below the peak. This expression was
chosen empirically based on the MC distribution shown in figure 10 and should
be considered only as a proxy to obtain an estimate of the signal population
outside the peak with no other physical meaning. The parameters $C_{1}$ and
$C_{2}$ were extracted from fitting, through an unbinned extended maximum
likelihood fit, MC data consisting only of $e^{-}e^{+}$ pairs in the region
$1.4-1.8$ MeV, giving: $C_{1}=(4.89\pm 0.46)\cdot 10^{-4}$ and $C_{2}=(1.99\pm
0.28)\cdot 10^{-4}$. Note that the analysis in Ferrario:2019kwg  disregarded
the shallow wings in the $e^{-}e^{+}$ spectrum, and instead assumed that the
signal is completely described by a Gaussian. This assumption is valid when
the signal population outside the peak is negligible compared to the
background, which is the case before application of the topological cut
(figure 9, left). For a modestly effective topological analysis, this holds
also after applying the cut. However, this consideration loses validity as the
topological cut becomes more effective, and a realistic fit should consist of
the sum of equations (4) and (5):

$f(E)=f_{bkg}(E)+f_{sig}(E)$ (6)

Figure 10: Energy spectrum of $e^{-}e^{+}$ pairs around the 1.6 MeV double
escape peak, fitted according to equation (5).

The procedure for estimating the signal efficiency and background acceptance
using the fit was as follows. The parameters $C_{1}$ and $C_{2}$ were
extracted from the MC fit to the true signal spectrum and assumed to hold also
for detector data, without being affected by the application of the
topological cut. This approximation, which reflects an assumption that the
signal efficiency of the the blob cut is energy-independent over the range
$1.4-1.8$ MeV, was confirmed to hold on MC data. For each value of
$E_{0,blob2}^{(i)}$ we then fitted the data with $f(E)$ (equation (6)), with
$A_{1}$, $A_{2}$ and $B_{1}$ as free parameters, before and after applying the
cut $E_{blob2}>E_{0,blob2}^{(i)}$ using an unbinned maximum likelihood fit
(the parameters $\mu$ and $\sigma$ were extracted once from a fit to the
Gaussian alone and kept constant in subsequent fits). The number of background
and signal events passing a given value of $E_{0,blob2}^{(i)}$ were calculated
by integrals over $\pm 3\sigma$ around the peak centroid:

$N_{bkg,i}=\int_{\mu-3\sigma}^{\mu+3\sigma}A_{1,i}\exp(-A_{2,i}E)dE$ (7)

$N_{sig,i}=\int_{\mu-3\sigma}^{\mu+3\sigma}B_{1,i}\left(\frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{\left(E-\mu\right)^{2}}{2\sigma^{2}}\right)+C_{1}\mathrm{erfc}\left(\frac{E-\mu}{\sqrt{2}\sigma}\right)+C_{2}\right)dE$
(8)

The estimate was found to be within 1% of the real number of events for both
signal and background populations in MC. The signal efficiency and background
acceptance for the cut $E_{blob2}>E_{0,blob2}^{(i)}$ were calculated as:

$\epsilon_{i}=N_{sig,i}/N_{sig,0}$ (9)

$b_{i}=N_{bkg,i}/N_{bkg,0}$ (10)

where the subscript “0” refers to no cut (i.e., $E_{0,blob2}=0$). The figure
of merit for the $i$-th threshold was calculated, as before, by
$f.o.m_{i}=\epsilon_{i}/\sqrt{b_{i}}$.

Figure 11: Experimental detector data events in the 1.6 MeV double escape peak
region after topological classification. Top two rows: $e^{-}e^{+}$
candidates; bottom two rows: background candidates.

Figure 11 shows four examples of detector data events reconstructed using RL
deconvolution. The first two rows show events classified as signal, while the
third and fourth rows display events classified as background (deconvolution
and classification were done using the choice of parameters described below).
Note that the second background event contains a blob-like satellite close to
its start point. This event passed the gross single-track cut, and, in the
classical analysis, would have likely been classified as signal. The refined
reconstruction offered by RL deconvolution allows separating this satellite
from the main track and correctly classify this event as background.

###  5.2 Results

The RL process and subsequent analysis involve several parameters, whose
values determine the blob identification quality which ultimately reflects on
the signal efficiency and background acceptance. These parameters can be
divided in two groups: those affecting the quality of the reconstructed image
and the end-point identification accuracy, and those which affect the energy
calculation inside the blobs. The first group is related to the RL process and
includes the SiPM charge threshold $q_{cut}$, number of RL iterations
$N_{iter}$ and final cleaning threshold $\epsilon_{cut}$. The second group is
comprised of the BFS characteristic parameters, namely the voxel size
$l_{voxel}$ and blob radius $R_{blob}$. For a given choice of parameters the
$f.o.m.$ attains a maximum for a particular value of $E_{0,blob2}$. By
carefully adjusting the values of these parameters one can try to maximize the
$f.o.m.$, and therefore the experimental sensitivity in the energy region
under investigation.

The optimization process is highly demanding in computing resources: beyond
being a multivariate problem, the determination of the $f.o.m.$ for every
parameter configuration requires analyzing tens of thousands of events. We
therefore adopted here a pragmatic approach, which allowed attaining a high
value for the $f.o.m.$, but without proving that this is a global maximum. The
analyzed dataset comprised $3.8\cdot 10^{4}$ MC and $1.9\cdot 10^{4}$ detector
data events passing the energy, fiducialization and gross single-track cuts.
We performed separate analyses to MC events using the true information
(equations (1)-(3)), MC events where signal and background estimates were made
using the fitting procedure (equations (7)-(10)), and experimental detector
data using the same fitting procedure. The optimal choice of parameters and
resultant $f.o.m.$ vary, to some extent, between MC and detector data.
Appendix B provides a detailed description of the optimization process. Here,
we present the final results of the analysis.

Our choice of parameters aimed to simultaneously achieve a high $f.o.m.$ for
both detector data and MC (using the fitting procedure) while looking for
optimal MC (true information) performance. We refer to the chosen set of
parameters as the “optimal configuration,” where: $q_{cut}=10$ PE,
$N_{iter}=75$, $\epsilon_{cut}=0.008$, $l_{voxel}=5$ mm and $R_{blob}=18$ mm.

Figure 12 shows 2D histograms of signal and background events from detector
data, binned according to their blob1 and blob2 energies, comparing the
classical analysis (A) to RL deconvolution with the optimal choice of
parameters (B). Each histogram contains both signal (electron-positron) and
background (Compton scatter) events, in the range $1.4-1.8$ MeV, where the
signal population occupies a higher region in the blob1-blob2 plane. The
horizontal white lines mark the choice of blob2 energy cut which maximizes the
$f.o.m$. Application of RL deconvolution clearly leads to an improved
separation between the two groups.

Figure 12: 2D histogram of blob1 and blob2 energies of 1.6 MeV double escape
peak events from detector data: (A) classical analysis; (B) RL deconvolution
with the optimal choice of parameters. Each histogram contains events passing
the fiducial cut and gross single-track cut with energies in the range
$1.4-1.8$ MeV, with both signal and background populations (with higher and
lower blob2 energies, respectively). The horizontal white lines mark the
choice of blob2 energy cut which maximizes the $f.o.m$.

Figure 13A shows the signal efficiency as a function of background rejection
($1-b$) for 1.6 MeV double escape peak events, for the optimal choice of
parameters. The figure includes the curves for detector data, MC (using the
fit and true information) and – for comparison – the classical analysis of
detector data events Ferrario:2019kwg . Increasing background rejection is
equivalent to moving the horizontal white lines in figure 12 upward. Figure
13B shows the $f.o.m.$ as a function of the threshold on blob2 energy for the
same datasets and parameters. For RL deconvolution the maximal (optimal)
$f.o.m.$ is obtained for $E_{0,blob2}=340$ keV.

Figure 13: (A) Signal efficiency vs. background rejection for 1.6 MeV double
escape peak events, for detector data, MC (using the fit) and MC (using the
true information), for the optimal choice of parameter configuration (marked
rectangles in figure 18 in appendix B). The curve resulting from the classical
analysis of data (from Ferrario:2019kwg ) is shown for comparison. (B) The
figure of merit for the optimal parameter configuration for detector data, MC
(fit and true) and classical analysis vs. the threshold on blob2 energy. The
maximal (optimal) $f.o.m.$ is for a blob2 energy cut at 340 keV.

Table 1 shows the signal efficiency, background acceptance and $f.o.m.$ for
detector data and MC (with the classification based on both the fit and true
information), for the optimal choice of parameters and optimal value of blob2
energy cut. For comparison the table also includes the results of the
classical analysis. The errors represent statistical uncertainties (standard
deviation); systematic effects were found to be dominated by the choice of
reconstruction parameters and are not further considered here.

Table 1: Signal efficiency and background acceptance for the optimal figure of merit. Dataset/analysis | Signal efficiency | Background acceptance | Figure of Merit
---|---|---|---
Classical | $71.6\pm 1.5\%$ | $20.6\pm 0.4\%$ | $1.58\pm 0.04$
Data | $56.6\pm 2.2\%$ | $3.7\pm 0.7\%$ | $2.94\pm 0.28$
MC (fit) | $59.4\pm 1.6\%$ | $4.0\pm 0.5\%$ | $2.98\pm 0.20$
MC (true) | $59.4\pm 1.0\%$ | $3.3\pm 0.4\%$ | $3.26\pm 0.20$

For detector data, the RL-based analysis using the optimal choice of
parameters provides a 5.6-fold reduction of background acceptance compared to
the classical analysis (overall topological background rejection factor of
$\sim 27$), accompanied by a relative reduction of signal efficiency by 21%.
According to the simulated data (true MC information), a 6.2-fold reduction of
background and a 17% relative reduction of signal is achieved with the chosen
configuration.

Figure 14: (A) Reconstructed energy for blob2 for deconvolution and classical
analysis. The “true” blob energy is found by integrating over a sphere of 18
mm radius centered on the true track end-point. (B) Distributions divided by
population, either signal or background.

The enhanced background rejection power is a consequence of the improved blob
reconstruction, both in positioning and radial extent, enabled by RL
deconvolution. The RL-reconstructed blob energy distribution is much closer to
the expected one than that derived from the classical analysis, as illustrated
in figure 14. Only blob2 is shown as the difference in performance is more
pronounced than for blob1 (which also exhibits an improvement). Figure 15,
left, compares the overall performance of RL deconvolution and the classical
analysis in end-finding. The histograms show the distributions of distances
between the reconstructed end positions and the true ones (extracted from MC)
for all events in the 1.6 MeV ROI, combining signal and background with no
distinction between blob1 and blob2. Employing RL deconvolution clearly
reduces the errors, leading to improved blob-based classification. In appendix
B we show that, as expected, high f.o.m. values are correlated with reduced
errors in end-finding. The right panel of figure 15 focuses on the error
distributions for RL-reconstructed events, separating them by event and blob
type. It shows that the best results are obtained for blob2 in background
events, where there is no Bragg peak. While there is considerable reduction in
error for the other cases compared to the classical analysis, in about 30% of
the events the error is larger than the blob radius (18 mm), indicating that
there is still ample room for improvement.

Figure 15: End-finding error distributions. Left: overall error distributions
obtained by RL-deconvolution and the classical analysis for 1.6 MeV double
escape peak events; right: distributions for RL-deconvolution divided into
signal and background events, with subdivision into blob1 and blob2.

##  6 Discussion

Richardson-Lucy deconvolution was shown in this work to be a highly effective
tool for enhancing image reconstruction in NEXT, leading to a major
improvement in topological background rejection. The application of the method
to detector data in the 1.6 MeV double escape peak of 208Tl, using a cut on
blob2 energy, yielded a background acceptance level $b=(3.7\pm 0.7)\%$ with a
signal efficiency $\epsilon=(56.6\pm 2.2)\%$ (figure of merit
$f.o.m.=\epsilon/\sqrt{b}=2.94\pm 0.28$), with similar results for detector
data and simulated events. This represents a drastic improvement relative to
previous results by the collaboration (20.6% background acceptance with 71.6%
signal efficiency, $f.o.m.=1.58$ Ferrario:2019kwg ) and greatly boosts the
background rejection realizable by the experiment. The obtained level of
background acceptance is similar to the value previously reported by the
Gotthard experiment using visual inspection of double escape peak events at
1.6 MeV Gotthard1998 . The new results are also considerably better than those
obtained using a deep convoluted neural network on non-deconvolved tracks to
classify double escape peak events, where the background acceptance was 10%
and signal efficiency was 65% ($f.o.m.=2.06$) Kekic:2020cne .

The primary effect of employing RL deconvolution is the attainment of refined
3D track images, which allows better identification of the track ends, and
therefore improved positioning of the blob centers and better estimates of the
energy they contain. Improved end finding generally results from the enhanced
resolving power offered by RL-deconvolution. Several illustrative examples for
this are given in appendix C.

The focus of this work was on the analysis of experimental detector data
recorded at the 1.6 MeV 208Tl double escape peak. However, we also probed the
possibility of implementing the RL-based method to MC signal and background
events in the $Q_{\beta\beta}$ ROI, using the NEXT-White detector MC model.
Preliminary analysis indicated that similar levels of background acceptance
and signal efficiency are expected in this energy range. However, since
NEXT-100 will operate under different conditions (15 bar, 112 cm maximal
drift, 15.6 mm SiPM spacing and a different EL and tracking plane geometry),
such results are only of indicative nature, and a full simulation is deferred
to a separate publication.

In spite of the improved event classification offered by the new method, the
full potential Richardson-Lucy deconvolution is yet to be exploited by the
collaboration, as the results presented in this work remain a first-approach
evaluation using existing tools (e.g., the BFS algorithm), which may not be
optimal for the fine-grained output of the RL procedure. In particular,
although figure 15 shows a clear improvement in track-end finding, the
distributions of deconvolved events have significant tails extending to large
errors, requiring the use of large blob radii for the analysis. Presently,
several ideas for potential improvement in event classification are under
study. These include improved algorithms for end finding, as well as the
potential use of Machine Learning approaches for the classification of high-
definition reconstructed events. In addition, RL deconvolution can be further
developed to a full 3D method instead of following the slice-by-slice approach
employed in this work, which could lead to enhanced image quality along the
drift direction. Lastly, image quality may greatly benefit from the use of
low-diffusion gas mixtures, such as Xe-He Felkai2018 ; McDonald2019 ;
Fernandes2020 , or Xe “spiked” with low concentrations of molecular additives
Azevedo2016 ; Henriques2017 ; Henriques2019 .

Prior to further improvement of the method, the collaboration is presently
evaluating the benefits offered by RL deconvolution in its current form, in
particular for the analysis of $2\nu\beta\beta$ events in NEXT-White. As an
example, figure 16 shows a 2.0 MeV double beta candidate (from NEXT-White
detector data) reconstructed by the RL procedure described above.

Figure 16: Deconvolved 2.0 MeV $2\nu\beta\beta$ candidate obtained during the
current data taking of NEXT-White.

###### Acknowledgements.

The NEXT Collaboration acknowledges support from the following agencies and
institutions: the European Research Council (ERC) under the Advanced Grant
339787-NEXT; the European Union’s Framework Programme for Research and
Innovation Horizon 2020 (2014–2020) under the Grant Agreements No. 674896,
690575 and 740055; the Ministerio de Economía y Competitividad and the
Ministerio de Ciencia, Innovación y Universidades of Spain under grants
FIS2014-53371-C04, RTI2018-095979, the Severo Ochoa Program grants
SEV-2014-0398 and CEX2018-000867-S, and the María de Maeztu Program
MDM-2016-0692; the Generalitat Valenciana under grants PROMETEO/2016/120 and
SEJI/2017/011; the Portuguese FCT under project PTDC/FIS-NUC/2525/2014 and
under projects UID/04559/2020 to fund the activities of LIBPhys-UC; the U.S.
Department of Energy under contracts No. DE-AC02-06CH11357 (Argonne National
Laboratory), DE-AC02-07CH11359 (Fermi National Accelerator Laboratory), DE-
FG02-13ER42020 (Texas A&M) and DE-SC0019223 / DE-SC0019054 (University of
Texas at Arlington); the University of Texas at Arlington (USA); and the Pazy
Foundation (Israel) under grants 877040 and 877041. DGD acknowledges Ramon y
Cajal program (Spain) under contract number RYC-2015-18820. JM-A acknowledges
support from Fundación Bancaria “la Caixa” (ID 100010434), grant code
LCF/BQ/PI19/11690012. AS acknowledges support from the Kreitman School of
Advanced Graduate Studies at Ben-Gurion University.

##  Appendix A Richardson-Lucy deconvolution

The Richardson-Lucy (RL) algorithm aims to recover, by deconvolution, an
underlying sharp image from an observed blurred and noisy one. The algorithm
is iterative, generating a sequence of improved approximations for the
underlying image using the (presumably known) point spread function (PSF) of
the imaging process. In this appendix we outline the mathematical procedure,
as employed on 2D images.

We denote by $W(x,y)$ the underlying sharp image and by $F(x,y)$ the PSF. In
the absence of noise, the ideal blurred image $H(x,y)$ is obtained as a
convolution of $W$ and $F$:

$H(x,y)=\iint
W(x^{\prime},y^{\prime})F(x-x^{\prime},y-y^{\prime})dx^{\prime}dy^{\prime}$
(11)

In principle, $W(x,y)$ could be recovered from $H(x,y)$ by solving this
integral equation. This could be done by discretization, converting $H$, $W$
and $F$ into matrices and the integral to a double summation. This generates a
system of linear equations with the elements of $W$ as the unknowns:

$\sum_{j}F_{i,j}W_{j}=H_{i}$ (12)

where the indices $i$ and $j$ refer to individual elements in $H$ and $W$,
combining the enumeration of both the row and column, and $F_{i,j}$ describes
the influence of the $j$-th element of $W$ on the $i$-th element of $H$.

In reality, because of the presence of noise, the actual observed image
$\tilde{H}(x,y)$ is different from the ideal one $H(x,y)$. The system of
linear equations then becomes:

$\sum_{j}F_{i,j}W_{j}=\tilde{H}_{i}$ (13)

Attempting to solve equations (13) generally yields poor results, with large
discontinuities in $\\{W_{j}\\}$, as well as non-physical negative values.
This occurs because the process tends to amplify short-wavelength errors in
$\tilde{H}$, which are characteristic of noisy images Lucy:1974yx .

The approach of the RL algorithm is different. It begins by noting that one
could formally write:

$W(x,y)=\iint
H(x^{\prime},y^{\prime})G(x-x^{\prime},y-y^{\prime})dx^{\prime}dy^{\prime}$
(14)

provided that we define the inverse kernel $G$ as:

$G(x-x^{\prime},y-y^{\prime})=\frac{W(x,y)F(x-x^{\prime},y-y^{\prime})}{H(x^{\prime},y^{\prime})}$
(15)

where $\iint F(x-x^{\prime},y-y^{\prime})dx^{\prime}dy^{\prime}=1$. Since $G$
depends on $W$, the direct calculation of $W$ from equation (14) is
impossible. However, the process may work iteratively, if one could provide
successively improved approximations for $G$, which, in turn, would rely on
successive estimates of $W$.

We begin with an initial estimate $W^{(0)}$ for $W$, where $W^{(0)}(x,y)$ is
generally taken to be a flat image. In the $r$-th iteration we calculate an
intermediate blurred image $H^{(r)}$ by:

$H^{(r)}(x,y)=\iint
W^{(r)}(x^{\prime},y^{\prime})F(x-x^{\prime},y-y^{\prime})dx^{\prime}dy^{\prime}$
(16)

This allows finding an estimate for $G$:

$G^{(r)}(x-x^{\prime},y-y^{\prime})=\frac{W^{(r)}(x,y)F(x-x^{\prime},y-y^{\prime})}{H^{(r)}(x^{\prime},y^{\prime})}$
(17)

The new estimate for $W$, $W^{(r+1)}$, is then calculated following equation
(14), with $\tilde{H}$ replacing $H$ and $G^{(r)}$ replacing $G$:

$\begin{split}W^{(r+1)}(x,y)&=\iint\tilde{H}(x^{\prime},y^{\prime})G^{(r)}(x-x^{\prime},y-y^{\prime})dx^{\prime}dy^{\prime}\\\
&=W^{(r)}(x,y)\iint\frac{\tilde{H}(x^{\prime},y^{\prime})}{H^{(r)}(x^{\prime},y^{\prime})}F(x-x^{\prime},y-y^{\prime})dx^{\prime}dy^{\prime}\end{split}$
(18)

The discussion in Lucy:1974yx  shows that if successive changes in $W^{(r)}$
are sufficiently small, in the limit $r\rightarrow\infty$ the scheme converges
to the solution of the corresponding maximum likelihood problem222The
possibility of track reconstruction using a maximum likelihood approach had
been previously explored by the NEXT Collaboration Simon:2017pck , but was
later disfavored for the method presented here.. It further shows that if the
number of elements in $H$ is equal or larger than those of $W$ (i.e., the
system of linear equations (13) is over-determined), this solution is unique.

##  Appendix B Parameter optimization

In this appendix we describe in detail the steps taken to optimize the choice
of parameters used for the RL process and subsequent analysis.

The first step in the analysis was to choose a value for the SiPM charge
threshold, $q_{cut}$. It determines how much of the signal is cleaned out
before starting the deconvolution process. Cutting too low may lead to the
inclusion of distant signals (reflected light or photons induced by
photoelectrons emitted from the gate), while cutting too high may bias and
distort the output. The impact of $q_{cut}$ on the $f.o.m.$ was studied for
several cut values, over the range $5-25$ PEs (in 2 $\mu$s). We performed this
scan keeping $N_{iter}=90$, $\epsilon_{cut}=0.008$ (in arbitrary units),
$l_{voxel}=5$ mm and $R_{blob}=21$ mm (these values were chosen as a
reasonable starting point after visual inspection of many events). For each
value of $q_{cut}$ we calculated the $f.o.m.$ as a function of blob2 energy
threshold and found its maximal value. Given the results, shown in figure 17,
we settled on a 10 PE cut due to a much better match between detector data and
both MC fitted data and true information.

Next, we considered the effect of the number of RL iterations and final
cleaning cut. Understanding the optimal point to stop applying RL iterations
is of prime importance. If not applied enough times, the reconstructed charge
distribution remains too blurry, which harms the blob energy estimation. On
the other hand, over-iterating can result in noisy artifacts and in breaking
up of the track to disconnected segments. These effects are strongly related
to the subsequent application of the cleaning cut. If applied correctly, it
can remove artifacts that appear in the iterating process. However, this is a
delicate parameter as a too high cut could lead again into track
fragmentation.

Figure 17: Variation of the maximal $f.o.m.$ achieved for both experimental
detector data and MC for different cuts on the input signal. The scan on
$q_{cut}$ was done with fixed values for the other parameters, as described in
the text.

Given the observed relation between the number of iterations and the cleaning
cut, we decided to scan and optimize both parameters simultaneously. The
number of iterations was varied in steps of 15, and the cleaning cut was
scanned over the range $0.005-0.017$. For this scan, we fixed $l_{voxel}=5$ mm
and $R_{blob}=21$ mm. For each configuration we found the maximal $f.o.m.$ as
a function of $E_{0,blob2}$. The results are shown in the top part of figure
18 for detector data, MC using the fit and MC using the true information. We
chose a configuration which displayed a high $f.o.m.$ for both data and MC,
including the true information, namely $N_{iter}=75$ and
$\epsilon_{cut}=0.008$, yielding $f.o.m.(data)=2.77$, $f.o.m.(MC\;fit)=2.89$
and $f.o.m.(MC\;true)=3.11$ (marked rectangles).

DataMCTrue

Figure 18: Parameter scan for optimizing the figure of merit, equation (3),
for events in the double escape peak ROI for data, MC using the fitting
procedure, and MC using the true information (i.e., whether the event contains
a positron or not). The result for the optimal configuration is marked in red.
Top row: maximal figure of merit for different combinations of the number of
RL deconvolution iterations and the final threshold for the cleaning cut,
keeping a voxel size of 5 mm and blob radius of 21 mm. Bottom row: maximal
figure of merit for different combinations of voxel size and blob radius, for
75 RL iterations and a cleaning cut of 0.008. The SiPM charge cut is 10 PE in
all cases.

With the above choice of $N_{iter}=75$ and $\epsilon_{cut}=0.008$ we moved to
testing the effect of voxel size and blob radius. If voxels are too small, the
track may be broken into disconnected segments, resulting in identifying
internal points as the track ends. On the other hand, if voxels are too large,
one loses the advantages of track refinement through RL deconvolution: the
track is “re-smeared” and its ends are shifted. The blob radius is intimately
related to the choice of voxel size. The blob must be large enough to contain
the full energy deposited by the electron as it approaches the Bragg peak
(which is shared between several voxels), but not too large, as this can lead
to the inclusion of energy outside of the peak. Keeping $q_{cut}=10$ PE,
$N_{iter}=75$ and $\epsilon_{cut}=0.008$, we scanned $l_{voxel}$ from 3 to 10
mm, and $R_{blob}$ from 12 to 27 mm. As before, for each parameter
configuration we scanned the threshold on blob2 energy and found the maximal
$f.o.m.$ The results are shown in the second row of figure 18, with rectangles
marking the final choice of parameters, namely $l_{voxel}=5$ mm and
$R_{blob}=18$ mm. We refer to this choice ($q_{cut}=10$ PE, $N_{iter}=75$,
$\epsilon_{cut}=0.008$, $l_{voxel}=5$ mm and $R_{blob}=18$ mm) as the
“optimal” choice of parameters.

Signal

AllBlob1Blob2

Background

Figure 19: Blob location error distance for the RL parameter scan for MC
events. The first column shows the overall value while the second and third
columns show the blob1 and blob2 distribution respectively. The result for the
optimal configuration is marked in red. Top row: signal events within the
double escape peak ROI. Bottom row: background events in the same ROI.

As a cross-check and validation of the choice of main RL parameters, namely
the cleaning cut and number of iterations, the distance between the true
location of the track ends and the reconstructed ones was computed and
evaluated (for MC events), generating error distributions as shown in figure
15. To quantify the analysis, we consider the error level below which the
distribution contains 68% of the events. We scanned the main RL parameters
over a reasonable range of values, and calculated, for each configuration, the
68% error level (to which we refer as the “$1\sigma$ location error”).

The results of the scan are displayed on figure 19, which shows the $1\sigma$
location error obtained for different values of the number of iterations and
cleaning cut for both blobs (first column), blob1 (second column) and blob2
(third column), where the scan results are shown separately for signal (first
row) and background events (second row). Considering both blobs together
(first column), the optimal choice of $N_{iter}$ and $\epsilon_{cut}$ (red
rectangle) leads to an overall error close to the minimum, which would have
been obtained for $N_{iter}=75$ and $\epsilon_{cut}=0.011$. Note that the
blob1 population shows a much smaller change between configurations while the
gradient in blob2 population is considerably stronger, especially in the case
of background events. Figure 20 shows the end-location $1\sigma$ error values
of figure 19 as a function of the best $f.o.m.$ value obtained for each
configuration using the true information (where the best $f.o.m.$ is
calculated for the optimal choice of blob2 energy threshold). The data display
a clear correlation between high $f.o.m.$ values and small errors in end-
finding. Again, blob1 location error stays roughly constant while the blob2
location error improvement is fully correlated with increasing $f.o.m.$
values. This indicates that the overall improvement in the f.o.m. is primarily
driven by the error in the location of blob2.

SignalBackground

Figure 20: Blob location error distance, defined in the text, for all
configurations of the RL parameter scan and as a function of the best $f.o.m.$
value achieved with such configuration. Left: location error for signal events
within the double escape peak ROI. Right: background events in the same ROI.
Figure 21: Comparing different interpolation approaches for NEXT-White data
and MC (using the fit), for the optimal choice of RL parameters. (A) Signal
efficiency vs. background rejection for 1.6 MeV double escape peak events. (B)
The figure of merit for each dataset as a function of the blob2 cut.

Lastly, a brief study of the impact of the interpolation method was performed.
The evaluation consisted of repeating the topological analysis using the
optimal configuration but with different interpolation approaches, to compare
the figure of merit. In addition to the bicubic interpolation, we evaluated a
linear interpolation and a nearest-neighbor approach, were the value assigned
to each $10\times 10$ mm2 bin is equal to the number of photons detected in
the closest SiPM. The results are shown in figure 21. While mostly all
configurations are compatible within error, the maximal performance is
achieved with the bicubic interpolation. This was expected as the parameters
were optimized using bicubic interpolation and it is possible that different
results can be achieved with the other approaches if optimizing the parameters
for those. However, this study is left for future work.

##  Appendix C Examples of successful RL-based classification vs. failure of
the classical analysis

In this appendix we discuss three representative examples illustrating the
reasons for improved classification of signal and background events by the RL-
based method. The examples shown in figures 22 and 23 are simulated MC
background events in the 1.4-1.8 MeV ROI, which are misclassified by the
classical analysis as signal, and correctly identified as background by the RL
process. For each event, the top row shows the three projections of the raw
sensor response, binned in $10\times 10\times 1.8$ mm3 voxels, with a charge
cut of 30 PE for the classical analysis. The bottom row shows the
corresponding projections after RL deconvolution (charge cut of 10 PE with 75
iterations). Yellow symbols designate the track ends found by the classical
analysis and purple symbols - the ones found by the RL process. Squares
represent the center of blob2 and circles those of blob1. The deconvolved
images further include the true track overlaid in green. The circles represent
the blobs used in both analysis methods.

Event 1 is a photoelectric absorption of a Compton scattered gamma,
accompanied with a delta electron, where the track is “folded” such that its
true start point is close to its main part. Since the classical analysis
rebins the SiPM hits in 15 mm voxels, the true start point (purple square) is
merged into the track “body”, and a distant internal point (yellow square) is
misidentified as the track extremity. Since the local ionization density near
this point is high, the energy contained in the classical blob2 centered at it
lies above the threshold value, and the event is identified as signal. The RL
process, on the other hand, identifies both ends with an error of a few mm,
and correctly places the center of blob2 close to the track starting point.

Event 2 is a double Compton scatter of a 208Tl 2615 keV gamma. The first
scatter (starting point of the main track, close to the purple square) gives
rise to an energetic Compton electron that creates the main track. The
scattered photon interacts $\sim 2$ cm from the first vertex, generating a
Compton electron of lower energy and a low-energy gamma, which interacts at a
short distance. The classical analysis fails to separate the two Compton
electron tracks and places the center of blob2 close to the Bragg peak of the
second electron. Since the energy contained in this blob is above threshold
the event is identified as signal. RL deconvolution, in contrast, successfully
resolves the two tracks and places blob2 center close to the first Compton
vertex.

Event 1: binned SiPM

Event 1: deconvolution

Event 2: binned SiPM

Event 2: deconvolution

Figure 22: Single electron events (background) in the pair production ROI
misclassified by the classical analysis as signal and correctly identified as
background by the RL process. For each event the top row is the binned SiPM
response, and the bottom row – the deconvolved projections. Yellow symbols
mark the track ends found by the classical analysis, purple – by the RL
process. Squares are blob2 centers, circles – those of blob1. The true track
is overlaid on the deconvolved images. Misclassification of Event 1 results
from a merger between the start point of the track and its main part. Event 2
(double Compton scatter) is misidentified as a single track because of the
limited resolving power of the classical analysis.

Event 1: binned SiPM

Event 1: deconvolution

Figure 23: An additional background event misidentified as signal by the
classical analysis and correctly identified as background by the RL process.
Symbols have the same meaning as in figure 22. Here the mistaken
classification results from the high charge cut used in the classical analysis
which removes $\sim 2$ cm of the start of the track.

Event 3 (figure 23) provides an example for the effect of reduced SiPM charge
threshold ($q_{cut}$) in the RL process. The 30 PE cut employed in the
classical analysis (after careful optimization, as discussed in
Ferrario:2019kwg ) eliminates $\sim 2$ cm of the low ionization-density tail
of this event. By chance, the classical blob2 energy lies above threshold,
which leads to the misidentification of this event as signal. The same track,
under the RL process, is fully reconstructed in the tail region, with proper
placement of blob2 energy at the true starting point of the photoelectron.

## References

  * (1) F. T. Avignone, S. R. Elliott and J. Engel, _Double beta decay, Majorana neutrinos, and neutrino mass_ , _Rev. Mod. Phys._ 80 (2008) 481.
  * (2) S. Davidson, E. Nardi and Y. Nir, _Leptogenesis_ , _Phys. Rep._ 466 (2008) 105 .
  * (3) M. Blennow, E. Fernandez-Martinez, J. Lopez-Pavon and J. Menéndez, _Neutrinoless double beta decay in seesaw models_ , _JHEP_ 2010 (2010) 96.
  * (4) S. Dell’Oro, S. Marcocci, M. Viel and F. Vissani, _Neutrinoless double beta decay: 2015 review_ , _Adv. High Energy Phys._ 2016 (2016) 2162659.
  * (5) M. J. Dolinski, A. W. Poon and W. Rodejohann, _Neutrinoless double-beta decay: status and prospects_ , _Annu. Rev. Nucl. Part. Sci._ 69 (2019) 219.
  * (6) GERDA collaboration, _Final results of gerda on the search for neutrinoless double- $\beta$ decay_, _Phys. Rev. Lett._ 125 (2020) 252502.
  * (7) Majorana collaboration, _Search for neutrinoless double- $\beta$ decay in ${}^{76}\mathrm{Ge}$ with 26 kg yr of exposure from the Majorana Demonstrator_, _Phys. Rev. C_ 100 (2019) 025501.
  * (8) LEGEND collaboration, _The large enriched germanium experiment for neutrinoless double beta decay (LEGEND)_ , _AIP Conf Proc_ 1894 (2017) 020027.
  * (9) KamLAND-Zen collaboration, _Search for Majorana neutrinos near the inverted mass hierarchy region with KamLAND-Zen_ , _Phys. Rev. Lett._ 117 (2016) 082503.
  * (10) I. Shimizu and M. Chen, _Double beta decay experiments with loaded liquid scintillator_ , _Front. Phys._ 7 (2019) 33.
  * (11) EXO-200 collaboration, _Search for neutrinoless double-beta decay with the upgraded EXO-200 detector_ , _Phys. Rev. Lett._ 120 (2018) 072701.
  * (12) nEXO collaboration, _nEXO Pre-Conceptual Design Report_ , 1805.11142.
  * (13) DARWIN collaboration, _Sensitivity of the darwin observatory to the neutrinoless double beta decay of ${}^{136}\mathrm{Xe}$_, _EPJC_ 80 (2020) 808.
  * (14) J. J. Gómez-Cadenas, _Status and prospects of the NEXT experiment for neutrinoless double beta decay searches_ , in _54th Rencontres de Moriond on Electroweak Interactions and Unified Theories_ , pp. 201–206, 2019, 1906.01743.
  * (15) PandaX-III collaboration, _PandaX-III: Searching for neutrinoless double beta decay with high pressure ${}^{136}\mathrm{Xe}$ gas time projection chambers_, _Science China Physics, Mechanics & Astronomy_ 60 (2017) 061011.
  * (16) CUORE collaboration, _Improved limit on neutrinoless double-beta decay in ${}^{130}\mathrm{Te}$ with CUORE_, _Phys. Rev. Lett._ 124 (2020) 122501.
  * (17) SNO+ collaboration, _Current status and future prospects of the SNO+ experiment_ , _Adv. High Energy Phys._ 2016 (2016) 6194250.
  * (18) CUPID collaboration, _CUPID pre-CDR_ , 1907.09376.
  * (19) CUPID-Mo collaboration, _New limit for neutrinoless double-beta decay of ${}^{100}\mathrm{Mo}$ from the cupid-mo experiment_, _Phys. Rev. Lett._ 126 (2021) 181802.
  * (20) AMORE collaboration, _First results from the AMoRE-Pilot neutrinoless double beta decay experiment_ , _EPJC_ 79 (2019) 791.
  * (21) D. Nygren, _High-pressure xenon gas electroluminescent TPC for $0\nu\beta\beta$-decay search_, _Nucl. Instrum. Meth. A_ 603 (2009) 337.
  * (22) A. Bolotnikov and B. Ramsey, _The spectroscopic properties of high-pressure xenon_ , _Nucl. Instrum. Meth. A_ 396 (1997) 360 .
  * (23) NEXT collaboration, _Energy calibration of the NEXT-White detector with 1% resolution near Q ββ of ${}^{136}\mathrm{Xe}$_, _JHEP_ 10 (2019) 230.
  * (24) XENON collaboration, _Energy resolution and linearity of XENON1T in the MeV energy range_ , _EPJC_ 80 (2020) 785.
  * (25) R. Luescher, J. Farine, F. Boehm, J. Busto, K. Gabathuler, G. Gervasio et al., _Search for $\beta\beta$ decay in ${}^{136}\mathrm{Xe}$: new results from the Gotthard experiment_, _Phys. Lett. B_ 434 (1998) 407 .
  * (26) NEXT collaboration, _The Next White (NEW) detector_ , _JINST_ 13 (2018) P12010.
  * (27) NEXT collaboration, _Sensitivity of a tonne-scale NEXT detector for neutrinoless double beta decay searches_ , 2005.06467.
  * (28) B. Jones, A. McDonald and D. Nygren, _Single molecule fluorescence imaging as a technique for barium tagging in neutrinoless double beta decay_ , _JINST_ 11 (2016) P12011.
  * (29) NEXT collaboration, _Demonstration of single-barium-ion sensitivity for neutrinoless double-beta decay using single-molecule fluorescence imaging_ , _Phys. Rev. Lett._ 120 (2018) 132504.
  * (30) N. Byrnes, F. W. Foss, B. Jones, A. McDonald, D. Nygren, P. Thapa et al., _Progress toward barium tagging in high pressure xenon gas with single molecule fluorescence imaging_ , _J. Phys: Conf. Ser._ 1312 (2019) 012001.
  * (31) P. Thapa, I. Arnquist, N. Byrnes, A. A. Denisenko, F. W. Foss, B. J. P. Jones et al., _Barium chemosensors with dry-phase fluorescence for neutrinoless double beta decay_ , _Sci. Rep._ 9 (2019) 15097.
  * (32) I. Rivilla, B. Aparicio, J. M. Bueno, D. Casanova, C. Tonnelé, Z. Freixa et al., _Fluorescent bicolour sensor for low-background neutrinoless double $\beta$ decay experiments_, _Nature_ 583 (2020) 48.
  * (33) NEXT collaboration, _Sensitivity of NEXT-100 to neutrinoless double beta decay_ , _JHEP_ 2016 (2016) 159.
  * (34) NEXT collaboration, _First proof of topological signature in the high pressure xenon gas TPC with electroluminescence amplification for the NEXT experiment_ , _JHEP_ 01 (2016) 104.
  * (35) NEXT collaboration, _Demonstration of the event identification capabilities of the NEXT-White detector_ , _JHEP_ 10 (2019) 052.
  * (36) NEXT collaboration, _Demonstration of background rejection using deep convolutional neural networks in the next experiment_ , _JHEP_ 2021 (2021) 189.
  * (37) W. H. Richardson, _Bayesian-based iterative method of image restoration_ , _J. Opt. Soc. Am._ 62 (1972) 55.
  * (38) L. B. Lucy, _An iterative technique for the rectification of observed distributions_ , _Astron. J._ 79 (1974) 745.
  * (39) E. D. C. Freitas, C. M. B. Monteiro, M. Ball, J. J. Gómez-Cadenas, J. A. M. Lopes, T. Lux et al., _Secondary scintillation yield in high-pressure xenon gas for neutrinoless double beta decay (0 $\nu$$\beta$$\beta$) search_, _Phys. Lett. B_ 684 (2010) 205.
  * (40) NEXT collaboration, _Calibration of the NEXT-White detector using ${}^{83m}\mathrm{Kr}$ decays_, _JINST_ 13 (2018) P10014.
  * (41) T. Cormen, C. Stein, R. R.L. and C. Leiserson, _Introduction to algorithms, 2nd ed_. McGraw-Hill Higher Education, U.S.A., 2001.
  * (42) J. Martín-Albo, _The NEXT experiment for neutrinoless double beta decay searches_ , Ph.D. thesis, Valencia University, IFIC, Valencia, Spain, 2015\.
  * (43) NEXT collaboration, _Electron drift properties in high pressure gaseous xenon_ , _JINST_ 13 (2018) P07013.
  * (44) S. van der Walt, J. L. Schönberger, J. Nunez-Iglesias, F. Boulogne, J. D. Warner, N. Yager et al., _scikit-image: image processing in Python_ , _PeerJ_ 2 (2014) e453.
  * (45) NEXT collaboration, _Helium-xenon mixtures to improve the topological signature in high pressure gas xenon TPCs_ , _Nucl. Instrum. Meth. A_ 905 (2018) 82 .
  * (46) NEXT collaboration, _Electron drift and longitudinal diffusion in high pressure xenon-helium gas mixtures_ , _JINST_ 14 (2019) P08009.
  * (47) NEXT collaboration, _Low-diffusion Xe-He gas mixtures for rare-event detection: electroluminescence yield_ , _JHEP_ 2020 (2020) 34.
  * (48) C. Azevedo, L. Fernandes, E. Freitas, D. Gonzalez-Diaz, F. Monrabal, C. Monteiro et al., _An homeopathic cure to pure xenon large diffusion_ , _JINST_ 11 (2016) C02007.
  * (49) NEXT collaboration, _Secondary scintillation yield of xenon with sub-percent levels of CO2 additive for rare-event detection_ , _Phys. Lett. B_ 773 (2017) 663 .
  * (50) NEXT collaboration, _Electroluminescence TPCs at the thermal diffusion limit_ , _JHEP_ 2019 (2019) 27.
  * (51) NEXT collaboration, _Application and performance of an ML-EM algorithm in NEXT_ , _JINST_ 12 (2017) P08009.