Multiple Monoenergetic Gamma Radiography (MMGR) with a compact superconducting cyclotron

Smuggling of special nuclear materials (SNM) and nuclear devices through borders and ports of entry constitutes a major risk to global security. Thus, technologies are needed to reliably screen the flow of commerce for the presence high-$Z$ materials such as uranium and plutonium. Here we present an experimental proof-of-concept of a technique which uses inelastic ($p,p'$) nuclear reactions to generate monoenergetic photons, which provide means to measure the areal density and the effective-$Z$ of an object with an accuracy which surpasses that achieved by current methods. We use an ION-12$^{ \text{SC}}$ superconducting 12 MeV proton cyclotron to produce 4.4, 6.1, 6.9, and 7.1 MeV photons from a variety of nuclear reactions. Using these photons in a transmission mode we show that we are able to accurately reconstruct the areal densities and effective-$Z$ ($Z_{\text{eff}}$) of a mock cargo. This methodology could enable mobile applications to screen commercial cargoes with high material specificity, providing a means of distinguishing benign materials from SNM such as uranium and plutonium.


I. INTRODUCTION
Within the field of nuclear security cargo security focuses on the screening of the flow of commerce for the presence of special nuclear materials (SNM) (such as plutonium and uranium) fully assembled weapons, as well as other non-nuclear contraband. The detection of SNM in particular is important in the context of nuclear terrorism. It is estimated that the economic impact alone in case of the detonation of a crude nuclear weapon in an urban setting can cause damage in excess of $1 trillion [1]. As many as 50,000 maritime intermodal containers enter the United States daily through a variety of maritime ports and other pathways [2]. Given the density of their content [3], their anonymity, and the speed at which these need to be screened, screening cargo containers for nuclear materials presents a major challenge to national and global security. In this paper we present a new radiographic concept that uses transmission of photons of precise and well known energies produced by a superconducting cyclotron to achieve accurate determination of an object's areal density and material type.
Modern passive detection systems in ports of entry have the capability to detect the radiation signals from spontaneous nuclear decay from some types of SNM. However, these signals have a low signal-to-background ratio when the source is shielded, rendering passive detection mostly ineffective against well-informed smugglers. In particular due to the low spontaneous fission rate of 235 U, highly enriched uranium is very difficult to detect * Currently with the MITRE Corporation, Bedford, Massachusetts 01730, USA. The authors affiliation with the MITRE Corporation is provided for identification purposes only and is not intended to convey or imply MITREs concurrence with, or support for, the positions, opinions, or viewpoints expressed in this work. † E-mail: aregjan@mit.edu via these methods. Therefore, scanning of cargoes that potentially contain shielded SNM requires active interrogation and transmission radiography techniques that utilize external radiation sources. Furthermore, while screening for radiation may reveal some SNM, it does not address the auxiliary -and equally important for the transportation industry and customs agencies -goal of detecting other types of contraband, such as smuggling of economic goods. Active interrogation involves directing nuclear radiation such as photons and neutrons to the object and subsequently measuring the secondary radiation, such as delayed and prompt neutrons from photofission to gain information about the content of the object [4][5][6]. Transmission radiography involves measuring the attenuated flux of the interrogating beam's transmitted particles, comparing it with the incident flux, and using that to infer the areal density and possibly the material type subtending a particular pixel. An additional strength of radiography, when compared to other active methods, is its ability to achieve imaging and thus allow detection of conventional contraband -which is the primary goal for many customs agencies. While radiographic techniques are primarily based on photons, radiographic techniques involving other particles do exist such as using neutron beams [7,8]. For a high-level discussion of other active interrogation methods see Runkle et al. [9].
Prior research on MMGR is discussed in Sec. I. It was conducted at the MIT Bates Research and Engineering Center and has been described in detail by Henderson et al. [10]. This research effort however had a significant disadvantage, in the form of a large accelerator and the use of 11 B(d, nγ)C 12 reaction which along with photons produces fast neutrons. The large size limits RFQ's fieldability, while the neutrons contribute to radiation dose. In this study we build upon the prior work by Henderson et al. [10], instead using a more fieldable platform. It uses a compact, relatively light and low-power superconducting cyclotron of a size of 2.5 × 2 × 2 m 3 [11,12].
Using photons from the neutron-less 12 C (p, p γ) 12 C and 16 O(p, p γ) 16 O reactions we show that it is possible to accurately reconstruct the areal density and Z eff of a variety of homogeneous and heterogeneous objects with Z eff values in the range of 13-92.
The attenuation of photons in materials is approximated by the following equation: Here A is the attenuation factor, I is the transmitted photon intensity, I 0 is the source intensity, x is the areal density, and µ is the total mass attenuation coefficient from photoelectric effect, Compton scattering (CS), and pair production (PP). As indicated by Eq. 1 radiography at a single energy depends on both X and µ, the later being a function of both energy and Z, resulting in an underdetermined equation of multiple unknowns. Thus measurements at multiple energies may allow for simulatneous determination of X and µ, the later allowing to infer Z. As more well-separated lines are available the system of equations becomes more overdetermined, allowing for higher accuracy in the inferences. Thus, one method to determine the effective atomic number (Z eff ) value and areal density (x) of the cargo material is to exploit the Z and energy dependence of µ and to take transmission measurements at two or more energies [10,13,14]. Above 4 MeV the major photon interaction mechanisms are CS and e + e − PP, and their contributions to the mass attenuation coefficient can be described as follows [10,15]: where µ cs and µ pp are the mass attenuation coefficient of CS and PP with the corresponding interaction cross sections σ cs and σ pp , N A is the Avogadro number, A is the atomic number of the material, m e is the rest mass of electron, and E is the energy of the photon. As discussed in Henderson et al., the f (E) term in the pair production cross section estimation is a function of energy with negligible dependence on Z except for the Bohr correction.
Since the Z/A ratio for most isotopes is in the range 0.4-0.5, the attenuation from CS does not significantly depend on Z. It primarily depends on the areal density x and energy of the photon. Similarly, the mass attenuation coefficient of pair production is linearly dependent on the Z value of the material. Thus the measurement of total µ can be used to infer Z. CS and PP dominate at different energy ranges. For high-Z materials CS dominates from 1 to 3.5 MeV, hence in that range the the attenuation is primarily dependent on areal density. PP is the dominant interaction above 5 MeV. See Knoll et al., Figure 2.20 for a more detailed description. Thus, a measurement at lower energies can allow for a determination of areal density, with a subsequent measurement at higher energies allowing to infer Z -as shown by Henderson et al., Eq. 2, the measurement of attenuation at two energies with clearly dominating processes could be used to infer Z: ln(A(E PP ))/ ln(A(E CS )) ∝ Z. Since the lowest energy line used in this study is 4.4 MeV, a more sophisticated analysis is required, however.

A. Compact superconducting cyclotron for MMGR
The accelerator used in this study is a 12 MeV isochronous proton cyclotron made by Ionetix [11]. It has an internal target and is designed to produce the 13 N radioisotope for positron emission tomography (PET). The cyclotron operates in continuous wave mode, thus reducing pulse pileup in the detectors. To create the needed photons, the internal target uses a graphite collimator and a static water pocket with a 50 µm-thick aluminum window. The impinging proton beam creates photons through (p, p γ reactions. Since the energy threshold of the 12 C(p, n) 12 N and 16 O(p, pn) 15 O reactions are 19.64 and 16.65 MeV respectively, the only sources of neutrons are the (p, n) and (p, pn) reactions on aluminum [17]. The strongest observed photon line is at 4.44 MeV from the 12 C(p, p γ) 12 C reaction. From oxygen (in water), photons are produced by the de-excitation of 16 O * created through (p, p ) reaction, with energies of 6.13, 6.92, and 7.12 MeV. Using existing cross-section data for the above reactions [18][19][20][21][22][23] and stopping power data tables [24] for proton, the thick target yields can be calculated by performing a numerical integration using the following equa- A where E th is the threshold energy of corresponding (p, p γ) reaction, E p is the energy of the protons (12 MeV), σ(E) is cross section of the the (p, p γ) reaction, dE dx is the energy loss of the proton, and ρ is the density of the material. Figure 1 plots the photon yield from different reactions at different proton energies, as determined by the numerical integration.  Protons are accelerated to approximately 12 MeV, to an orbit with a radius of 14 cm. Before striking the internal water target the beam is collimated by a broad graphite collimator. While for most (p, p γ) reactions the photons are emitted approximately isotropically, in some reactions there is a factor of two difference between some emission angles [18,[20][21][22]25]. Some of the emitted photons then traverse the 5 cm steel wall of the cyclotron and pass through the 10.2 cm × 83.3 cm opening in the concrete wall. In order to minimize neutron activation and background radiation outside of the accelerator room, a block of 20.6 cm thick borated high density polyethylene (HDPE) is placed at the opening 251.8 cm from the target. The photons are further collimated by a 5.1 cm × 20.3 cm lead collimator. The photon beam then transits through the mock cargo and is measured by detector 1 (Det 1) at a distance of 399 cm away from the proton target. Another detector (Det 0) is positioned 45.7 cm to the side of Det 1, which is equivalent to an angle of 6.54 degrees relative to the photon beam direction. Det 0 thus measures the photons directly from the target, and as such serves as a normalization in the analysis. The detectors used to measure the transmitted spectra consist of a 3.81 cm × 3.81 cm cyllindrical LaBr 3 (Ce) scintillator, allowing high energy resolution of ∼3% (FWHM) and fast processing (0.016 µs primary decay time) at the energies of interest [26,27]. The detector pulses are digitized and processed using CAEN V1725 waveform digitizer operating in digital pulse processing pulse shape discrimination (DPP-PSD) mode [28]. The waveform digitizer is controlled using the ADAQ (AIMS Data AcQuisition) framework, producing data files for subsequent analysis [29]. The digitizer recorded the arrival time and energy of each pulse to form a transmitted energy spectrum. For both detectors, the energy trigger threshold are set to approximately 0.8 MeV, filtering the majority of low energy background counts and the 511 keV signals. The integration windows are set at 0.35 µs which allowed for the full capture of the scintillation light pulses from the LaBr 3 detector.

C. Experimental Cargoes (Materials) and dose measured
A variety of homogeneous and heterogeneous (mixed) materials were used as the mock cargo in this experiment . The homogeneous materials included aluminum, copper, tin, lead, and depleted uranium, with a range of areal densities. To determine the Z eff value of the heterogeneous cargo, we first compute the effective µ at 4.44, 6.13, 6.92, and 7.12 MeV for the mixtures using the following formula: where W i is the mass fraction of material i in the mixture and µ i (E) is the mass attenuation coefficient of material i at energy E as found in XCOM database [30]. The calculated µ(E) of the mixture are then compared to the list of homogeneous materials in the NIST database at the corresponding photon energy E. The Z values of the homogeneous material with the best matching mass attenuation coefficient at the four energies are then averaged and used as the actual Z eff for the heterogeneous material. Each of the radiography experiments performed consists of two measurements. The first is a 5 minute run with no material (open air) while the second is a 60 minute run with a mock cargo material. During the analysis the 60 minute run with the mock cargo is divided into 10 equal subsets of data. The reconstruction analysis is performed on individual subsets, allowing for empirical determination of the fluctuations and accuracy in the reconstructions. Over a course of 19 different experiments, the average current measured on the water target and the graphite collimator are 6.0 µA and 0.56 µA, respectively.

A. Spectral Analysis
The transmission data is analyzed in ROOT, an object oriented data analysis library [31]. Each spectrum is calibrated using the following procedures: the 2.21, 3.00, 3.42, 3.93, 4.44, 5.11, 5.62, 6.13, 6.41, 6.61, 6.92, and 7.12 MeV peaks are fitted with a Gaussian function using the TFit class; the results of these fits are then used in a second order polynomial calibration fit; a peak search and background subtraction is then performed on calibrated spectra using the ROOT TSpectrum class. An example of a resulting spectrum can be seen in Figure 3. The TSpectrum background estimation did not include the 7.12 and 6.92 MeV peaks and their first escapes, as no lines are present above those peaks, making the region encompassing them background-free. For a detailed description of the TSpectrum refer to Ref [32].
After calibration and background subtraction, seven peaks are fitted with a Gaussian at 3.42, 3.93, 4.44, 5.11, 5.62, 6.41, and 7.12 MeV. The numbers of detected counts at 3.42, 3.93, 4.44, 5.11, and 5.62 peaks are determined by tallying counts within two standard deviations from the mean of each fitted peak. The 3.42, 3.93, and 4.44 MeV counts, which constitute the photopeak and escape peask of the 4.44 MeV photons, are summed as a single effective 4.44 MeV tally. For the 6.13 MeV counts, only the escape peaks at 5.11 and 5.62 MeV are included: the 6.13 MeV peak consist of both 6.13 MeV photopeak and the 7.12 MeV double escape peak, and there is no simple way of analyzing their contribution individually. An example of the tally regions are highlighted in Figure 3. The region highlighted in purple indicates the summation of total counts in the 6.275 ≤ E ≤ 7.242 MeV region. These counts originated from the 6.92 and 7.12 MeV photons alone.

B. Reconstruction Algorithm
After obtaining the transmitted counts from the spectra in the three regions of interest, the counts measured by Det 1 (on-axis) is first normalized by counts in Det 0 (off-axis) as follows to minimize the various time-dependent systematic effects: C0(E) , where C 0 (E) and C 1 (E) are the recorded counts from Det 0 and Det 1 at energy E. These systematic effects could originate from beam focusing, vacuum quality, and the temperature of internal parts of the cyclotron ultimately leading to variations in beam position and energy, which in its turn leads to changes in photon yields over time. The normalized counts are then used to determine the transmission ratio: where C cargo (E) and C air (E) are the normalized counts at energy E for mock cargo and open beam. This ratio also achieves the cancellation of geometry-dependent systematic effects. As a key step of the reconstruction this ratio is then compared to calculated ratios predicted from a variety of possible values of areal density and Z. For Z of 1 to 100 and areal density x of 1 to 150 g/cm 2 a table of ratios R calc (E, Z, x) is calculated at each combination of Z and x using Eq. 1 and table of µ at Ref. [30]. With the computed table of attenuation and measured transmission rations at different energies, a figure of merit F is defined: where E 1 is 4.44 MeV and E 2 is 6.13 MeV, and σ[] represents the statistical uncertainty of the nominator, as determined by a propagation of errors returned by the TSpectrum class. In the summation term is primarily sensitive to Z eff . The second term is additionally sensitive to the areal density. The process of the reconstruction then searches for a global minimum in F space and identifies the corresponding values of Z and x as the reconstructed values. Figure 4 shows a heat map of calculated F for the case of a copper cargo with areal density of 60.3 g/cm 2 , with reconstructed values of Z eff = 26 and x = 60 g/cm 2 very close to the actual numbers.

IV. RESULTS
As mentioned in Sec. II C, each of the 60 minute experimental transmission datum are equally divided into ten smaller data subsets. The radiography analysis described in Sec. III A is then performed on each of the data sets. The mean and standard deviation of the reconstructed Z and x are computed across the the data set, and these are compared to the actual known values of the mock cargo. These reconstruction results are plotted in Figure 5 and 6 for homogeneous and heterogeneous cargoes, respectively. In most cases the reconstructed values are within the statistical uncertainty of the actual values. The large size of the error bars for the higher Z and higher x materials is a reflection of limited statistics of the transmitted counts, due to increased attenuation at those values. It is now possible to compare the reconstructed and actual values of Z eff and x in a χ 2 test, as a way of determining whether the deviations are due to statistics or some residual, unaccounted-for systematic effects. The reduced χ 2 of Z eff and x between the average reconstructed values and the actual values for all 19 experiments are 0.90 and 1.28, respectively. These correspond to p-values of 0.58 and 0.19, respectively, indicating that the deviations are dominated by statistics. As seen in Figure 5, there is significant separation between low-Z material such as aluminum, medium-Z material such as tin, and high-Z material such as lead. Furthermore, there is a distinct separation between materials with different x. The difference of reconstructed Z between lead and uranium are small, however, partially due to the insufficient photon transmission statistics in the tests with high-Z materials.

V. CONCLUSIONS AND FUTURE WORK
In the experiments presented in this study we demonstrated the feasibility of using monoenergetic photons generated from nuclear reactions to accurately reconstruct the areal density and Z eff of a variety of homogeneous and heterogeneous mock cargoes with a broad range of Z. The precision of the measurements was limited by the statistical precision of the collected data. Future work should focus on using accelerators with higher energies, as this will increase the yield of the photons, as shown in calculations described in Figure 1. For example, increasing the proton beam energy from 12 MeV to 19 MeV would more than double the photon yields, while being only slightly above the neutron production thresholds. While the experiment used beam currents of 5 µA, higher beam currents are necessary to achieve a discrimination between lead and uranium: based on the results in this study we estimate that a 300 µA, 12 MeV proton beam is required for a 5σ discrimination between Pb and U. Finally, advances in high temperature superconductors may in the future allow for cyclotrons with stronger magnetic fields, leading to smaller sizes, lower cost, and thus increased applicability.

ACKNOWLEDGMENTS
This work was supported in part by the U.S. Department of Homeland Security Domestic Nuclear Detection Office under a competitively awarded collaborative research ARI-LA Award, ECCS-1348328. B.S.H. gratefully acknowledges the support of the Stanton Foundation Nuclear Security Fellowship. The research also benefited from generous funding by Andiscern Corporation, which additionally provided the accelerator used in the study. The authors are grateful to Richard C. Lanza, who developed some of the initial ideas behind this work, for his support, encouragement, and valuable advice. The authors wish to thank Steven J. Jepeal and Prof. Zachary Hartwig for the collaborative effort involving the ION-12 SC that made this work possible.