A new method for imaging nuclear threats using cosmic ray muons

Muon tomography is a technique that uses cosmic ray muons to generate three dimensional images of volumes using information contained in the Coulomb scattering of the muons. Advantages of this technique are the ability of cosmic rays to penetrate significant overburden and the absence of any additional dose delivered to subjects under study above the natural cosmic ray flux. Disadvantages include the relatively long exposure times and poor position resolution and complex algorithms needed for reconstruction. Here we demonstrate a new method for obtaining improved position resolution and statistical precision for objects with spherical symmetry.


Introduction
Cosmic ray muon scattering measurements provide a method for imaging the internal structure of objects using the information contained in the naturally occurring cosmic ray flux found at the surface of the earth [1][2][3]. The use of the signal from Coulomb scattering enables three dimensional imaging of complex scenes with high sensitivity to objects composed of materials with high atomic charge (Z). This makes this technique uniquely suitable for detecting and measuring the properties of atomic explosives that may be inaccessible with other techniques.

LA-UR-13-22644
If a nuclear material has been detected it is important to be able to measure details of its construction in order to correctly evaluate the threat. The flux of cosmic rays is low so it takes very long exposures to produce images with high resolution. In this paper we show how one can take advantage of spherical symmetry to improve the statistical precision of muon imaging.

Three Dimensional Imaging
Although cosmic rays are highly penetrating and can image through considerable overburden, the flux is limited. The times required detecting the presence of quantities of uranium or plutonium necessary to create a nuclear explosion are on the order of minutes, and the times needed to image these devices with ~2 cm resolution are on the order of hours.
Imaging with cosmic rays is based on measuring the multiple Coulomb scattering of the muons. The dominant part of the multiple scattering polar-angulardistribution is Gaussian: the Fermi approximation, where  is the polar angle and  0 is the mean multiple scattering angle, which is given approximately by:

2)
The muon momentum and velocity are p and  respectively, l is the material thickness, and X 0 is the radiation length for the material. This equation needs to be convolved with the cosmic ray momentum spectrum in order to describe the angular distribution.
There are several algorithms for generating tomographic images using the input and the output trajectories of the cosmic rays described by Schultz et al. [2] Here we use a very simple method with a scene composed of voxels . For each voxel a one-dimensional histogram of scattering angle in 1 mrad steps from 0 to 200 mrad is created. The histogram is incremented for every cosmic ray for which the incident cosmic ray intercepts the voxel and for which the distance between the input and the output cosmic rays in the plane of the voxel is less than 2 cm. For large scattering angles this requirement associates measured scattering events with defined voxels in which the most of the scattering occurred.
The scattering distribution for each voxel is fitted with a model that uses seven momentum groups, [4] p i , to approximate the muon spectrum with, 3) to approximate the muon energy distribution.
The model has been calibrated with data taken through three thicknesses of lead, 5.08, 10.16 and 15.24 cm. The amplitudes, A i , of each energy group, as well as the intrinsic angular resolution and a fixed number of radiation lengths due to the drift tubes and other structural elements of the muon detectors were fitted to minimize the logarithm of the likelihood function assuming the data were describe by a Poison distribution. This model does not account for changes in the shape of the muon spectrum due to stopping. A maximum likelihood fit to the set of lead data is shown in Error! Reference source not found.. Also shown is the decomposition of one of the data sets into its momentum groups. Images were constructed by fitting the angular distribution for each voxel to obtain the average number of radiation lengths of material that the ensemble of histogram entries has traversed. This reconstruction algorithm is scalable to large data sets, is simple to compute, and provides near optimal use of the scattering information. However, it doesn't optimally use the vertex information.
We have imaged several spherically symmetric objects: nested spherical shells of copper and tantalum, the copper shell alone, and a hollow lead ball. The outer radii were 6.5, 4.5, and 10 cm and the inner radii were 4.5, 1, and 2.5 cm for the copper, tantalum and lead LA-UR-13-22644 shells respectively. Cartesian slices through the three dimensional tomograms, centered on the object, are shown in Figure 2.

One Dimension Imaging
The center of each object (x c ,y c ,z c ) was estimated by finding the centroid of the signal from the object using the Cartesian slices shown in Figure 2. The were used for a onedimensional reconstruction. A trajectory was defined by the line x(s) where: Here (x 0 , y 0 , x 0 ) is a point on the line with direction cosines (x', y', x') . The point of minimum distance between (x c , y c , z c ) and (x 0 , y 0 , x 0 ) is given by s=s 0 : The radius of closest approach is: A two dimensional histogram of scattering angle versus radius is shown in Figure 3. A spherically symmetric object can be described by a set of shells at r i with thickness and of a material with radiation length, X 0i , and weighted density  vi /X 0i .
The radiation length weighted path length, L i , as a function of r i can be obtained from the data shown in the 2-dimensional histogram by using the multi-group fitting technique described above (Equation 3) for each radial bin. The fit to the lead data shown in Figure 1 has been corrected by 12% to account for the average 1/cos() increase in thickness of the lead in the planar geometry and then used here.
The L i are related to the set of radiation length weighted volume densities,  vj , by a path length vector P ij (see Figure 4): The path length vector is the length a particle at r i traverses through shell j:

8)
This can be solved for the  vj /X 0j using the regularization techniques described in Press. [5] This technique dampens the on-axis noise that arises in conventional Abel inversions [6] which is important here because of the poor statistics, especially at small r.
LA-UR-13-22644 The results for the three objects, each with 24 hours of exposure, are shown in Figure 5.
The  vj /X 0i for each of the materials studied here, copper, tantalum, and lead, are with 10% of the tabulated values. [7] One can easily distinguish the void inside each of the shells, even the 1 cm radius void at the center of the tantalum shell. An analysis of the width of the edges in Figure 5 give a position resolution of 3 mm. It is worth noting that these objects are difficult to study with conventional x-ray radiography. While the cavity can be observed with relatively crude collimation techniques, quantifying the cavity density takes sophisticated anti-scatter techniques. [8]

Conclusion
Cosmic ray scattering data were taken on a set of spherically symmetric objects. The data were analyzed assuming spherical symmetry. The data were stored in a two dimensional histogram of scattering angle vs. radius. The angular distributions were fitted by a sum of Gaussians whose amplitudes were fixed by fits to data taken on a set of planar objects. This resulted in one-dimensional plots of thicknesses in radiation lengths for each of the objects. These were analyzed with a regularized Abel inversion technique yielding radiation-lengthweighted volume densities. These results allowed a quantitative evaluation of the material composition of the objects.