2.1 NRF principle

The NRF cross section for absorption via the resonant energy level E_r can be expressed by the Breit–Wigner distribution [41]:

where Γ is the width of the level at E_r, Γ₀ is the partial width for transitions between E_r and the ground state, ħ is the Planck constant, and c is speed of light.

In practice, the NRF cross section should be calculated taking into account Doppler broadening. If the true Voigt profile is approximated as a Gaussian profile, Eq. (1) then becomes [42]

here Δ = $E_{\text{r}}\sqrt{k_{\text{B}}T/M{c}^2}$ is the Doppler width, k_B is the Boltzmann constant, T is the absolute temperature, and M is the mass number of the nucleus. According to Eq. 2, one can calculate the integrated NRF cross section σ_int for ²³⁵U and ²³⁸U. As shown in Table 1, the σ_int values are consistent with experimental data for ²³⁵U and ²³⁸U [41, 42]. The NRF lines at 1734 keV (²³⁵U) and 1782 keV (²³⁸U) have NRF cross sections of 22.0 and 21.9 eV b, respectively. Considering their relatively large σ_int values, these two separate NRF lines are priorities for the identification of uranium, which is selected as a typical SNM here.

Because of conservation of energy and momentum, a free nucleus undergoing NRF will recoil with kinetic energy E_rec, which is determined by the Compton-like formula

where θ is the scattering angle of the photon relative to its incident direction.

NRF is generally considered to occur only between states that differ by two or fewer units of angular momentum. The angular distribution of NRF γ rays is analogous to that of γ-ray cascades. For an NRF interaction of transitions Ja(L₁)Jb(L₂)Jc, where L₁ and L₂ are the multipole orders of excitation and de-excitation, respectively, the angular distribution W(θ) can be written as [43]

where P_2n(cosθ) is the Legendre polynomial expansion, and A_2n is given by

where F_2n(L₁JaJb) and F_2n(L₂JcJb) are constants that depend on the spin states of the transitions and photon multipolarities [44]. For the resonant state at 1782 keV (²³⁸U), the NRF follows a transition sequence of 0 → 1 → 0, whose angular correlation can be expressed as W(θ)=0.75×(1+cos²θ). By contrast, at 1734 keV (²³⁵U), W(θ) depends on the spin, J = 9/2 or 11/2. Because this state can de-excite to the first excited state and ground state of ²³⁵U, several multipolarity combinations are obtained according to the spin selection rule (see Fig. 2). However, it is still impossible to obtain an exact expression of W(θ) because their mixing ratios remain unknown. For simplicity, we employ an isotropic W(θ) for NRF γ-ray emissions in the simulations. In fact, a non-isotropic angular distribution would contribute at most a ~10% fluctuation to the NRF yields in our configuration (see Fig. 1). More details are given in Sec. 4.

Figure 2:Full-screen View

(Color online) W(θ) values for the NRF transition from the 1734 keV resonant state, which has hypothetical spin 9/2 (red) or 11/2 (blue), to the ground state with spin 7/2 (a) and to the first excited state with spin 9/2 (b). In panel (a), three possible multipolarity combinations are allowed: only dipolar (solid lines), only quadrupolar (dotted lines), and one dipolar transition and one quadrupolar transition (dashed lines). In panel (b), four hypothetical multipolarity combinations are allowed: only dipolar (solid lines), only quadrupolar (dotted lines), dipolar-quadrupolar (dot-dashed lines), and quadrupolar-dipolar (dashed lines).

2.2 Scattering NRF spectroscopy

To realize SNM identification and isotope ratio prediction, 1D sNRF spectroscopy is applied. As shown in Fig. 1, a quasi-monochromatic γ-ray beam impinges on the target to be interrogated, causing resonant (NRF) and non-resonant (Compton scattering, pair production, and photoelectric absorption) interactions. The backscattered NRF γ rays are measured by four high-purity germanium (HPGe) detectors (scattering detectors) located at 135° from the beam direction in order to take advantage of the decreasing intensity of non-resonantly backscattered radiation. The horizontal position (x) is varied from -15 to 15 mm in eight steps of 3.75 mm each. A total of eight sNRF γ-ray spectra are obtained.

In sNRF spectroscopy, one can use the NRF cross section σ_NRF(E) and angular distribution W(θ) to construct a semi-analytical expression for the expected NRF counts. For a photon beam of incident flux I(E) interacting with the target, a small part of the photon flux near the resonant energy E_r will undergo resonant (NRF) and non-resonant (atomic) interactions. The resulting NRF yield then produces a double-differential rate of NRF detections in the infinitesimal solid angle dΩ,

where E and E’ are the energy of the incident photons and scattered NRF photons, respectively; L is the thickness of the irradiated target; ϵ(E’) is the intrinsic photopeak detection efficiency; μ_NRF(E) = Nσ_NRF(E) denotes the linear attenuation coefficient, with N being the number density of interrogated isotopes; and μ_eff(E,E’) is the effective attenuation coefficient, which is given by

Here μ_nr(E) and μ_nr(E’) are the non-resonant attenuation coefficients of the incident photons and NRF photons, respectively.

2.3 Transmission NRF spectroscopy

After sNRF spectroscopy is performed, a tNRF-CT technique is applied to perform tomographic imaging (see Fig. 1). The flux of the γ-ray beam transmitted through the target is preferentially attenuated (notched [45]) around the resonant energy E_r because the NRF cross section is much larger than those of the non-resonant interactions. This notched γ-ray beam further impinges on a witness target composed of suspicious isotopes so that the remaining γ rays may undergo NRF in the witness target. Another array of four HPGe detectors (transmission detectors) are located at 135° to record the NRF photons produced at this stage. The resonant attenuation inside the interrogated object is then evaluated. The γ rays transmitted through the witness target are diagnosed by a LaBr₃(Ce) detector to evaluate the non-resonant attenuation. To obtain the CT images, the interrogated object is translated horizontally (x) from -15 to 15 mm (with a step length of 3.75 mm) and rotated (by θ_r) from 0 to 180 (with a step length of 22.5); consequently, a total of 64 sets of spectra are obtained. In addition, a set of spectra without the interrogated object is obtained.

The attenuation factor of on-resonance γ rays at (x, θ_r) can be expressed as

where (μ/ρ)_ave is the average mass attenuation coefficient of the CT target (i.e., the interrogated target) on the incident beam path, and L is the diameter of the CT target. σ_NRF is the NRF reaction cross section of the isotope of interest, and N_t(x, θ_r) is the isotope number density on the γ-ray incident path. For the off-resonance γ rays, σ_NRF is negligible; thus, the attenuation factor of the off-resonance γ rays is

The NRF resonant attenuation factor at (x,θ_r) can be derived as follows [37]:

where -ln(ε_ON)) and -ln(ε_OFF) are the attenuation factors of the on-resonance and off-resonance γ rays, respectively. C_ON(x, θ_r) and C_ON,blank denote the NRF yields recorded by the transmission detectors with and without the CT target, respectively. C_OFF(x, θ_r) and C_OFF,blank are the integration yields of the spectral region of interest (ROI) recorded by the LaBr₃(Ce) detector with and without the CT target, respectively. Note that ε_OFF is an approximate estimate of the atomic attenuation effect of on-resonance γ rays in Eq. 10 when a narrow-band γ-ray beam is used. Consequently, the NRF attenuation factor depends only on N_t(x, θ_r), which is required to reconstruct the CT images of SNMs.

2.4 Simulation algorithm

To model the NRF process in this study, we developed a new class, G4NRF, in the Geant4 toolkit [46, 47]. The pure virtual method G4VUserPhysicslist::ConstructProcess() was implemented in the simulation, and the method AddDiscreteProcess() was used to register the NRF process. Introducing a customized NRF process into the simulation requires the implementation of two features. First, the cross sections for the interaction must be provided; second, the final state resulting from the interaction must be determined. A series of NRF cross sections was calculated using Eq. 2. Information on the final states was obtained using Eq. 4. The transitions to the ground states and first excited states of ^235,238U are considered. The HPGe detectors have an energy resolution of 0.1% (RMS), which can be achieved using present detector technology. The Ge crystals are 10 cm in diameter and 10 cm in length. The full-energy peak efficiency of each HPGe detector was also simulated with the Geant4 toolkit.

3.1 Isotope identification by sNRF signature

In sNRF spectroscopy, the target is irradiated by a photon beam with a Gaussian energy distribution [centroid energy of 1.76 MeV and energy spread of 3% in standard deviation (SD)] and a photon intensity of 10¹⁰ photons per second, which can be readily delivered by a state-of-the-art LCS γ-ray source [48, 49]. Among the eight energy spectra obtained in the 1D scan, the sNRF signatures at 1687 keV (²³⁵U), 1734 keV (²³⁵U), 1737 keV (²³⁸U), and 1782 keV (²³⁸U) appear only in the spectra obtained at x = 1.8 and -1.8 mm. The presence of these sNRF signatures gives a preliminary estimate of the SNM isotopic composition of the interrogated target. This result can potentially reveal a 1D map of SNM isotopes, as reported in the literature [50, 51]. Moreover, these sNRF signatures can potentially be used for the tomographic imaging of multiple isotopes, which is an interesting issue to study.

Figure 3 shows typical energy spectra of γ rays recorded by the scattering detectors at scan points of x = 1.8 and 5.6 mm. The NRF signals are simultaneously observed in the former and disappear in the latter. The NRF γ-ray peak at 1687 keV is caused by the transition from the resonant state of ²³⁵U at 1734 keV to the 9/2^- excited state at 46 keV. The peak near 1734 keV is caused by the transition from the 1734 keV level in ²³⁵U to the ground state and the transition from the 1782 keV state in ²³⁸U to the first excited state at 45 keV (with photon emission at 1737 keV). Note that these two closely spaced NRF lines cannot be well discriminated owing to spectral broadening resulting from imperfect detector resolution. The NRF peaks at 1687, 1734 (or 1737), and 1782 keV are then fitted with four Gaussian distributions on top of an exponentially decaying continuum background. The fitting function for these NRF peaks is written as

Figure 3:Full-screen View

(Color online) (a) Simulated sNRF γ-ray spectra at scan points (x, θ_r) = (1.8 mm, 0) (red line) and (5.6 mm, 0) (black line) recorded by scattering detectors located at 135. There are clear photopeaks at 1687, 1734, and 1782 keV. The density of the uranium rod is set to 19 g/cm³. The isotopic composition of the uranium rod is set to 80% ²³⁵U and 20% ²³⁸U. Single and double escape peaks are not observed because the low-probability escapes are obscured by low-energy photons deposited in the detectors with large volume. (b) Magnified spectrum showing the NRF signal region of 1650–1800 keV. Solid lines represent the fitting curves of four Gaussian peaks plus an exponential background fit.

where c₁ and c₂ describe the shape of the background, and ak, Ek, and σk are the area, mean, and SD fit parameters of the k^th peak. The fitting curve is shown in Fig. 3 (b). The corresponding NRF yields, Y₁₆₈₇, Y_1734,1737, and Y₁₇₈₂, were obtained. Because the branching ratio (denoted as b₁) of the 1734 keV transition to the 1687 keV transition is 100:60(20), the NRF yield for the 1734 keV transition was further deduced as $Y_{1734}^{\text{ded}}=b_1{Y}_{1687}$. Similarly, the yield for the 1737 keV transition, $Y_{1737}^{ded}$, was also deduced. Moreover, the NRF yield for the overlapping peak near 1734 keV was estimated as $Y_{\mathrm{1734,1737}}^{\text{ded}}$ =$Y_{1734}^{\text{ded}}+Y_{1737}^{\text{ded}}$, as shown in Table 2. The extracted $Y_{\mathrm{1734,1737}}^{\text{ded}}$ values agree well with Y_1734,1737, indicating that the effect of branching ratios was implemented appropriately in our simulations. The significance of the sNRF signals of ²³⁵U can be expressed in units of σ as Δ_sNRF = S/δS, where S is the NRF peak yield, and δS is the corresponding statistical error. The result shows that a peak significance of ${\Delta }_{\text{sNRF}}^0$ = 48.9σ can be obtained within an sNRF scan time of $t_{\text{sNRF}}^0$ = 8 s.

The isotope ratio of ²³⁵U to ²³⁸U, ²³⁵U/²³⁸U, is related to the NRF yield ratio, ²³⁵Y/²³⁸Y, by the following equation:

where ²³⁵Y/²³⁸Y is the peak yield ratio, and I(E_238U)/I(E_235U) is the ratio of the incident γ-ray intensity. ^ε238U/^ε235U is the detection efficiency ratio of the HPGe detector calculated using the MCLCSS code [52] in Geant4. br_238U/br_235U is the ratio of the absolute branching ratios. ${\displaystyle \int {\sigma }_{{238}_{\text{U}}}}(E)\text{d}E/{\displaystyle \int {\sigma }_{{235}_{\text{U}}}}(E)\text{d}E$ is given by Eq. 2. W_238U(θ)/W_235U(θ) is the angular momentum ratio of ²³⁸U and ²³⁵U NRF emission [43]. According to Eq. 12, the abundance ratios of ²³⁵U/²³⁸U at x = 1.8 and -1.8 mm for three ²³⁵U enrichments were calculated on the basis of the NRF yields at 1687 and 1782 keV, as shown in Table 3. The predicted isotope ratios are consistent with the theoretical values within the uncertainty. Note, however, that the attenuation of γ rays with different energies as they penetrate the wrapping materials should be considered in Eq. 12 to improve the effectiveness of the isotope ratio prediction for thicker shielding.

Figure 4 shows the simulated NRF yields of the 1782 keV line when the HPGe detectors are located at angles of 90, 115, and 135. It is shown that the NRF yields increase with detection angle because the W(θ) value for the 1782 keV line at 135 is larger than those at the other two angles. Consequently, a detection angle of 135 is employed in our study (see Fig. 1).

Figure 4:Full-screen View

(Color online) The 1782-keV NRF yield as a function of target areal density at detection angles of 90, 115, and 135. Only the statistical uncertainty is considered here. The isotopic composition of the uranium rod is 80% ²³⁵U and 20% ²³⁸U.

3.2 Tomographic imaging using tNRF signature

In tNRF imaging, the uranium rod is ideally assumed to be composed of pure ²³⁵U, and its default density is set to 19 g/cm³ to reduce the computational requirements. The interrogating γ-ray beam has a Gaussian distribution (a centroid energy of 1734 keV and an energy spread of 1% in SD) and a photon intensity of 10¹⁰ photons per second. Figure 5 shows typical spectra obtained by the transmission detectors and LaBr₃(Ce) detectors at measurement points of (1.8 mm, 0), (1.8 mm, 22.5), and (1.8 mm, 45), where the uranium rod, lead rod, and air rod, respectively, are in the path of the interrogating γ-ray beam. As shown in Fig. 5 (a), the energy spectra recorded by the LaBr₃(Ce) detector are different because they depend on the atomic attenuation coefficients of the penetrated materials. The 1734 keV peak intensity recorded by the HPGe detectors at the scan point (1.8 mm, 0) is significantly lower than those at the other two scan points. The reason is that the intensity of the γ-ray beam transmitted through the CT target decreases with θ_r owing to strong resonant absorption.

Figure 5:Full-screen View

(Color online)(a) γ-ray spectra detected by the LaBr₃(Ce) detector at the points (1.8 mm, 0) (blue solid line), (1.8 mm, 22.5) (magenta dashed line), and (1.8 mm, 45) (red dotted line). The energy spectrum of the incident γ-ray beam (black solid line) is also shown. The shaded area represents the ROI at 1694–1774 keV. (b) Magnified spectra showing the tNRF signal of ²³⁵U around 1734 keV recorded by the transmission detectors when the interrogated object was located at (1.8 mm, 0) (blue solid line), (1.8 mm, 22.5) (magenta dashed line), and (1.8 mm, 45) (red dotted line).

As mentioned above, 65 sets of spectra were obtained. On the basis of the NRF peak yields at 1687 and 1734 keV, the ϵ_ON(x, θr) values for these scan points are obtained. Similarly, the values of ϵ_OFF(x, θr) are obtained on the basis of the integration of the spectral ROI recorded by the LaBr₃(Ce) detector [see Fig. 5 (a)]. These extracted values are further incorporated into the simultaneous algebraic reconstruction techniques (SART) algorithm [53], which is suitable for the reconstruction of high-quality images with limited observation angles. Fig. 6 shows images of -ln(ε_OFF, -ln(ε_ON), and -ln(ε_NRF) reconstructed by SART. The -ln(ε_OFF values decrease with the density of the materials; this behavior is similar to that of a conventional X-ray CT image. In addition, the signal of the uranium rod is clearly enhanced in the -ln(ε_ON) image, although the lead rod and air rod are still visible. In the -ln(ε_NRF) image, the signals produced by the air and lead rod do not appear, and the contrast of the uranium rod is higher than that in the -ln(ε_ON) image. These results demonstrate that the tNRF imaging method can be used to spatially discriminate a suspicious SNM isotope. The significance of the tNRF signals of ²³⁵U can be expressed in units of σ as Δ_tNRF = S/δS, where S and δS are the average $|-\ln ({\text{ε}}_{\text{NRF}})|$ value of four pixels in the uranium rod region and that of the remaining pixels, respectively. A tNRF imaging time of $t_{\text{tNRF}}^0$ 65 s is expected to yield a significance of ${\Delta }_{\text{tNRF}}^0$ = 8.0σ. This relatively long time is correlated with the photon flux, SNM concentration, and weak response of tNRF imaging to SNMs.

Figure 6:Full-screen View

(Color online)(a) Image of non-resonant atomic attenuation factor -ln(ε_OFF). (b) Image of on-resonant attenuation factor -ln(ε_ON). (c) Image of nuclear resonance attenuation factor -ln(ε_NRF). (d) Geometry and density of the interrogated target.

To evaluate the dependence of the attenuation factors on SNM density, we reconstructed CT images of ²³⁵U at artificial target densities of 14, 19, and 24 cm³. Reconstructed images similar to those in Fig. 6 were obtained. The average -ln(ε_ON), -ln(ε_OFF, and -ln(ε_NRF) values over the uranium rod region were then extracted, as shown in Fig. 7. One can see that the extracted values increase with increasing ²³⁵U density, which is consistent with the theoretical predictions. Subsequently, using Eq. 10, we also calculated the -ln(ε_NRF) values of tNRF images of the SNMs ²³⁸U, ^{239, 240}Pu, and ²³⁷Np (see Table 4). The expected -ln(ε_NRF) values for ²³⁸U (2468 keV) and ²⁴⁰Pu (2152 keV) are significantly larger than that for ²³⁵U (1734 keV), indicating excellent potential for the use of tNRF-CT for these SNMs. However, the -ln(ε_NRF) values for ²³⁹Pu and ²³⁷Np are smaller than that for ²³⁵U. This result suggests that the interrogations of both ²³⁹Pu and ²³⁷Np are questionable; a higher beam flux or longer imaging time would be required to obtain better results.

Figure 7:Full-screen View

(Color online) Dependence of attenuation factors -ln(ε_ON), -ln(ε_OFF, and -ln(ε_NRF) on the ²³⁵U density. The theoretical predictions given by Eqs. 8, 9, and 10 are also shown for comparison.