Momentum-Dependent Cosmic Ray Muon Computed Tomography Using a Fieldable Muon Spectrometer

: Cosmic ray muon tomography has been recently explored as a non-destructive technique for monitoring or imaging dense well-shielded objects, classically not achievable with traditional tomographic methods. As a recent example of technology transition from high-energy physics to real-world engineering applications, cosmic ray muon tomography has been used with various levels of success in nuclear nonproliferation. However, present muon detection systems have no momentum measurement capabilities and recently developed muon-based radiographic techniques rely only on muon tracking. This unavoidably reduces resolution and requires longer measurement times thus limiting the widespread use of cosmic ray muon tomography. Measurement of cosmic ray muon momenta has the potential to signiﬁcantly improve the efﬁciency and resolution of cosmic ray muon tomography. In this paper, we propose and explore the use of momentum-dependent cosmic ray muon tomography using multi-layer gas Cherenkov radiators, a new concept for measuring muon momentum in the ﬁeld. The muon momentum measurements are coupled with a momentum-dependent imaging algorithm (mPoCA) and image reconstructions are presented to demonstrate the beneﬁts of measuring momentum in cosmic ray muon tomography.


Introduction
Monitoring large and dense objects, e.g., cargo containers or spent nuclear fuel dry casks, or large facilities, e.g., nuclear reactors, using cosmic ray muons could assess the stored content and provide information on whether it agrees with the declared content, maximize non-proliferation, and minimize potential security threats [1][2][3][4][5]. Conventional radiographic methods are limited, e.g., X-rays and gamma rays, and result in less than a few centimeters of penetration in lead [1]. Increasing the photon energy does not increase penetrating distance due to an increase in production of electron-positron pairs. Similarly, electrons cannot penetrate far into matter due to large radiation emission (Bremsstrahlung). For example, 1 GeV electrons can penetrate lead only within a few centimeters and commercially available portable electron accelerators can generate electrons with energies less than 10 MeV. Protons have been shown to penetrate several meters in dense materials and their potential has been demonstrated for large object imaging [2]. However, production of such subatomic particles requires large and immovable facilities, e.g., accelerators. Lastly, neutrinos can penetrate several kilometers in matter but their detection is impractical outside the laboratory, requiring kiloton size detectors.
A non-conventional class of ionizing particles, cosmic ray muons, presents certain advantages over the abovementioned charged and uncharged radiation. Cosmic ray muons are positively or negatively charged particles with~200 times the mass of electrons, generated naturally in the upper layers of the atmosphere, and reaching the earth with a mean energy of 3-4 GeV at a rate of~10,000 particles m −2 min −1 [6]. Muons can penetrate large and dense materials and their range increases with increasing energy. Figure 1 shows the range of muons as a function of atomic number and energy. Muons have the largest penetrating capabilities even in dense materials such as uranium and lead. In addition, muons are freely available and no radiological sources are required, resulting in a zero artificial radiation dose. Although the existence of cosmic muons and the theory which describes their transport and Multiple Coulomb Scattering (MCS) in matter is well understood, it is only recently that the applicability of muon scattering for imaging has been investigated [1,7]. Muon tomography is based on MCS and connects the expected muon scattering angle with the material atomic number and muon momentum. According to MCS, the distribution of muon scattering angles can be approximated using a Gaussian distribution and the characteristic parameters of the distribution can be determined from muon momentum and material properties. Currently, cosmic ray muon tomography is solely based on scattering angle measurement without the capability for muon momentum measurement. Recent efforts at INFN have developed a large prototype based on MCS coupled with muon momentum reconstruction algorithms [8]. Although these efforts show promise, no portable detectors exist that can measure muon momentum in the field. This requires longer measurement times to compensate for the reduction in image resolution. To address this limitation, a Cherenkovbased muon spectrometer is currently under development for measuring muon momentum in the field [9,10]. Using pressurized gas Cherenkov radiators and depending on the pressure level, we can create a set of increasing muon momentum thresholds. Each radiator is triggered only when the incoming muon momentum exceeds the muon momentum threshold level. As a result, the signals from each radiator depend on the incoming muon momentum as not all Cherenkov radiators can generate a signal. By analyzing which radiator is triggered, we can then estimate the actual muon momentum. In addition, we extend the widely used reconstruction algorithm, Point of Closest Approach (PoCA) to include momentum into a new momentum dependent PoCA (mPoCA) algorithm. In the end, results comparing the two image reconstruction algorithms, PoCA and mPoCA, are presented to show the benefits of momentum measurement in cosmic ray muon applications.

Overview of Muon Tomography Detection Capabilities
Cosmic ray muon tomography is based on the principle of multiple Coulomb scattering (MCS). The estimated muon scattering angle distribution can be approximated using a Gaussian distribution with zero mean and standard deviation [11]: where θ is the muon scattering angle, σ θ is the width of Gaussian scattering angle distribution, X is the thickness of the scattering medium, and X 0 is the radiation length. β µ is the ratio of the muon's speed to the speed of light (c) and p µ is the muon momentum. Early efforts have mainly focused on muon detectors for cargo scanning applications and proof of principle studies [12][13][14][15] and several detector types have been proposed for muon tomography. Existing detectors measure the locations and directions of incoming and outgoing muons and include drift tubes, scintillators, resistive plate chambers, and solid state detectors. Muon detectors are comprised of multiple parallel planes of positionsensitive chambers although cylindrical geometries have also been proposed [16][17][18][19]. A typical arrangement is shown in Figure 2. Cox et al. [20] summarized the requirements for muon detector development. Spatial resolution in the order of sub-mm, coincidence timing in the order of nanoseconds, and picosecond timing resolution would be required to differentiate between sub-GeV and GeV muon energies for future applications. They noted that efficient detectors with a large area and sensitivity over hundreds of GeV would be needed to optimize statistics and reduce the time taken to collect an adequate number of muons. Muon events are recorded only when muons pass through all detector planes within a time window of a few nanoseconds. Only muon trajectories that fall within the superimposed grid boundaries are useful for reconstruction. Density maps can be generated with millimeter to centimeter spatial resolution and different materials can be discriminated depending on their atomic number using reconstruction algorithms. The discrimination capability and image quality depend on muon flux and are usually limited by incomplete information about muon path and momentum [21]. Minimizing the number of muons provides substantial time savings but can worsen the image quality, whereas a higher number of measured muons can improve image contrast and resolution but also increases measurement time. Recent efforts have shown that muons in the order of millions are required for imaging large-scale objects such as shipping containers [17,22,23] or dry casks [19,24,25]. This translates to measurement time in the order of hours or even days, depending on the detector characteristics, i.e., size and orientation. Table 1 provides a summary of existing muon detectors.  [28] Plastic scintillator strips 3 × 6 m 2 3 × 6 × 2.8 m 3 3.5 mm (3.5 mrad) >90% INFN [29] Drift wire chambers 3 × 2.5 m 2 3 × 2.5 × 1.6 m 3 200 µm >90% CRIPT [16] Plastic scintillator bars 2.0 × 2.0 m 2 1.6 × 2.0 × 2.0 m 3 3 mm (6 mrad) >99.5% RPC [22] Resistive plate chambers 0.5 × 0.5 m 2 0.09 × 0.5 × 0.5 m 3 300-900 µm >95% The PoCA image reconstruction algorithm has been widely used in cosmic ray muon tomography with applications in cargo scanning, spent fuel cask imaging, and reactor monitoring. The algorithm assumes that muons are scattered at a single point. This ignores multiple Coulomb scattering resulting in resolution limitations and image blur. The PoCA algorithm, due to its minimal computing requirements, is typically applied for reconstructing large objects with thousands of voxels, within a reasonable amount of time. PoCA assigns a point at the minimum distance between the incoming and outgoing trajectories. The intersection of the two trajectories is described mathematically as follows: where p, Q are points on the incoming trajectory L 1 and the outgoing trajectory L 2 , respectively. In 3D space, where trajectories may not intersect, the algorithm calculates the middle point of the shortest line between trajectories and the scattering angle, θ, is assigned to a voxel on a 3D grid. At the estimated point of scatter the following value is assigned: where θ x and θ y are the projected angles onto an x-y plane with respect to a positive xand y-axis, respectively ( Figure 2). As more muons pass through the object volume, each pixel (or voxel) is filled with scattering angles. The final value of each voxel is the average scattering angle divided by the voxel size L: where N is the number of scattering angle values in a voxel. Dividing by the voxel size is needed to normalize the voxel value for geometries where the voxel size is not the same everywhere. Materials with a higher atomic number will generate higher voxel values as a result of increased scattering in this location. The computational complexity of the algorithm depends on the number of muons passing through the object volume and scales with O(N) [24]. Other reconstruction algorithms that have been explored for use in cosmic ray muon tomography belong to the class of filtered back-projection algorithms (FBP). When a large number of evenly spaced projections is not available, a necessary condition for FBP, or significant scattering takes place along the ray path, a solution can be provided using algebraic reconstruction techniques (ART). ART is computationally expensive but conceptually simpler than FBP. In this algebraic approach, a grid is superimposed on top of the object. Each grid pixel contains the unknown object function and a set of algebraic equations is set up for the determination of the unknown function. A necessary condition for ART is the correct estimation of the ray paths between the source and detector. Simulation and modeling in cosmic ray muon tomography is typically undertaken using GEANT4 (GEometry ANd Tracking) [30][31][32][33], a Monte Carlo code for transport of high energy particles in matter developed at CERN. As schematically shown in Figure 3, GEANT4 requires parameters related to object geometry, muon energy, and angular distributions.

Fieldable Muon Spectrometer and Momentum-Dependent PoCA
Recent simulations showed that a significant number of muons (>10 6 muons) would be needed to adequately image large objects, e.g., identify missing spent nuclear fuel in dry casks [19,24,25]. This was confirmed by a recent experimental effort by LANL [26] which resulted in~200 h of measurement time in order to register only 1.62 × 10 5 muons, about 100 times lower than theoretically possible. The simulated results did not account for any additional measurement time needed to address noise and/or detector measurement uncertainties. Minimizing time measurement requirements would require maximizing the number of muons used for imaging to muons measured, i.e., imaging efficiency. This motivates the use of muon momentum as an additional feature to improve cosmic ray muon tomography beyond what is currently feasible.

Fieldable Muon Spectrometer
When the velocity of a muon surpasses the speed of light in an optically transparent medium, or a radiator, Cherenkov light is emitted (Figure 4). Therefore, there is a threshold muon momentum level for muons to induce Cherenkov radiation in a given radiator and it is derived as follows [34]: where p th is the threshold muon momentum, n is the refractive index of the medium, and m µ c 2 is the muon rest mass energy. In gas radiators, the refractive index depends on pressure and temperature and it can be approximated by the Lorentz-Lorenz equation: where A m is the molar refractivity, p is the gas pressure, T is the absolute temperature, and R is the universal gas constant. Therefore, the threshold muon momentum is given by: Equation (8) shows the relation between the Cherenkov threshold momentum for muons and gas pressure. We choose CO 2 because it can be pressurized (up to 5.7 MPa) without condensation. Continuous threshold momentum levels with an interval of 1.0 GeV/c are designed by pressurizing CO 2 gas radiators, p th = 1.0, 2.0, . . . , 5.0 GeV/c ( Figure 4). We used a SiO 2 solid radiator to achieve the lowest threshold momentum level (p th = 0.1 GeV/c) because it is not possible to obtain with pressurized CO 2 at room temperature. As a result, a signal will be generated only when the momentum threshold is less than the actual muon momentum. This applies even in the case where a muon interacts with all radiators. The actual muon momentum can be estimated by measuring the Cherenkov signals generated or not generated in each radiator. Overview of Cherenkov muon spectrometer using one solid (SiO 2 ) and five pressurized CO 2 radiators (top) [35]. The expected yields of Cherenkov and scintillation photons of each radiator (italic numbers represent radiator IDs) as a function of muon energy (lower left) and their wavelength distribution (lower right) are shown [35].

Momentum-Dependent PoCA
To include momentum in a PoCA algorithm, the measured muon momentum is additionally tagged in a PoCA point along with a scattering angle and this momentumdependent PoCA imaging algorithm is named mPoCA. Therefore, each muon event includes five values: 3D Cartesian components of a PoCA point (P x , P y , and P z ), the reconstructed spatial muon scattering angle, θ space , and the measured muon momentum, p µ , where N µ is the total number of muon events. All measured muon momenta are classified by a finite number of momentum levels depending on the number of radiators in the muon spectrometer. Therefore, X i (Equation (9)) can be decomposed with respect to momentum groups: where L X is the length of X j , N p is the number of muon momentum groups, and V lmn is the PoCA voxel when l, m, and n are integers. Because the expected muon scattering angle distribution depends on muon momentum, each X j T (j = 1, 2, . . . , N p ) will have a unique scattering angle distribution. This decomposition process was not possible without muon momentum measurement.

GEANT4 Simulations
The feasibility of the proposed fieldable muon spectrometer was evaluated using GEANT4. To generate muons from the actual measured muon spectrum, a "Muon Event Generator" is used based on a phenomenological model that captures the main characteristics of the experimentally measured spectrum coupled with a set of statistical algorithms [36,37]. The muons generated have zenith angles in the range 0-90 • and energies in the range 1-100 GeV. The "Muon Event Generator" generates the necessary user-defined histograms for use with the GEANT4 G4GeneralParticleSource macro file. Hence, a stochastic muon transport simulation using GEANT4 was performed to evaluate the characteristics of the proposed muon spectrometer.
To measure muon momentum, we modeled a Cherenkov muon spectrometer in GEANT4 which consisted of one glass (1 cm), five CO 2 gas radiators (10 cm), and thin photonabsorbing black-aluminum foils (0.1 cm) to minimize photon losses to the environment. To measure the produced optical photons, an array of silicon photomultiplier (SiPM) was placed in each radiator. The overall length was 51.7 cm and the surface area 0.2 × 0.2 m 2 . SiO 2 with 1.45 refractive index was used for the glass radiator. Pressurized CO 2 gases were used as gas radiators to provide sequentially increasing refractive indices and increasing momentum threshold levels. Some CO 2 gas radiators needed to be depressurized to achieve p th = 4 and 5 GeV/c. Any optical photon events were terminated once they arrived either at the photon absorbers or outside of "world" volume during the simulations. All major physics are included in the GEANT4 reference physics list, QGSP_BERT [33].

Results
A proposed overview implementation in nuclear safeguards is shown in Figure 5. The spectrometer is installed between objects and muon trackers to measure momentum for all muon events. Because the performance of a muon spectrometer depends on the number of radiators, we assessed the performance by increasing the number of radiators to 10 and 100 radiators and extending the momentum range from 0.1 to 10.0 GeV/c. We performed four simulations with a different number of radiators using 10 4 muon samples in each case: (i) 10 2 radiators (fine measurement resolution, ±0.05 GeV/c) without noise, (ii) 10 2 radiators with noise (i.e., including scintillation and transition photon emission), (iii) 10 radiators (coarse resolution, ±0.5 GeV/c) without noise, and (iv) 10 radiators with noise. When the instrumental uncertainty is not considered, the cosmic ray muon spectrum was successfully reconstructed using 10 2 radiators ( Figure 6). Although the increased number of radiators improved resolution, it negatively impacted the signal-to-noise ratio due to the decreased Cherenkov signals. When noise was added, the spectrum was slightly shifted to the right because noise (additional photon signal) increases the probability of false classification. We could improve the accuracy by suppressing the noise level because the momentum overestimation pattern was predictable. For the simulation results with 10 radiators, the cosmic ray muon spectrum was adequately reconstructed even if there were fewer radiators.  To statistically evaluate our muon spectrometer, we calculated the classification rate (CR), which is the probability that the system correctly identifies the actual muon momentum. The computed classification rates using 10 4 mono-energetic muons for various discrimination levels to remove noise, i.e., 0, 1, and 2, are shown in Figure 6. These linear discriminators were used to remove any noise signal from the total signal. When muon momentum was less than 1.0 GeV/c, the overall CR was less than 60% because the scintillation photons dominated in the glass radiator. When 1.0 < p th < 4.0 GeV/c, all gas radiators showed high CRs because the predictable noise level (mostly from the scintillation) was discriminated out from the total signal. This situation was further improved with a linear discriminator. However, when 4.0 < p th < 6.0 GeV/c, CRs began to decrease. This is attributed to low Cherenkov photon yields due to the depressurized CO 2 (~a half atmosphere pressure) gas. When p th > 6.0 GeV/c, on the other hand, CR rebounded because Cherenkov photon yields gradually increase as muon momentum increases. However, it is not efficient to use any linear discriminator in high p th (>5.0 GeV/c).
To explore the benefits of muon momentum measurement in nuclear safeguards, e.g., special nuclear material (SNM) monitoring, three spherical special nuclear material samples, i.e., significant amounts (SQ) of high-enriched uranium (HEU), low-enriched uranium (LEU), and plutonium (Pu), were studied using Monte Carlo simulations where they were surrounded by various thicknesses of lead shielding. 5 × 10 3 muons were used to build scattering angle variance distributions [38,39]. In Figure 7, the mean scattering angle variance of distributions for HEU, LEU, and Pu are plotted as a function of thickness of lead shielding with 2σ errors when the muon momentum follows a cosmic ray muon spectrum with p µ = 0.2-100.0 GeV/c (left column) and p µ = 3 GeV/c (right column). Without measuring muon momentum, the scattering angle variance distributions for the three types of SNMs begin to overlap as the thickness of lead shielding increases. On the other hand, they are clearly distinguishable with muon momentum knowledge even when the thickness of lead shielding is 30 cm. The following simulation results show the reconstructed images of spherical LEU samples surrounded by 5 cm thick lead shielding (Figure 8). We used 10 5 muon samples in GEANT4 to simulate Monte Carlo muon transports and muon tomographic images were reconstructed using a (i) PoCA algorithm with a cosmic ray muon spectrum when p µ = 0.2-100.0 GeV/c, (ii) PoCA algorithm with mono-energetic muon, p µ = 3 GeV/c, and (iii) mPoCA algorithm where materials were classified in three groups, M1, M2, and M3, for SNMs, shielding materials, and others, respectively. In the mPoCA algorithm, the most probable material group was recorded instead of scattering angle values (Equation (5)) in each voxel. The results show that the imaging resolution is significantly improved and boundaries between SNM (LEU) and lead shielding (5 cm) are visually distinguished when the number of noise voxels in the air gap decreases.

Conclusions
Cosmic ray muons are relativistic particles with energies in the order of GeV that can play an important role in non-proliferation and nuclear material monitoring applications because they present unique advantages over existing methods. Research on muon detectors and imaging approaches for real-time monitoring and imaging of large objects using muon scattering, displacement, and transmission are constantly improving and expanding the cosmic ray muon applications in various fields. In the present work, we evaluated a Cherenkov muon spectrometer using multi-layer pressurized gas Cherenkov radiators. Six sequential threshold momentum levels were set using a solid radiator (SiO 2 ) and pressurized CO 2 gases, p th = 0.1, 1.0, 2.0, 3.0, 4.0, and 5.0 GeV/c. Then, we performed detailed GEANT4 simulations to evaluate the functionality of the proposed fieldable muon spectrometer. We were able to measure muon momentum with a resolution of ±0.5 GeV/c (a mean CR is 87%) within a momentum range of 0.1 to 10.0 GeV/c. We showed the potential of measuring muon momentum in nuclear nonproliferation, e.g., special nuclear material monitoring, by comparing muon scattering distributions with and without muon momentum measurement. Finally, we extended the widely used reconstruction algorithm, PoCA, to include momentum into a new momentum-dependent PoCA (mPoCA) algorithm. With the mPoCA, imaging resolution improved and noise decreased by a factor of 3 to 4.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author upon reasonable request.