## Abstract

We introduce a neutron-gamma emission tomography (NGET) technique for rapid detection, three-dimensional imaging, and characterization of special nuclear materials like weapons-grade plutonium and uranium. The technique is adapted from fundamental nuclear physics research and represents a previously unexplored approach to the detection and imaging of small quantities of these materials. The method is demonstrated on a radiation portal monitor prototype system based on fast organic scintillators, measuring the characteristic fast time and energy correlations between particles emitted in nuclear fission processes. The use of these correlations in real time in conjunction with modern machine learning techniques provides unprecedented imaging efficiency and high spatial resolution. This imaging modality addresses global security threats from terrorism and the proliferation of nuclear weapons. It also provides enhanced capabilities for addressing different nuclear accident scenarios and for environmental radiological surveying.

## INTRODUCTION

While special nuclear materials (SNMs) like plutonium and uranium form the basis for enormous technological developments in the last century, they have also led to some of the most serious threats to humankind and its very existence. Although the danger of a global military conflict involving nuclear weapons of mass destruction is at the center of public awareness in this context, the potential consequences from acts of terrorism involving these materials may not fall far behind. The stewardship and control of SNM is therefore of major importance for protecting the safety of present and future generations, and it is one of the main concerns of international collaboration within the International Atomic Energy Agency (IAEA). Key to these efforts is the ability to efficiently and quickly detect, localize, and characterize SNM in different national security scenarios as well as in routine inspection of passengers, goods, vehicles, etc. at, e.g., sea- and airports and transport logistics hubs. This requires advanced sensor-based detection of threat objects, including radiation detection systems. Among the most important materials of interest are weapons-grade uranium, weapons-grade plutonium (WGP), and radioactive sources that could be used for radiological dispersion devices. An effective defense against these and similar threats must take advantage of the latest developments in radiation detection technology. Much-needed advances in this field include not only enhanced sensitivity of detection but also the possibility to perform efficient localization of SNM in screening operations of, e.g., objects, persons, and vehicles in the field. Besides detection systems related to national security, there are several other areas that would benefit greatly from the possibility to efficiently image SNM. These include nuclear emergency response, nuclear safeguards, and environmental surveying.

The detection of SNM in radiation portal monitors (RPMs) is one of the critical links in the nuclear security chain. RPMs are designed to alarm if the measured radiation flux, primarily of neutrons and γ-rays, exceeds predefined thresholds when persons, vehicles, packages, or other objects are passing through them. Once a source is detected, it would be highly beneficial to rapidly discern its precise location. Current RPM systems do not have these capabilities. Nuclear emergency responders and safeguards inspectors would also benefit substantially from the ability to rapidly verify the location of a source without close manual intervention.

Because the presence of neutrons is a key indicator of materials undergoing spontaneous or induced fission, neutron detectors are crucial for the detection of SNM. Thermal neutron gaseous proportional counters based on ^{3}He are currently the most common detectors used in nuclear safeguards and nuclear security systems, such as in RPMs (*1*, *2*). However, because of the global shortage of ^{3}He, there is an increasing interest in developing ^{3}He-free systems (*3*). In addition to proportional counter systems replacing ^{3}He with other high thermal-neutron absorption cross-section materials, there is also an increasing focus on fast-neutron detection based on organic scintillators. While thermal neutrons typically do not carry any useful information of their initial energy or emission point, fast-neutron detection (i.e., the detection of neutrons with typical kinetic energies in the megaelectronvolt range that are unperturbed on their way from the source) opens new possibilities for imaging SNM. Several types of neutron imagers have been developed, including coded-aperture imagers (*4*, *5*), time-encoded imagers (*6*), and neutron scatter cameras (*7*–*13*). While neutron scatter cameras can be made more compact than time-encoded or coded-aperture neutron imagers, their efficiency suffers greatly from the requirement of detecting two consecutive neutron scatters. Furthermore, a fundamental complication in the design of neutron scatter cameras is the need to maximize the probability that the incident neutron scatters only once in the scatter detector and subsequently escapes so that it may interact in another detector element. If the scatter detector is too thick, additional scatters become likely, corrupting the energy measurement and making it impossible to reconstruct the neutron incident direction. Attempts to circumvent this problem would lead to reduced detection efficiency or to complex and costly detection systems.

Although the main focus of previous studies has been on neutron detection, there are potential advantages associated with the additional detection of γ-rays as a means to identify the presence of SNMs, with very high sensitivity. These are mainly the significantly higher multiplicity of prompt fission γ-rays and their short flight time from source to detector. Most of the γ-rays in fission are emitted in prompt cascades depopulating short-lived (typically picoseconds or shorter) excited states in the fission products, with a multiplicity distribution that can extend significantly beyond an average of 5 to 10 (*14*–*16*). Therefore, in particular, fast γ-neutron correlations deserve further attention in the development of more sensitive methods for nuclear safeguards and security applications (*17*). These techniques have been used in fundamental nuclear physics experiments for decades [see (*18*) for a recent example]. An additional advantage of fast γ-neutron coincidence detection compared with detection of neutrons alone is the energy independence of the photon time of flight, which reduces the correlation time spread between the detected particles. It was previously proposed to use time-of-flight (TOF) information for event-by-event three-dimensional (3D) imaging of SNM in neutron scatter cameras (*19*) by combining the fast coincidence detection of a neutron scattering between two detector elements with the detection of a correlated photon in a third detector element. However, for these triple-coincidence systems, the detection efficiency will be notably reduced compared with the already low efficiency of neutron scatter imaging and it may not be useful for applications where a rapid assessment is needed. The velocity distributions of neutrons emitted in the fission processes are well characterized for the most common nuclides present in SNM. This information is used in the present work to deduce the position of a source from measured fast γ-neutron TOF correlations.

We here present a new approach, neutron-gamma emission tomography (NGET), to 3D localization of SNM using fast γ-neutron coincidence detection applicable to national security measures, nuclear emergency response, nuclear safeguards, environmental radiological surveying, and related fields. While the 3D imaging method is here demonstrated within the framework of an RPM prototype system, it is versatile and can be adapted to a wide variety of different detection geometries and applications.

Two different techniques are presented, both exploiting the properties of correlated γ-neutron pairs originating from the same fission events. We first show results from an implementation of this approach based on estimating a region of origin from the measured energy deposited by the neutron and the γ-neutron time difference, event by event. Only the detection of a single neutron interaction is required, and the neutron energy is estimated on the basis of sampling the approximately known kinetic energy distribution of fission neutrons above the detected recoil energy, which is the minimum kinetic energy that could be carried by the incident neutron. The method has similarities with emission tomographic methods used in medical imaging, such as positron emission tomography (PET) (*20*, *21*). PET imaging uses the fact that 511-keV photon pairs from positron annihilation are strongly correlated, directionally in space and in time, to deduce the distribution of positron-emitting isotopes within the field of view. A PET detector system with good enough (∼subnanosecond) time resolution can, additionally, exploit the relative TOF information for the detected photon pairs to improve on its 3D image resolution (*22*). Differently from positron annihilation, the physics of the nuclear fission process does not provide easily deduced direct directional correlations between the emitted photons and neutrons. However, neutrons being massive particles have a velocity that is directly related to their kinetic energy. Furthermore, because neutrons emitted from spontaneous fission of SNM have rather well-established energy distributions, often approximated by a Watt spectrum (*23*) with parameters depending weakly on the nuclide, it is possible to estimate the probability that a detected neutron had a certain incident velocity based only on a partial energy measurement or without measuring the energy at all. Using these estimates, the measured relative TOF between the neutron and the photon in a correlated neutron-photon pair can then be translated into information about their point of origin, i.e., the position of the radiation source. To our knowledge, it is the first time a physical model of the emission process is directly applied to tomographic imaging. Similarly to a neutron scatter camera, this technique then relies on standard image reconstruction techniques. We here present results obtained using a deconvolution algorithm based on Bayes’ theorem (*24*).

In another implementation of NGET imaging, we rely on cumulative distributions of correlated γ-neutron time differences recorded between different pairs of radiation sensors for rapid, online reconstruction of the source position. A particular advantage of this cumulative NGET (CNGET) technique is that it only requires the difference in detection time between correlated particles. Hence, it is enough that the detector elements have a fast (nanosecond) time response and that they are sensitive to both γ-rays and fast neutrons. The CNGET technique therefore does not, a priori, require the usually more expensive ability to discriminate between neutrons and γ-rays, nor does it require measurement of particle energies. This technique relies on machine learning algorithms or statistical iterative methods to accurately localize SNM.

By localizing unknown sources in 3D, both implementations of the NGET imaging method also provide the basis for identifying, quantifying, and characterizing small amounts of SNM.

## RESULTS

In this work, rapid and accurate localization of SNM samples with NGET is demonstrated using an RPM prototype system developed at the Royal Institute of Technology (KTH) (*25*). A basic requirement on the detection system is an array of radiation sensors with fast time response, sensitive to fast neutrons and γ-rays. Most favorably, the sensors are simultaneously sensitive to both fast neutrons and γ-rays and have a high capability to discriminate between the two types of particles. This is the case for the organic scintillation detectors used in the present work, which provide efficient γ-neutron discrimination and fast timing of the order of 1 ns. The RPM prototype system consisted of an array of eight 127-mm-diameter by 127-mm-length cylindrical liquid organic scintillator cells arranged in two detector assemblies, displaced by 1 m from each other (see the ‘Discussion’ section). Measurements have been carried out using a californium-252 (Cf-252) radioactive source with mass 3.2 × 10^{−9} g encapsulated in a 4.6-mm-diameter × 6-mm cylindrical ceramic casing. Cf-252 has spontaneous fission as a 3% decay branch (the remaining 97% being due to α decay) and an average prompt-neutron multiplicity of 3.76 per fission (*26*). The source is equivalent to around 100 g of WGP (7% plutonium-240 and 93% plutonium-241) in terms of its approximately 1900 s^{−1} fission rate. This would correspond to a sphere of radius approximately 1-cm metallic plutonium. Further details of the setup and the measurements are described in Materials and Methods.

### NGET imaging performance

We first present results from applying the novel NGET imaging principle with the KTH RPM prototype using the technique based on event-by-event mapping of possible points of origin using probability density functions (PDFs). The PDFs are derived from the measured neutron recoil energy and relative detection time of γ-neutron pairs, which are correlated from the same fission event. The principle is illustrated schematically in Fig. 1. For each correlated detection of a neutron and a γ-ray in two specific detector elements, the neutron and γ-ray flight paths are well defined as a function of the possible position of a fission event in space. Knowledge about the prompt fission neutron energy spectrum shape can then be used to extract the probability distribution for the unknown source position given the observation of a certain time difference between γ-ray and neutron detection. This probability distribution can be further constrained using information on the detected neutron recoil energy. The measured neutron-recoil energy is, in most cases, due to one or more elastic scatterings on protons in the detector medium. Neutrons and protons have approximately the same mass, so the scattering process is similar to classical billiards. For each energy deposition in the detector material from a recoiling proton, the incident neutron energy must be higher or equal to this value. The detector energy response to a single neutron interaction is an approximately known nonlinear function of the deposited energy. The full detector response as a function of incident neutron energy, which is highly dependent on the incident neutron energy and the geometrical dimensions of the detector, can also be calculated approximately or measured. Hence, known energy distributions of fission neutrons can, in combination with the detector response function, be used to estimate the probability distribution for the incident neutron energy. Each γ-neutron coincidence event then gives rise to a unique PDF, depending on which detector elements were involved, the measured time difference between the γ-ray and neutron interactions, and the detected neutron energy. By combining the information from several such events, a combined PDF for the source location can be deduced. In this way, the uncertainty of the source location is successively reduced for each correlated γ-neutron event that is detected. By combining information from several, geometrically displaced, sensor pairs, a 3D localization in space is achieved.

The algorithm developed for this purpose uses Bayesian inference and is described in detail in Materials and Methods. The result of the algorithm is the PDF of finding the source as a function of the position. The application of the algorithm to experimental data from the KTH RPM prototype system with the source in different positions is illustrated in Fig. 2. After 10 s of measurement, corresponding to between 100 and 200 collected events, the source location PDF has an SD of about 4 cm. A longer measurement time of 30 s reduces the SD by about a factor of 2. For longer measurement times, the average radial SD of the point spread function that is achieved with this method for the current detector configuration is around 1 cm.

The Bayesian inference algorithm for source localization assumes that all detected events originate from the same point and it can therefore not simultaneously identify the positions of multiple sources. For this purpose, a different algorithm is applied, which has a somewhat lower spatial resolution for a single source but is capable of imaging multiple sources. The application of this algorithm to the measurement of two Cf-252 sources, one with approximately three times lower activity, is shown in Fig. 3.

### NGET imaging based on cumulative time difference distributions (CNGET)

CNGET 3D imaging is a simplified approach to full NGET providing rapid and accurate 3D localization of SNM based solely on measuring the accumulated γ-neutron arrival time difference distributions for different combinations of detector elements. In the present detector system, there are eight different sensor elements and therefore 8 × 7 = 56 unique combinations of γ-neutron time difference distributions. Assuming a point source of time-correlated neutrons and γ-rays such as a small amount of SNM, each time distribution contains unique information regarding its position. For many detection geometries, such as the present one, it may be assumed that the photon is detected before the neutron because the mean velocity of neutrons from nuclear fission is roughly an order of magnitude smaller than the speed of light. Therefore, the method does not require γ-neutron discrimination capabilities, nor does it require measuring the particle energies. The CNGET technique uses the full set of time difference distributions, which are accumulated continuously during the measurement time, each corresponding to one unique pair of detector elements to localize the source. The measurement time will vary depending on the application, from short (seconds to minutes) time scales in the case of security and emergency scenarios to longer time scales, e.g., for spent fuel verification in nuclear safeguards. Figure 4 shows examples of γ-neutron time difference distributions obtained from measuring on the same Cf-252 radioactive source as described above, placed in different positions. As can be seen in the figure, these distributions have different characteristics depending on the position of the source. Hence, they can be used to provide the location of the source in space with the aid of, e.g., machine learning.

In the present study, we have trained an artificial neural network (ANN) algorithm using the deep learning toolbox in the MATLAB (*27*) computational software package. Figure 5 shows the results obtained for the Cf-252 source described above when placed in three different positions, each point measured during 1 min.

## DISCUSSION

To ensure public safety and international security, the risk that nuclear or other radioactive materials could be subject to illegal trafficking and used in criminal or other unauthorized acts must be minimized. This requires versatile detection systems that can rapidly detect SNM and other radiological threats in real time. The possibility to rapidly and precisely locate SNM, e.g., at border controls, postal mail centers, or transport hubs, would provide an important tool for counterterrorism measures and nonproliferation efforts. 3D imaging of SNM also has several other important potential applications, e.g., for nuclear emergency response, spent fuel verification in nuclear safeguards, and environmental surveying.

In the recent work of Steinberger *et al.* (*13*), a state-of-the-art compact 2D radiation imaging system for SNM was presented on the basis of the principle of neutron scatter imaging with organic scintillation detectors. The here presented NGET and CNGET techniques enable localization of SNM in 3D with higher accuracy and orders of magnitude higher efficiency. Although there is a substantial difference between the detector system used in the present work and that used by Steinberger *et al.* in terms of detection geometry, in the total volume of the active detector medium, etc., an important reason for the higher performance obtained here is the use of the fundamental properties of γ-fast neutron correlations from nuclear fission, in particular the emission of multiple, prompt γ-rays with a high average multiplicity (*25*).

It is noteworthy that the RPM detector system used in the present work has not been optimized for efficiency, nor has it been designed for imaging purposes. Substantially improved spatial resolution and response uniformity could be obtained by using a set of more uniformly distributed detectors or adopting a higher detector granularity (i.e., smaller detector cells). Still, the results are remarkably good with an average position uncertainty of 4.2 cm for a 10-s measurement of this relatively weak source, approximately 35% of the Cf-252 source activity specified in the ANSI N42.35-2016 industry standard (*28*) for RPM sensitivity characterization with sources moving at a speed of 1.2 m/s. Furthermore, while the imaging performance is here demonstrated with stationary sources, the method could readily be adapted to linearly moving objects with the aid of an optical tracking system.

In this work, we have discussed the use of γ-fast neutron correlations from fission for rapid localization and imaging of SNM. An additional source of correlated neutrons and γ-rays from these materials are (α, nγ) reactions on light elements [see (*29*)]. These reactions are possible when actinides are in close contact with light elements, e.g., as fluoride or oxide compounds commonly used in the processing and fabrication of nuclear fuels. Many of the most important actinides present in common SNM undergo spontaneous α-decay. The emitted α-particles may induce (α,n) reactions with light nuclei when their kinetic energy is sufficient to overcome the Coulomb barrier. In a fraction of these reactions, the residual nucleus is left in an excited state, typically with a short, subnanosecond lifetime, which may decay with the subsequent emission of one or more γ-rays. The time correlations between the emitted neutrons and γ-rays are typically similar to those for spontaneous fission, although the multiplicity of neutrons and, in particular, γ-rays in these reactions are lower. Although the image reconstruction efficiency is expected to be lower for these reactions, they may provide an important means to locate the presence of some nuclides of interest for which the spontaneous fission branch is insignificant. The neutron energy spectra from (α,n) reactions at the relevant incident α-particle energies (*30*) can have substantial overlap with those from spontaneous fission. Nevertheless, the effects of (α, nγ) reactions on the performance of the NGET imaging algorithm and its optimization require further investigation.

Scattering in the object of interest or in surrounding materials is a potential problem for all radiation imaging methods. One of the challenges in applications, such as nuclear security, safeguards, and emergency response, is the possible presence of shielding materials that may not only attenuate the radiation from the source but also distort the time and energy correlations in the detected radiation. For example, RPMs may need to deal with radiation sources that are embedded in cargo and/or clandestinely transported within shielding materials. Materials surrounding a source reduce the measurable flux, thereby increasing the minimum detectable activity. They may also lead to bias in the localization due to inhomogeneous attenuation and scattering of neutrons resulting in different energies and flight times. Similar effects are present for neutron scatter cameras and warrant further attention.

Both NGET techniques presented here can be used to rapidly localize small amounts of SNM with high accuracy. This opens the possibility for a paradigm shift in national security applications, nuclear safeguards, and environmental surveying and for nuclear emergency responders. In addition, the CNGET technique only relies on the fast timing properties of the sensor elements and can therefore be implemented with relatively modest technical modifications into standardized RPM designs using conventional plastic scintillators. On the other hand, the image resolution improves more rapidly as a function of measurement statistics (i.e., measurement time) with full event-by-event NGET imaging due to the additional information on the neutron kinetic energy that is included in the image reconstruction. An additional advantage is that the event mapping technique uses an analytical probabilistic approach to the reconstruction of the source location as opposed to empirical machine learning methods.

## MATERIALS AND METHODS

### Detection system

The 3D imaging capabilities of NGET are developed and illustrated within the framework of an RPM system. However, the method is applicable to any radiation detector system containing a minimum of two geometrically displaced sensor elements with a fast (ns) time response and with at least one detector element sensitive to fast neutrons and at least one sensor element sensitive to γ-rays. Some organic scintillators are well suited for this purpose as they are sensitive to and may discriminate between fast neutrons and γ-rays using pulse shape analysis (PSA). They also have a fast time response (typically around 1 ns or better) as well as limited spectroscopic capabilities. The RPM prototype system constructed at KTH (*25*) that we use for illustrating the new imaging method consists of an array of sensors placed in two vertical pillars with a geometry as specified in the ANSI N42.35-2016 industry standard (*28*), together with the associated electronics. In the current configuration, a total of eight 127-mm-diameter by 127-mm-length cylindrical EJ-309 scintillator (*31*) cells are arranged in two two-by-two detection assemblies, one in each pillar. A computer-aided design (CAD) model of the detector system including its mechanical support structure is shown in Fig. 6. The EJ-309 scintillator was chosen for its excellent γ-neutron discrimination properties (*32*), which enabled the RPM to detect single neutrons with a low level of contamination from background γ-rays (*25*). γ-Neutron discrimination is not strictly needed for imaging because correlated neutrons and γ-rays are typically well separated by their relative detection times. However, the results are somewhat improved by identifying the neutrons using PSA and measuring their recoil energy in event-by-event mode.

It is worthwhile to point out that the imaging capabilities of the system are not dependent on the object being inside the portal, nor that the system contains two separate pillars of sensor elements. Each detection assembly has its own inherent imaging capabilities, and different, more compact sensor geometries are also possible, providing single-view imaging. The attainable image resolution will depend on the detailed geometry of the detection system.

The EJ309 scintillator was chosen as detection medium as it features a number of chemical properties recommending it for use in environmentally challenging conditions outside the laboratory. These properties include a high flash point, low vapor pressure, and low chemical toxicity (*31*). However, there are other organic scintillator materials with excellent and complementary properties like trans-Stilbene (trans-1,2-diphenylethene) or plastic scintillators (*33*, *34*) available that may be equally well or better suited for different applications.

Each detector cell is optically coupled to a Hamamatsu R1250 photomultiplier tube. The four photomultiplier tubes in each detection assembly are powered by a four-channel CAEN DT5533 high-voltage power supply. The photomultiplier tube anode pulses are acquired and digitized with an eight-channel CAEN DT5730 digitizer board featuring a 2-Vpp dynamic range, 14-bit resolution, and 500-MHz sampling rate. The time-synchronized data are transferred using a universal serial bus (USB) 2.0 link to a personal computer (PC) running on a Linux or Windows platform. Besides the raw data transfer of detector waveforms, the digitizer also features firmware programmable PSA using two field programmable gate arrays (FPGAs) of type Intel/Altera Cyclone EP4CE30, each handling four of the analog-to-digital converter data streams. The digitized signals from the RPM detector modules can be processed in real time within the digitizer FPGAs to extract charge integrals, pulse shape information, and timing information. The data are transmitted via USB to an industrial PC, where further processing is performed. The PSA algorithm for distinguishing γ-rays from neutrons is based on the charge comparison method (*35*), whereas the time-of-arrival information is extracted using a digital constant fraction discrimination algorithm. Energy information for each trace is deduced using a moving window deconvolution algorithm. The processed data stream contains information on the individual detected neutrons and γ-rays (energy, time, and the sensor element that registered the hit) as well as time-averaged single-neutron and γ rates, fast coincidence rates for γ-neutron, neutron-neutron, and γ-γ events with adjustable coincidence time windows and thresholds. Typical coincidence time windows are 0 to 100 ns for γ-γ and neutron-neutron events and 10 to 100 ns for γ-neutron events. For different detection systems, the optimal coincidence time windows will scale with the dimensions of the objects of interest.

### Methods

*Source localization using Bayesian inference*. The NGET algorithm designed to find the location of a fission source is based on Bayes’ theorem, which gives a way to estimate the probability of a hypothesis *H* to be true, given the occurrence of an event *E**P*(*E*∣*H*) is the probability of the event to occur given that the hypothesis is true, *P*(*H*) is the probability that the hypothesis is true, and *P*(*E*) is the total probability that the event occurs.

We define the event *E* as the detection of a γ-neutron coincidence, with a time difference Δ*t*_{nγ} between γ-ray and neutron detection and a detected neutron-recoil energy, as measured from the light output *L*. For each voxel in space, we form the hypothesis *H _{i}* that the event originated from the center of the voxel

From the neutron TOF, its kinetic energy, *E _{n}*, is calculated. The probability that a neutron with this kinetic energy is emitted is taken from evaluated prompt fission neutron spectral data (

*36*). The PDF of possible event origins forms a fuzzy semispherical shell around the neutron detector position in space (see the illustration in Fig. 7).

The predominant interaction of neutrons in the organic scintillator is due to elastic scattering off protons. The neutron mean free path in the detector medium is on the order of a few centimeters for typical fission neutron energies. Therefore, the probability to have more than one scatter within a detector element can be relatively large, even for modest detection volumes. Each scattering interaction transfers a portion of the neutron’s kinetic energy to the recoiling proton. A fraction of this energy is converted into fluorescent light in the scintillator. In the measurements, the scintillation light is converted to an electrical current pulse in the photomultiplier tube. The charge contained in this pulse is then integrated to provide a measure of the deposited energy.

Qualitatively, the light output recorded from neutron interactions in the detector medium allows us to determine a lower limit for the kinetic energy of the incident neutron. However, because of the nonlinear conversion of recoil energy into scintillation light and the stochastic nature of the scattering process, there is a large variation in light output from the scintillator for events with a given incident neutron energy. A quantitative estimate of the probability that a neutron of a given kinetic energy produces a certain light output is given by the light-output response function *R*(*L*∣*E _{n}*). We have calculated it using the Monte-Carlo code Geant4 (

*37*) to account for the possibility of several scattering events inside the detector by the same neutron as well as other types of interactions. In the Monte-Carlo code, we use the proton light-output function from Enqvist

*et al.*(

*38*) and also take the detector energy resolution into account. An example of the light-output response function is shown in fig. S1, where it is compared with experimental data from Enqvist

*et al.*(

*38*). This information is then used when calculating the PDF for each event.

Given an event *E*, the probability of the hypothesis *H _{i}* to be true is calculated as

*E*) is the prompt fission neutron spectrum, δ

_{n}*E*is the neutron kinetic energy resolution, and the index

_{n}*i*refers to a specific voxel. The total probability that the event occurs is then

The probability that the hypothesis is true is initially unknown and the quantity that we are looking for. As an initial guess, it is taken to be uniformly distributed in the space defined by the voxels*n* is the number of voxels. The probability that the hypothesis is true is updated after each event. If we would be sure that all observed events are from the source that we are looking for, it would be natural to use*k* indexes the event number. In this case, the convergence is very fast; however, any spurious event due to accidental coincidence with background or scattered neutrons has a detrimental effect on the final result. To make the convergence more reliable, we therefore update the image according to*P*_{k − 1}(*E*) in the weighted sum; this quantity is small for events whose PDFs have a small overlap with the current guess for the source location probability. Thereby, it gives a way to suppress events that do not agree with most events.

The result of the algorithm is the PDF of finding the source as a function of the position. The application of the algorithm to experimental data from the RPM with the source in different positions is illustrated in Fig. 2. After 10 s of measurement with the detection geometry and source used in the experiments (1900 fissions/s), the source location PDF has an SD of about 4 cm. The absolute accuracy of the most probable location for long measurement times with the present detector geometry is given by a point spread function with an average SD of about 1 cm in a sensitive volume of 1.0 m × 1.0 m × 1.8 m. The achievable accuracy is highly dependent on the number of radiation sensors and their geometrical placement. A larger area of sensitivity and improved accuracy could be achieved with a refined detection geometry.

*Imaging multiple and distributed sources*. The basic assumption made in the source localization algorithm is that all events originate from the same source position. If more than one source is present, the algorithm will, at best, find only the most intense one. However, image reconstruction can be applied in a similar way to the Bayesian source localization algorithm. Instead of treating each event separately, accumulated data on Δ*t*_{nγ} and *L* is used. From the accumulated data, an image is deconvoluted using the iterative Bayesian unfolding algorithm of D’Agostini (*24*). The number of coincidences originating from the voxel *i* is estimated as*n _{E}* is the number of events. The response function

*P*(

*E*∣

*H*) is the same as in the source localization algorithm described above. The deconvolution algorithm requires a first guess of the image

_{i}*P*(

*H*) as starting point. The closer the initial guess is to the true image, the better the result will be. In Fig. 3, we demonstrate the functionality of the method in the worst-case scenario, i.e., no prior knowledge of the image. As starting point, the prior is then taken to be uniformly distributed in space. The image intensity resulting from Eq. 10 is used as the prior in the next iteration and so on.

When applying an iterative image deconvolution, it is necessary to use a convergence criterion to determine when to stop the iteration. We use here the same approach as Steinberger *et al.* (*13*). The variance in image difference from two consecutive iterations is used to characterize how nonuniformly the image changed*n _{C}*(

*H*)}

_{i}*is the vector of voxel intensities for iteration number*

_{k}*k*. In fig. S2, this parameter is plotted as a function of the iteration number, for the data presented in Fig. 3. We have chosen the convergence criterion

*NGET imaging based on cumulative time difference distributions (CNGET)*. For a system of *N* detector elements, there are *N*(*N* − 1) unique time difference distributions, taking into account which detector element detected the first in a time-correlated pair of particles emitted from a fission event. The CNGET technique uses this complete set of cumulative time difference distributions to determine the location of SNM in 3D using machine learning. With this technique, one may, similarly to event-by-event–based NGET, analyze these time difference distributions continuously as they are updated during a measurement, successively improving on the accuracy of the localization. Although not discussed in the present work, the method can also be applied to imaging of multiple or distributed sources using iterative image reconstruction methods.

Machine learning based on a four-layer feed-forward ANN was used to analyze the data in the present work. The Neural Network Fitting App (nftool) from the Deep Learning Toolbox in the MATLAB computational software was used for training and testing. Nftool can solve the data-fitting problems using feed-forward networks with two or more layers. From this app, a MATLAB script is generated to modify the training process. In the present work, training data were constructed as a matrix of input column vectors. Each vector included the first, second, third, and fourth moment extracted from the γ-neutron time difference distributions, for all 56 detector combinations. The target vectors were the actual positions of the Cf-252 source for each measurement. During the training process, the training and test data were randomly divided into two sets so that 90% of the data were used for training and 10% were used as a completely independent test of network generalization. The ANN used a sigmoid transfer function in the three hidden layers and a linear transfer function in the output layer. The number of hidden neurons was set to 10 in each of the three hidden layers. The best prediction results were achieved when the network was trained using a Bayesian regularization algorithm. During the training process, the performance of the network was monitored. The regression plot, displayed in fig. S3, displays the network outputs for training and test data. If the network outputs would be equal to the targets, the data would lie along a 45° line for the perfect fit and the *R* value would be equal to 1.0. For our problem, the fit was objectively good, with *R* values higher than 0.92 for all sets. The *R* value for the test set will, after some training time, reach its highest value and slowly start to decrease due to overfitting. This is the sign to stop the training process.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/7/21/eabg3032/DC1

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**We thank D. Chernikova for valuable support.

**Funding:**This work has received funding from the Swedish Radiation Safety Authority under contract nos. SSM2016-3954, SSM2018-4393, and SSM2018-4972 and the Carl Trygger Foundation for Scientific Research under contract no. CTS14:89.

**Author contributions:**B.C. conceived and supervised the project. J.P. and A.G. calibrated the detectors and carried out the measurements. J.P. developed the ANN and carried out the CNGET data analysis. A.G. developed the source localization algorithms based on Bayesian Inference and carried out the event-by-event data analysis. All the authors participated in discussions of the results and contributed to the writing of the manuscript.

**Competing interests:**B.C. is an inventor on patent applications related to this work filed by KTH Holding AB, KTH Royal Institute of Technology (Nos. US 17/255143, EP 3811121, CN 201980043247.2). The authors declare no other competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Data used in generating the main figures are available at Zenodo (DOI: 10.5281/zenodo.4642453).

- Copyright © 2021 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).