Fissile material detection using neutron time-correlations from photofission

The detection of special nuclear materials (SNM) in commercial cargoes is a major objective in the field of nuclear security. In this work we investigate the use of two-neutron time-correlations from photo-fission using the Prompt Neutrons from Photofission (PNPF) detectors in Passport Systems Inc.'s (PSI) Shielded Nuclear Alarm Resolution (SNAR) platform~\cite{pnpf} for the purpose of detecting $\sim$5~kg quantities of fissionable materials in seconds. The goal of this effort was to extend the secondary scan mode of this system to differentiate fissile materials, such as highly enriched uranium, from fissionable materials, such as low enriched and depleted uranium (LEU and DU). Experiments were performed using a variety of material samples, and data were analyzed using the variance-over-mean technique referred to as $Y_{2F}$ or Feynman-$\alpha$. Results were compared to computational models to improve our ability to predict system performance for distinguishing fissile materials. Simulations were then combined with empirical formulas to generate receiver operating characteristics (ROC) curves for a variety of shielding scenarios. We show that a 10 second screening with a 200~$\mu$A 9~MeV X-ray beam is sufficient to differentiate kilogram quantities of HEU from DU in various shielding scenarios in a standard cargo container.

The detection of special nuclear materials (SNM) in commercial cargoes is a major objective in the field of nuclear security. In this work we investigate the use of two-neutron time-correlations from photo-fission using the Prompt Neutrons from Photofission (PNPF) detectors in Passport Systems Inc.'s (PSI) Shielded Nuclear Alarm Resolution (SNAR) platform for the purpose of detecting ∼5 kg quantities of fissionable materials in seconds. The goal of this effort was to extend the secondary scan mode of this system to differentiate fissile materials, such as highly enriched uranium, from fissionable materials, such as low enriched and depleted uranium (LEU and DU). Experiments were performed using a variety of material samples, and data were analyzed using the variance-over-mean technique referred to as Y2F or Feynman-α. Results were compared to computational models to improve our ability to predict system performance for distinguishing fissile materials. Simulations were then combined with empirical formulas to generate receiver operating characteristics (ROC) curves for a variety of shielding scenarios. We show that a 10 second screening with a 200 µA 9 MeV X-ray beam is sufficient to differentiate kilogram quantities of HEU from DU in various shielding scenarios in a standard cargo container.

I. INTRODUCTION
Fissile materials refer to materials which, due to their nuclear structure, allow for sustained fission chains. The two most common isotopes which form the basis of fissile materials are 235 U and 239 Pu. Methods for detecting the presence of fissile materials support the goals of national and international programs for nuclear-nonproliferation. An extensive literature exists on potential identification schemes using active interrogation including neutron probes and delayed and prompt neutron signals [1][2][3][4][5], photon probes and delayed neutron signals [6][7][8][9][10][11][12], photon and neutron probes and delayed neutron signals [13][14][15], and photon probes and fission product radiation [16,17]. A general review of the various concepts can be found in Ref. [18]. All of these methods use single-particle signals. Multi-particle schemes, usually two-and three-neutron signals have also been pursued, with photon [19,20] and neutron [21][22][23][24][25] beams as well as in passive interrogation [26][27][28][29][30][31][32], where ambient radiations are used to induce the signals. Here we report on an active interrogation method using 9 MeV Bremsstrahlung photons to induce the emission of time-correlated, prompt neutrons measured within the Prompt Neutrons from Photofission (PNPF) system developed at Passport Systems Inc. (PSI) [33][34][35].

II. TIME-CORRELATION OF FISSION-CHAIN NEUTRONS
A fission chain reaction occurs in fissile material when the neutrons from each fission diffuse through the material and induce subsequent fissions in other nuclei. The neutrons emitted from this process are highly correlated, resulting in neutron count distributions that deviate significantly from Poisson distributions produced by random neutron events.
The theory of fission-chain correlations was initially developed by Richard Feynman while at Los Alamos [36]. The aim of that research was to describe the neutron fluctuations in a reactor pile where the measured neutrons originate from fission chains and random decays. To measure deviations from the Poisson expectation, Feynman defined a normalized second moment, where n is the measured neutron count per unit time, and n and n 2 represent averages over a series of time-gates. Subtracting the mean from the variance enforces Y 2F = 0 for the special case of a Poisson distribution, and division by the mean ensures that the quantity is independent of rate. The subscript "2F", not used in the original paper, has been added by others to credit Feynman and to denote that this quantity is derived from the second moment. This statistic has been used extensively in the study of nuclear reactor cores [37][38][39][40][41], and it has been extended for time varying sources of the fluctuations, e.g. in accelerator driven cores [42][43][44][45][46][47][48]. In these systems the number of neutrons generated spontaneously from radioactive decay greatly exceeds those generated by the introduced beam.
For neutron-induced fission-chains from an idealized point source, a time-dependent expression for Y 2F has been developed by Snyderman and Prasad [54][55][56], where, λ = inverse of the neutron correlation time, a convolution of fission-chain and neutron transit times, T = time-gate during which neutrons are counted, = neutron detection efficiency, M = neutron multiplication, 1 1−k , where k = pν is k-effective and p is the neutron-induced fission probability, ν = the first moment, equal to the mean number for the neutron distribution from a single induced fission, ν 2 = the second combinatorial moment, half the variance for the neutron distribution from a single fission.
Eq. 2 illustrates utility of Y 2F , whcih increases with the square of the multiplication and is therefore sensitive to the enrichment if the geometry of the object and efficiency of the detector are known. One can construct higher order moments (i.e. Y 3F,4F ) that depend on higher orders of the multiplication, but a corresponding dependence on higher orders of the efficiency makes these statistics more suitable for large acceptance detectors and/or long integration times due to the reduced coincidence rates. It is also important to note that Eq. 2 is appropriate for neutron beams. In Sec. V we will modify this expression for photon-induced fission-chains in extended objects.

III. EXPERIMENTAL SETUP
To study the feasibility of using photo-induced fissions chains to identify and characterize fissile materials we use a subset of the components of the PSI SNAR platform. The photon beam is generated from electron bremsstrahlung on a water-cooled radiator. The electron current is adjustable and uniform with high duty-cycle. For all experiments a 9 MeV electron kinetic energy was used. The beam currents are adjustable from 100 to 500 µAmps. For the results presented here, the beam current was set to 200 µAmps to avoid the non-Poisson effects of the data acquisition (DAQ) dead time losses. This corresponds to a rate of 2 × 10 12 Hz for photons with energies between 2 and 9 MeV, incident on an approximately 10 × 10 cm square spot size at the target location. The photon beam energy spectrum for the 9 MeV Bremsstrahlung radiator is shown in Fig. 1. The full SNAR configuration is shown in Fig. 2. The PNPF detectors consist of EJ-309 5 inch diameter liquid scintillator detectors coupled to 5" Hamamatsu photomutilpier tubes arranged in two sets of 2 × 8 arrays placed on opposite sides of the cargo container. The PNPF arrays have ∼5 cm thick high density polyethylene (HDPE) inserts placed in between the detectors to reduce adjacent detector cross-talk that can lead to an artificial Y 2F signal. Before and after tests show that these inserts reduce the cross-talk component by a factor of four.
Measurements of Y 2F were made for a set of objects designed to study the relationship between Y 2F and neutron multiplication. These objects, listed in Table 1, were constructed from discs of depleted uranium (DU) and highly enriched uranium (HEU), and blocks of low enriched uranium (LEU). The objects were arranged in geometries to achieve a higher range of multiplication values. Fig. 3 shows the stacking arrangement for the three HEU discs. The highest multiplication values were achieved with an interleaved stack of LEU and HEU with HDPE moderators placed above and below the stack to provide neutron reflection. Multiplication values were calculated using MCNP 6.2 [57] run in the kcode mode. The control benchmark was obtained with a beryllium block, in order to generate a large quantity of uncorrelated photo-neutrons anticipated to have no measurable Y 2F value. Data were collected from a series of 600-second exposures, which were then combined for each object during the data analysis. Approximately ten exposures were collected for each uranium configuration, with additional exposures taken for the beryllium object. List of objects exposed to the 9 MeV Bremsstrahlung beam within the Passport Systems SNAR facility.

A. Neutron Identification
The main component of the detection system consists of an array of EJ-309 scintillators, coupled to Hamamatsu photomultiplier tubes (PMT). The choice of the PMTs was based on their < 3 ns rise time, thus retaining the sensitivity to the difference between the two main fast components of the light production from the triplet-triplet annihilation. The increase in triplet-triplet annihilation from the much larger ionization density of the neutronrecoiled proton tracks results in a distinct delay in light output. Neutrons can be distinguished from photons by their relative fractions of fast and slow light output using various Pulse Shape Descrimination (PSD) techniques, and their deposited energies can be quantified. The left panel of Fig. 4 shows a histogram of PSD vs. deposited energy, in electron equivalents, for a set of detector events. The distribution is then split into individual energy bins, and the PSD distribution for every energy bin is fitted with two Gaussians, as can be seen in the right panel of Fig. 4, which shows the 2D cut used to separate neutrons from photons. However, this is not sufficient to suppress the very large photon population, created by the intense 200 µA photon beam. Therefore, more sophisticated pulse shape discrimination algorithms were developed to reject the photon pileups, which otherwise would contribute a significant background to the neutron population [35]. Overall, the combined pulse shape analysis allows a suppression of the photon population by a factor of ×10 7 , limiting the beam-on background count to only the cosmogenic neutrons.

B. Multiplicity analysis of neutron events
The experimental data were collected in a series of 600-second exposures, with some exposures ending prematurely due to disruptions within the data acquisition system due to the extended run times and need for data synchronization. The identified neutrons in each run were arranged in time-ordered sequence to facilitate binning in time-gates of varying duration. For clarity, we rewrite Eq. 1 above with explicit averages, where sums are taken over sequential time-gates, n refers to the bin-number corresponding to the neutrons detected within a specific time-gate, and b n is the normalized probability for detecting n neutrons. The errors are calculated using the full covariance matrix for a multinomial distribution [58], The Y 2F values and errors for each exposure were calculated separately and checked for outliers before combining the exposures to calculate a final value and error for each object. The Y 2F measurements are listed in Table 2 and plotted in Fig. 5 along with a set of corrected Y 2F values that will be discussed in the next section. We note that the uncorrected values illustrate the time dependence expected from Eq. 2. One can also see from the uncorrected beryllium measurements that a significant non-fission component contributes to the uncorrected Y 2F values.

C. Cross-talk Cuts and Corrections
Neutrons that deposit energy in more than one detector will produce a correlated signal that contributes to the Y 2F term. This effect is shown clearly in Fig. 6. The left panel shows the distribution of pairs of neutron hits for one DU run as a function of the coincidence time and channel separation for channels in the same array column. Adjacent channels within a single column are separated by a distance of 18.5 cm. For ∆chan = 1, 90% of these channels are separated by this distance, with the remaining 10% separated by a distance of 2 meters or greater. ∆chan = 2, 80% are separated by a distance of 37 cm. The corresponding distribution for a GEANT4 [59] simulation is shown on the the right. The DAQ trigger logic rejects new signals within 256 ns after a given trigger, leaving zero entries ∆chan = 0 row of Fig. 6 creating the empty bins in the top row that is obscured from view. A significant cross-talk effect is evident in the vertically adjacent and next-nearest neighbor channels (∆chan = 1, 2). These enhancements are also visible in the GEANT4 simulation shown in the right panel. Additional time-dependent structure in the data is attributed to the neutron-correlations, which were not included in this simulation.
We have corrected for this cross-talk effect on adjacent channels by fitting the time-dependence to a linear function at larger times and extrapolating into the cross-talk region, as shown in the left panel of Fig. 7. The fit was performed for time differences greater than 100 ns, and the integral of the fit was used to calculate the expected signal in the absence of cross-talk. A similar analysis was performed for next-nearest channels. The ratio of actual counts to extrapolated counts for adjacent and next-nearest channels is shown in the right panel Fig. 7. Note that for ∆ch = 2 the correction returns to unity for times that are smaller than the transit time between the next-nearest neighbor channels. Note that this correction was only used for time differences less than 100 ns, and was not used in the regions above where deviations from unity are dominated by statistical fluctuations. A cross-talk correction was also calculated for the horizontal neighbors separated by 2 meters (∆ch = 8) and was found to be negligible.
The Y 2F values corrected for cross-talk are listed on the second row of each time-gate in Table 2  as solid symbols in Fig. 5. Note that the corrected Y 2F values for beryllium are consistent with zero as expected. The difference between Y 2F values for the DU and HEU stacked discs, already visible before the application of the correction, is still easily discernable in Fig. 5. Although the Y 2F difference grows larger with increasing time-gate, the most significant difference for this measurement, approaching 5-sigma, occurs in the region of 50-100 ns. A time-gate of 100 ns is applied to the analysis of all subsequent measurements.

V. RESULTS AND SIMULATION BENCHMARKS
The Y 2F values for 100 ns time-gates and cross-talk corrections for the complete set of objects listed in Table 1 are plotted in Fig. 8  increase for the composite systems. These results are compared to MCNP simulations of the PSI beam incident on the fissile objects. The detector acceptance was approximated by applying a ±15-degree angular cut above and below the horizontal. The angular cut matches the angle subtended by the PNPF detector array, and it is required to account for the polar-angle dependence of the neutron emission rates, which varies with object and is strongest for the LEU+HEU+HDPE composite objects. A minimum energy cut of 1.5 MeV is also applied to the neutrons to approximate the online energy cut implemented in the DAQ. Neutrons within this acceptance region were used to calculated the Y 2F in the same manner as the data, using 100 ns time-gates, but without the need for cross-talk corrections. A geometric efficiency factor of 0.25% was applied to the simulated Y 2F values. This value was determined empirically and is comparable to the 0.243% geometric detector efficiency calculated in GEANT4. The MCNP results reproduce the overall trends of the data, although the simulations under-predict the decrease for DU and over-predict the increase for the LEU+HEU composite object. The highest multiplication objects consisting of the interleaved LEU and HEU objects with HDPE moderation are reasonably well-reproduced by the MCNP simulations. This complex dependence of the Y 2F values that was observed in the data and qualitatively reproduced by the simulation led us to reconsider some of the assumptions within Eq. 2. One important modification is to account for the fact that in photon-induced fissions, the initial fission differs from those of subsequent neutron-induced fissions in the chain. This is analogous to the modification required for spontaneous fission-chains developed by Prasad and Snyderman [55]. The initial photo-fission is especially important for small multiplication values, where the neutrons from the initial fission form a non-negligible contribution to Y 2F . We introduce ν γ for the mean number of neutrons produced in a photo-fission, and ν γ2 for the half-variance. We also account for the non-fission photo-production and absorption in the extended object by adding the terms f o for the fraction of neutrons emitted from fission over the total number of neutrons produced by all nuclear reactions and o for the neutron absorption within the object. The neutron detection efficiency is denoted separately as d . Accounting for these additional factors yields the following expression for Y 2F for photon-induced fissions in an extended object, The mean number of neutrons per photon-fission, ν g has been measured by Caldwell et al. [60] and evaluated by Chadwick [61], however, the value for ν γ2 has not yet been measured. To estimate this value within the simulation, a  TABLE 3. First and second combinatorial moments for photon and neutron induced fission for various levels of uranium enrichment. Neutron induced moments are calculated from Holden and Zucker tables for 1-2 MeV neutrons reported in the documentation for the MCNP fission package. Photon induced moments are calculated from the photo-fission model using the PSI beam energy spectrum.
separate simulation of the MCNP photo-fission package was performed using the PSI energy spectrum as input. The mean and half-variance values for photon-induced fissions obtained in this way are listed in Table 3. We also list the corresponding values for neutron-induced fissions as reported by Zucker and Holden [62]. To compare our simulations to Eq. 5, which does not include energy dependence, we remove all energy cuts from the analysis, and we lengthen the Y 2F time-gate for this simulation analysis to 1 ms, to encompass the correlation-time for all neutrons, including thermal. The values for f o and o were taken from the simulations and a perfect detector efficiency ( d = 1) was assumed. The full comparison between this MCNP simulation and Eq. 5 is shown in Fig. 9. Eq. 5 accounts for the full multiplication dependence observed in the simulations, but consistently under-predicts the Y 2F values by 20%. One possible explanation for this discrepancy is the lack of energy dependent terms in Eq. 5. It is also possible that our estimate for the unmeasured ν γ2 differs from the value used in the full simulation. The Y 2F values can be very sensitive to the value of ν γ2 for objects with low multiplication. To illustrate the sensitivity of Y 2F to the value of ν γ2 , we vary its value by ±10%, as shown by the upper and lower bands.

VI. SYSTEM PERFORMANCE PREDICTIONS
With the MCNP simulations benchmarked by experimental results and qualitatively supported by an underlying theory, additional simulations of shorter exposures with larger masses were performed to predict PNPF system performance for distinguishing between DU and HEU in shielded configurations. For this study we performed simulations for a set of 24 solid spheres of DU and HEU spanning a range of masses from 10 grams up to 50 kilograms. The masses, radii, and multiplication values are given in Table 4. MCNP simulations were performed for these objects using a PSI beam of 200 µA beam for a 10 second exposure. The emitted neutrons were collected over the entire acceptance and multiplied by the overall geometric efficiency factor of 0.25% to match the experiment. A lower minimum energy cut of E>1MeV was used, based upon potential improvements to the neutron identification, and a nominal time-gate of 100 ns was used for calculated Y 2F values.
Although the Y 2F values are calculated directly from the simulations, we chose to use experimental data to extrapolate errors. Fig. 10 shows a uniform, linear dependence in the Y 2F error as a function of N −1/2 counts , independent of object type. The data are well fit to the following functional form, We note that this formula does not extrapolate to zero error for zero counts, consistent with the fact that the errors are not strictly Poisson-distributed. This scaling is used to rescale our simulations for different objects and shielding configurations to predict the dependence on exposure time. Using the simulated Y 2F values and errors extrapolated from data, we use a Gaussian classifier to establish Receiver Operating Characteristics (ROC) curves for the purpose of distinguishing between HEU and DU objects. Using the HEU as the alarm scenario and the DU as the null scenario we can assign as "positive" value based on a Gaussian probability threshold: where Y exp is the measured Y 2F and σ is the extrapolated error. The ROC curves for different unshielded object masses are shown in Fig. 11 for spherical objects of 2, 5, 10, and 20 kg. For every object the ROC values were determined for the current PSI detector configuration (solid lines) and for one with 4× increased geometric efficiency. The analysis shows that while detection of small objects, e.g. 2kg, may be difficult, for larger objects the signal becomes large enough for a fast detection. This becomes especially true for a larger detector array. Note that the positive detection of 20 kg sphere with the 4× augmented detector array approaches 100% with a very low rate for false positives.

Scaling Y 2F
The scaling of the uncertainty Y 2F is described by Eq. 6. This is shown for the data in Fig. 5. The fractional uncertainty can also be plotted to be used in Figure 5: Data for the 300s and full duration runs. Red diamonds are the HEU points, blue diamonds the DU points. The uncertainty scales with 1/ p N counts independent of the e ciency or Y 2F value as expected from Eq. 6.
scaling the Monte Carlo to data. Fig. 6 has the fractional uncertainty plotted for 1/ p N counts and the linear fits for the DU and HEU samples, separately. The slopes of the fits to the data points is related to the e ciencies of the In addition to characterizing performance for unshielded objects, the system was also studied for shielded configurations, where the 5 kg objects were surrounded by iron and HDPE blocks of varying thickness. The first shielding scenario has the primary effect of reducing the incident photon flux, thus reducing statistics in a fixed measurement time. The second scenario's impact primarily consists of moderating the fast neutrons, and thus reduces neutron statistics. Figures 12 shows the ROC curves for the 5 kg object for several thickness of full density iron (left) and HDPE shielding (right).

VII. CONCLUSIONS
The objective of the experimental effort described in this work was to extend the already existing Prompt Neutrons from Photofission (PNPF) technique, described in Ref. [35], to distinguish between fissionable (e.g. DU) and fissile (e.g. HEU, WGPu) materials by using neutron multiplicity analysis. The goal of the initial stages of the experimental program was to develop the necessary statistical and mathematical formalism, the Monte Carlo simulation infrastructure, and to experimentally demonstrate the feasibility of the technique. The results showed a 5-sigma difference in the signal for the two types of objects, proving the feasibility of this methodology.
The next stages of the program focused on additional experiments, involving more complex target configurations using DU, LEU, and HEU. The Y 2F measurements were compared to simulations, which in turn were used to develop a theoretical-numerical model providing a basis for extrapolations to larger objects and shorter exposure times. The results showed that while the Y 2F signal is marginal for practical differentiation of the small objects used in the experiment, for 5kg spherical sizes the neutron multiplication M is large enough to differentiate fissile material types in shielded cargo configurations with reasonable scan times.
While both the feasibility and practicality of the methodology was demonstrated, significant additional work remains to be done. Measurements on a wider range of object sizes, both large and small, would be valuable for further understanding the dynamics of Y 2F . Also high statistics measurement of photo-fission neutron distributions with thin-foil U-235 and U-238 targets with large acceptance detector that would enable precise and accurate determination of ν γ2 would significantly reduce uncertainties for future system performance studies.