1 Introduction

Neutron spectrometry has been extensively studied [1-3]. A variety of neutron spectrometry techniques were reviewed by Brooks et al. [4], and were classified into seven groups based on the principles employed, namely, (1) measurement of energies of recoiled nuclei to infer about the properties of scattered neutrons, (2) measurement of energies of charged particles produced during neutron-induced nuclear interactions, (3) measurement of neutron velocities, (4) threshold methods based on determination of minimum neutron energy through identification of neutron-induced effects such as radioactivity, a specific gamma-ray energy or a phase transition, (5) employment of results from neutron-diffraction-based methods, (6) measurement of time-distribution of slowing down of high-energy neutrons in a medium and (7) determination of neutron-energy distribution via unfolding the readings from the detector based on energy-dependent responses to neutrons.

The Bonner Sphere Spectrometer (BSS) was first introduced by Bramblett, Ewing and Bonner [5]. In this system, a thermal neutron detector is enclosed by a series of spheres made up of hydrogenous materials acting as neutron moderators. To simplify the storage and setup [6], a new nested-shell system known as the Nested Neutron Spectrometer (NNS) (Detec, Quebec) was designed [7]. The NNS consisted of seven cylindrical high-density polyethylene (HDPE) shell moderators assembled in a Russian doll fashion. For measurements with the NNS, the count rates recorded by the bare detector and those recorded upon addition of each HDPE shell moderator, i.e., a total of eight, were obtained. Since the detector’s energy-dependent responses to neutrons upon addition of each shell moderator are known, the neutron energy spectrum (NES) can be determined from these count rates using "spectrum unfolding".

For NNS applications, the only spectrum unfolding tool commercially available is the NNS2.8.1 software supplied by the NNS vendor, which is based on the STAY’SL algorithm [8] that uses the least squares (LS) method to unfold the NES. However, the LS method necessitates a relatively close-to-real guessed spectrum to produce a correct unfolded NES [9], which will present difficulties under scenarios when the NES is entirely unknown, such as those for photoneutrons created in radiotherapy rooms, neutrons generated in nuclear reactors or unwanted neutrons created in spent fuel [10]. In view of these difficulties, Maglieri et al. [6] used the maximum-likelihood expectation-maximization (MLEM) method to unfold the NES at various locations in the radiotherapy room. It was correctly pointed out by Maglieri that the MLEM method did not require close-to-real guessed spectra to be able to generate a correct NES in contrast to the LS method [9,11]. However, the effect of initial guess spectra on the output unfolded NES was not studied. Further studies would be pertinent, and in particular would provide the users information on the effects of input variables on the final unfolded NES.

2 Materials and methods

As described above, the NNS consisted of seven HDPE shell moderators assembled around the central ³He probe. The ³He probe was a cylindrical rod-shaped proportional counter containing 2 atm of ³He. Table 1 summarizes the dimensions of the PE moderators. Addition of each HDPE shell moderator will lead to a different energy-dependent response of the detector to neutrons and thus a different count rate, so in total eight different count rates will be obtained from the eight configurations of the NNS, which are schematically depicted in Fig. 1. Our task is therefore to reconstruct the NES from these eight different count rates.

Fig. 1Full-screen View

Count rate measurements under eight different configurations of the NNS with (a) no, (b) 1, (c) 2, (d) 3, (e) 4, (f) 5, (g) 6 and (h) 7 high-density polyethylene shell moderators.

The response function for the NNS system comprising eight HDPE shell moderators was obtained specifically for the current design (i.e., cylindrical geometry). The spatial dependencies that might arise from the cylindrical geometry were already taken into account when the response functions were determined. Since the method described in the present work focused on the present geometry and design, it was not necessary to further the spatial dependencies.

2.1 Maximum-likelihood expectation-maximization (MLEM) method

The NNS consisted of seven circular-cylindrical high-density polyethylene (HDPE) shell moderators assembled in a Russian doll fashion. The ³He detector was accommodated at the center of the innermost shell moderator. Fig. 2 schematically shows the experimental setup.

Fig. 2Full-screen View

Schematic diagram of cross-sections of the NNS system, showing seven circular-cylindrical high-density polyethylene (HDPE) shell moderators (not to scale), with ³He detector (not to scale) accommodated at the center of the innermost moderator.

We denote the response function of the detector and the counts experimentally recorded by the detector upon adding the i^th moderator as R_i(E) and m_i, respectively, where E is the neutron energy. The response function was obtained from the vendor for the design and geometry of the NNS system. However, here, the response functions from the manufacturer was used, while the NES represented by n(E) was the unknown to be found through unfolding. The mathematical relationship among m_i, R_i(E) and n(E) was governed by the well-known Fredholm integral of the first kind:

$m_i\left(E\right)=\underset{E_{\text{min}} }{\overset{E_{\text{max}}}{{\displaystyle \int }}}{R}_i\left(E\right)n\left(E\right)\text{d}E$ (1)

where R_i(E) was referred to as the kernel. Eq. (1) can be written in a simplified form as

$m=Rn$ (2)

where m and n represent the data and solution functions, respectively, and R is a linear operator. The MLEM method maximizes the likelihood of obtaining the NES when convergence is achieved. To achieve this, the iterative form of the MLEM method is used as

$n_j^{k+1}=\frac{n_j^k}{{{\displaystyle \sum }}_{i=1}^N{r}_{ij}}{\displaystyle \sum }_{i=1}^N{r}_{ij}\frac{m_i}{{{\displaystyle \sum }}_{b=1}^J{r}_{ib}{n}_b^k}$ (3)

where n_j^k and n_b^k are the current and starting spectra, respectively, and r_ij and m_i are the response function of the detector and the count measurements, respectively. To perform the unfolding and to determine the Poisson fluctuations associated with the unfolded NES, a flexible yet fast numerical algorithm is required, which is shown below:

1: Start NRUunfold code

2: load measured counts: m, initial spectrum: n_b, cutoff: N_max, NNS normalization factor: f_nns

3: initialize i

4: while i < N_max, do step 5 to 12

5: compute d = transpose(m)×f_nns

6: determine a = d/(R×n_b)

7: compute c = transpose(R)×a

8: construct 8×1 array of ones: u₁

9: update the initial spectrum, n_b = (n_b•×c)/(transpose(R)×u₁)

10: if a = u₁ then break

11: update i, i = i+1

12: end

13: obtain the final spectrum, S_out

14: construct 52×52 array of zeros: z₁

15: for j = 1,2,…,N_spect

16: generate pseudo measurement (m_pseudo) from Poisson distribution with mean: d

17: initialize k

18: while k < N_max, do step 19 to 26

19: compute d = transpose(m_pseudo)×f_nns

20: determine a = d/ (R×n_b)

21: compute c = transpose(R)×a

22: construct 8×1 array of ones: u₁

23: update the initial spectrum, n_b = (n_b•×c)/(transpose(R)×u₁)

24: if a = u₁ then break

25: update k, k = k+1

26: end

27: obtain the final pseudo spectrum, S_pseudo

28: end

29: compute spectra variance, v = (1/52)×sum[(s_pseudo – s_out)²]

30: compute standard deviation, σ = √v

Please note that the operation (•×) represent the scalar multiplication, whereby (×) represent the matrix multiplication. The steps 5 to 12 in the presented algorithm represent the standard MLEM method, in which the ratio of data to convolution of changing guess spectrum is computed through the use of response matrix, which means the change in the convoluted spectrum is bounded by the response matrix. Here the measured counts (m) was an array with dimension of [1×8], the matrix d had a dimension of [8×1] and the supplied initial guess spectrum (n_b) was a matrix with dimension of [52×1]. The normalization factor (f_nns) was supplied by the manufacturer to correct for the response function (also provided by the manufacturer and used in the present work). The response functions provided by the manufacturer could be found in Ref. [6], which were employed in the present work for different moderators (labeled with different numbers) as shown in in Fig. 3. These response functions were obtained specifically for the design of the current NNS system. Here, we employed the 52-bin energy spectrum to facilitate comparisons between our results and those obtained using the STAY’SL algorithm, the latter also having used the 52-bin energy spectrum. It is remarked here that the response functions were supplied by the vendor in the form of a 52-bin energy spectrum, so changing the energy spectrum will require re-computation of the response functions and will lead to complications, which are out of the scope of the present work. Large variations in the unfolded spectra were expected when a 52-bin energy spectrum was derived from 8 measured count rates. Here, we used an error constraint of 0.01 in addition to the hard iteration cutoff (proposed in ref. [6]) to damp these variations in the unfolded neutron energy spectra.

Fig. 3Full-screen View

Response functions for different moderators (labeled with different numbers) for the Nested Neutron Spectrometer (see Ref. [6]).

2.2 Initial spectra

The main benefit of using the MLEM method (over the LS method) is its non-reliance on the initial guessed spectrum. The convergence of our MLEM algorithm was tested for a variety of initial guessed spectra obtained for the present 52-bin energy unfolding code, including (1) step function, (2) uniform function, (3) chopped-sine function, (4) exponential function, (5) absolute value of the Dirichlet function and (6) random noise spectrum. The Dirichlet function is often referred to as the "periodic sinc" or "aliased sinc" function:

$\left|D\left(x\right)\right|=\{\begin{array}{c}\frac{\sin \left(\frac{Nx}{2}\right)}{Nsin\left(\frac{x}{2}\right)}, x\ne 2\pi k\\ {\left(-1\right)}^{k\left(N-1\right)}, x=2\pi k\end{array}$ (4)

These spectra are shown in Fig. 4.

Fig. 4Full-screen View

Initial guessed spectra supplied to the NRUunfold code in the present study: (a) step function, (b) uniform function, (c) chopped-sine function, (d) exponential function, (e) absolute value of Dirichlet function and (f) random noise spectrum. (Note: the energy structure is based on 52 bins of equal size on a logarithmic scale, ranging from 1.12×10^-9 to 15.9 MeV).

3 Benchmarking of NRUunfold code

To benchmark the present model, we used the measured count rates supplied by the vendor for the heavy water moderated ²⁵²Cf fission neutron source [7], which were previously unfolded using the STAY’SL program [8]. The NRUunfold code generated spectra compatible with the STAY'SL output for most initial spectra, but some exceptionally bizarre initial guess (i.e., Dirichlet and random noise) gave unreasonable output spectra as shown in Fig. 5. The unfolding task involving the random noise spectrum was more tedious mainly due to its stochastic nature. The obtained results showed the effect of initial guess spectra on the output unfolded NES. The peak locations appeared somewhat better at lower energies and some (not very significant) shifts were noticed at higher energy peaks for the ²⁵²Cf neutron source, especially for very complex initial guess spectra such as the Dirichlet spectrum. This could be attributed to the large number of features in the heavy water moderated ²⁵²Cf neutron energy spectrum (when compared with, for example, ²⁴¹Am-Be) since it only used 7 to 8 bins of the 52 bins in the current energy range, which made the unfolding task tedious and in turn introduced some deviations. Furthermore, through employment of a 52-bin energy grid, such deviations were expected for complex neutron energy spectra such as the one from a heavy water moderated ²⁵²Cf neutron source since a larger number of bins would be used for unfolding these neutron spectra.

Fig. 5Full-screen View

Comparison among the unfolding results from STAY’SL and form the present NRUunfold code with different initial spectra for the heavy water moderated ²⁵²Cf fission neutron source. (Note: the energy structure is based on 52 bins of equal size on a logarithmic scale, ranging from 1.12×10^-9 to 15.9 MeV).

The energy ranges of thermal fluences and fast fluences for the heavy water moderated ²⁵²Cf source were determined to be 10^-8-10^-7 and ~1-10 MeV, respectively. The main benefit of the proposed algorithm was that the obtained spectra were comparable to that supplied by the software of the vendor, with the advantage of not being biased by the pre-knowledge of the spectrum to be calculated.

4 Computation scheme

Fig. 6 shows the normalized count rates from the Am-Be neutron source obtained using the response functions supplied by the manufacturer. These normalized count rates were used for sensitivity studies on NES unfolding (remark: these were not experimentally obtained). The normalized unfolded NES for the Am-Be neutron source was scored for a variety of initial guessed spectra. The NES unfolding was performed for various count rates obtained through multiplication of an integer n with the normalized Am-Be count rates (n×Φ_norm) shown in Fig. 6. The value of n was set as 100 for all initial guessed spectra shown in Fig. 4. The variations in the final unfolded NES and Poisson fluctuations were studied for each case and for two different iteration cutoffs (i.e., 10,000 and 20,000 cycles).

Fig. 6Full-screen View

Normalized count rates from the Am-Be neutron source obtained from the response functions supplied by the manufacturer.

5 Experimental measurements

The energy spectrum of a 10 mCi ²⁴¹Am-Be neutron source was experimentally determined in the radiation laboratory of the Nuclear Radiation Unit (NRU), Department of Physics, City University of Hong Kong. The geometry of the radiation laboratory is shown in Fig. 7 (with the source and NNS location shown). The concrete in the walls, ceiling and floor were the potential scattering sites for the neutron spectrometry measurements. The geometric relationship between the sensitive detector volume and the source was kept constant for all employed configurations with different shells. Interference from γ-rays was assumed to be negligible and not considered in the present measurements.

Fig. 7Full-screen View

Geometry of radiation laboratory in which the energy spectrum of a 10 mCi ²⁴¹Am-Be neutron source was experimentally determined using a Nested Neutron Spectrometer (NNS).

6 Unfolding results and discussion

The normalized unfolded NES obtained using various initial guessed spectra for the Am-Be source are shown in Fig. 8. The NRUunfold program could produce decent NES for most of the employed initial guess spectra, but some exceptionally bizarre initial spectra (Dirichlet and random noise) gave unreasonable output spectra. In contrast, the LS method employed in the unfolding software supplied by the NNS vendor (named NNS2.8.1) heavily relied on the supplied initial guessed spectrum.

Fig. 8Full-screen View

Comparison among unfolding results for Am-Be neutron source using NRUunfold code with six different initial spectra and using STAY’SL with a very precise guessed spectrum. Note that STAY’SL did not converge using the six initial spectra supplied to NRUunfold. (Note: the energy structure is based on 52 bins of equal size on a logarithmic scale, ranging from 1.12×10^-9 to 15.9 MeV).

The Chi-squared statistical test was used to assess the goodness of fit of the unfolded neutron energy spectra (using different initial guessed spectra) to the true spectrum. In the present Chi-squared test, the uncertainty associated with each bin was considered. The critical value for the test at p < 0.05 was 68.67. The results of the Chi-squared test for unfolding ²⁴¹Am-Be and ²⁵²Cf spectra are shown in Fig. 9.

Fig. 9Full-screen View

Chi-squared values for unfolding ²⁴¹Am-Be and ²⁵²Cf spectra from different initial guessed spectra with uncertainties considered.

The results obtained using the Chi-squared statistical test demonstrated that the unfolded spectra were compatible with the true spectrum. Nevertheless, the unfolded NES obtained using the Dirichlet and random noise initial guess spectra were unreasonable considering some important NES features (see Fig. 5 and 8). As such, the result spectra should be examined to judge whether the unfolding was satisfactory.

Comparisons between the unfolded NES and associated Poisson fluctuations for the normalized Am-Be counts (see Fig. 7) at 100×Φ_norm with 10,000 and 20,000 cutoff iterations for the six-different initial guessed spectra are shown in Fig. 10. The shape of the unfolded NES accurately described the neutron energy distribution emitted from the Am-Be source. Moreover, different initial guessed spectra led to similar unfolded NES for the same iteration cycle cutoff. Increasing the iteration cycle cutoff from 10,000 to 20,000 cycles slightly increased the Poisson fluctuations in the unfolded NES but the shape and the fluence rates were similar for both iteration cutoffs. It is remarked that both small number of counts (i.e., poor statistics) and large number of iterations significantly increased the noise in the final unfolded NES. We noted that the noise associated with unfolded neutron energy spectra would become noticeable for more than 20,000 iterations. Depending on the measured count rates, the introduced noise would cause fluctuations in the unfolded neutron energy spectra at various energies, which led to deviations of unfolded neutron energy spectra from realistic spectra.

Fig. 10Full-screen View

Unfolded NES with associated Poisson uncertainties (shaded in blue) at different neutron energies using (a) step, (b) uniform, (c) chopped-sine, (d) exponential, (e) Dirichlet functions and (f) random noise spectrum as initial guessed spectra with different iteration cycle cutoffs. (Note: the energy structure is based on 52 bins of equal size on a logarithmic scale, ranging from 1.12×10^-9 to 15.9 MeV).

Furthermore, the present computer program was used to unfold the experimentally obtained count rates using the NNS as shown in Table 2 for a 10 mCi Am-Be neutron source. The computed Am-Be count rates obtained using the manufacturer’s response functions were used for NES unfolding under the ideal situation, where no uncontrollable factors were involved in contrast to realistic experimental measurements. The Chi-squared statistical test was then carried out to ensure the goodness of the model in unfolding the neutron energy spectra: this provided information for assessing the accuracy of the present model before using the current algorithm for practical applications. We employed 10,000 iteration cycles, and the NNS normalization factor (f_nns) of 1.02 provided by the vendor. The NES unfolded using NRUunfold code with six different initial spectra are shown in Fig. 11 for different initial guessed spectra.

Fig. 11Full-screen View

NES unfolded using NRUunfold code with six different initial spectra with associated Poisson uncertainties (shaded in blue) at different neutron energies for a 10 mCi Am-Be neutron source. NES unfolded using STAY’SL with a very precise guessed spectrum is also shown for comparison. Note that STAY’SL did not converge using the six initial spectra supplied to NRUunfold. (Note: the energy structure is based on 52 bins of equal size on a logarithmic scale, ranging from 1.12×10^-9 to 15.9 MeV).

A large peak at E > ~1 MeV was displayed, which was the main energy peak for the Am-Be neutron source (cf. Fig. 8). Two relatively small peaks in the energy range from 10^-2 to 1 MeV and from 10^-8 to 10^-6 MeV were observed, which were likely due to reduction in neutron energies caused by neutron scattering in the laboratory. The peak in the range from 10^-8 to 10^-6 MeV could be referred to as the "thermal peak" and capturing of this peak demonstrated that the present computer program could reveal the neutron energy distribution ranging from thermal to fast neutrons. The potential biases from room scattering might be attributed to the concrete in the walls, ceiling and the floor. In contrast, the STAY’SL algorithm gave a very large thermal peak in the energy range from 10^-8 to 10^-6 MeV, with peak fluence even higher than that for the main energy peak (E > ~ 1 MeV). On the other hand, the Poisson fluctuations were more substantial when compared with the results (100×Φ_norm) shown for the Am-Be source in Fig. 10 for 10,000 or even 20,000 iteration cycles, which demonstrated that the low number of neutron counts significantly increased the uncertainty of the unfolded NES even for low iteration cutoffs.

7 Conclusion

In the present work, a generic and versatile MLEM algorithm was introduced, which could produce decent unfolded NES for most initial guess spectra. We tested six different initial guessed spectra, i.e., (1) step function, (2) uniform function, (3) chopped-sine function, (4) exponential function, (5) Dirichlet function and (6) random noise spectrum. Remarkably, decent unfolded NES neutron spectra were obtained for most of the guessed spectra. However, for the discontinuous Dirichlet function and the random noise spectrum, the unfolded NES were rather unreasonable, which suggested that an exceptionally bizarre initial guess might not lead to a problem. The present model was tested using neutron count rates generated from the manufacturer’s response functions and also those experimentally measured for a 10 mCi Am-Be neutron source. Our results also revealed the effects of input variables (e.g., guessed spectra, count rates and number of iterations etc.) on the final unfolded NES. The program is particularly useful when there is no prior knowledge on the NES, such as those for photoneutrons created in radiotherapy rooms and neutrons generated in nuclear reactors. Although NNS was used to demonstrate the feasibility of the method, the concept was applicable to the BSS and all nested-shell systems including NNS. The current program will also be useful for further enhancement and optimization of nested-shell systems, e.g., in designing the thickness of shells of these nested-shell systems. The present method can be used in various applications and it has the potential to alternate conventional unfolding methods which need precise initial guesses. However, more exhaustive studies should be carried out to confirm the general applicability of the method, e.g., for innumerable neutron spectra for wide energy ranges (thermal to GeV) in the world.