Under extreme conditions of high temperature and high pressure, the behavior of materials becomes extremely complex, because they are almost in a fluid state. For a long time, x-ray imaging has had an important role in diagnosing the state of objects. However, very high-energy x-rays (average energy 4 MeV) are needed to penetrate an object with a high areal density (∼200 g/cm²). Owing to the strong scattering of the x-rays and their short mean free path, the signal is so weak that it cannot be used for high-precision diagnosis.1,2

High-energy proton radiography (pRad), first developed in the 1990s in the USA,3 has obvious advantages over high-energy x-ray radiography in terms of penetration, density resolution, and multiple imaging.4–6 Multi-frame detectors have been developed for use in an 800 MeV pRad system.7,8 Hydrodynamics experiments have been performed using an 800 MeV pRad system, and changes in wave velocity and material density have been analyzed in dynamic experiments on this system.9 In Russia, both 800 MeV and 50 GeV pRad systems have been constructed, their spatial and temporal resolutions have been analyzed, and they have been used to perform a series of hydrodynamics experiments.10–15 For example, a 50 GeV pRad system was found to have a spatial resolution of 250 µm for an object with an areal density of 300 g/cm²,11 while a magnifying 800 MeV pRad system was found to have a spatial resolution of 115 µm for an object with an areal density of 17 g/cm².14 Another magnifying pRad system was designed at the Los Alamos National Laboratory to reduce dispersion.16 A 20 GeV pRad system was constructed at the Brookhaven National Laboratory and was found to have a spatial resolution of better than 200 µm at the center of the French Test Object, with a reconstructed density error of less than 1%.6

A pRad system was investigated in clinical trials in Germany, and the results showed that it could achieve a spatial resolution of better than 200 µm at a low radiation dose.17 A method for designing magnifying pRad systems using thick lenses has been proposed.18 Identical and magnifying 11 MeV pRad systems have been constructed by the China Academy of Engineering Physics. In the case of the identical pRad system, the spatial resolution of a thin object was found to be better than 100 µm. With the magnifying pRad system, a spatial resolution of 30 µm was obtained for a 10 μm-thick aluminum plate. The error in the density reconstruction of an object with an areal density of 27 mg/cm² was less than 2.3%.19–21 A high-energy pRad system has been applied to medical imaging. The spatial resolution was improved by using an inverse collimator.22,23 A high-energy electronic imaging system has been applied as a diagnostic tool for high-energy-density physics, with high spatial and temporal resolution.24,25 A 3.5–4.5 GeV pRad system has been constructed at GSI Helmholtzzentrum für Schwerionenforschung GmbH (GSI) and can achieve a spatial resolution of 30 µm and a temporal resolution of 10 ns.26 A parallel beam pRad system has been proposed to reduce image blur.27

A pRad system has quadrupole magnetic lens and collimators. It utilizes nuclear reactions and multiple Coulomb scattering (MCS) to achieve imaging. Point-to-point imaging is realized through the magnetic lens. The proton flux is controlled by the collimators, which can be used to adjust the image contrast and thus determine the density of single or mixed materials.28 Although the use of pRad for material diagnosis has been proposed,29 concrete information about the composition of a multi-material object cannot be obtained directly. This paper proposes a multi-material diagnostic method (MMDM) that can be used to find the material distribution. Furthermore, the application of this method to the diagnosis of a detonation test is presented.

In a fluid, the density of the matter generally changes with time. In a single-material diagnostic method (SMDM), the effect of density has to be eliminated, which can be realized by double imaging in pRad.

In pRad, the proton flux decays because of nuclear reactions and MCS with the angular collimators. The flux can be expressed as28I=I0⁡exp−∫μρdl1−exp−θcut22θ02,(1)where I₀ is the initial proton flux, μ is the absorption coefficient for the nuclear reactions, ρ is the density of the object, θ₀ is the root mean square of MCS, and θ_cut is the angle cut of the collimator. θ₀ can be represented asθ0=k∫ρX0dl,k=14.1MeVβpc.(2)Here, X₀ is the radiation length, β = v/c, where v is the proton velocity and c is the speed of light, and p is the proton momentum. To reduce blurring of the image, the momentum dispersion should be as small as possible.3 Therefore, k is considered to be a constant in this imaging system.

Let θ_c1 and θ_c2 denote the angle cuts of the first and second collimators. If θ_c1 ≥ 3θ₀, then the proton flux through the first angular collimator, I₁, which is the proton flux after the nuclear reactions, is given byI1=I0⁡exp−∫μρdl.(3)If θ_c2 ∼ θ₀, then the proton flux through the second angular collimator, I₂, which is the proton flux after MCS, can be obtained:I2=I11−exp−θc222θ02.(4)From Eqs. (3) and (4), we can obtain the following two relations for the absorption coefficient and the radiation length:∫μρdl=−ln(I10),I10=I1I0,(5)θc222θ02=−ln(1−I21),I21=I2I1.(6)

If the object is composed of a single material, then the characteristic parameter of this material, ω, can be obtained from Eqs. (2), (5), and (6):ω=μM=ln(I10)ln(1−I21),(7)whereM=2k2θc22X0.(8)If the second collimator is an inverse collimator, then the proton flux through it isI2′=I1⁡exp−θc222θ02,(9)and the characteristic parameter can again be obtained in the same way:ω=θc22μX02k2=ln(I10)ln(I21′),I21′=I2′I1.(10)

The instantaneous distribution of the absolute density of an object generally has to be obtained through fluid mechanics. Thus, the mass absorption coefficient or the distribution of the MCS parameter must be known. If the object is made of a single material, then the parameters required can be obtained by an SMDM. Certainly, they can also be determined from the initial state. However, if an object is composed of multiple materials, then the boundaries of the various materials may change, and the materials may even become mixed together. In this situation, it is difficult to identify the materials with an SMDM, and only the average values of the physical parameters can be calculated. Thus, SMDMs do not have wide application in fluid dynamics. We propose a method to solve some of these problems.

The relative proton flux through the scanned object can be obtained by a discretization of Eqs. (5), (6), and (8):Tfi=∑jμjρjGij,Tfi=−ln(I10),(11)Tsi=∑jMjρjGij,Tsi=1−ln(1−I21).(12)Here Tf_i and Ts_i are the relative proton fluxes due to the ith ray, μ_j is the absorption coefficient of the material in the jth voxel, ρ_j is the density for the jth voxel, M_j is related to the radiation length of the material in the jth voxel, and G_ij is the ijth element of the geometric matrix G. The characteristic parameter can be calculated by solving Eqs. (11) and (12), although these equations are usually underdetermined.

Spherically symmetric objects are often used in detonations.30,31 The geometric matrix for such an object can be divided into grids according to the number of detector pixels. Moreover, it is invertible. Therefore, we getμjρj=(G−1⋅Tf)j,(13)Mjρj=(G−1⋅Ts)j.(14)Thus, the characteristic parameter for the material in each voxel can be obtained:ωj=μjMj=(G−1⋅Tf)j(G−1⋅Ts)j.(15)This allows the different material components in the object to be identified. Moreover, the distribution of the characteristic parameter is also given by Eq. (15).

In this section, we simulate a 20 GeV pRad system with the Geant4 software toolkit.32 The quadrupole gradient of the system is 8 T/m, the quadrupole length is 2.0 m, and the drift length is 3.403 m.33 The tracks of protons through the two-stage pRad system are shown in Fig. 1(a). Protons with the same MCS angle can be focused on the Fourier plane. Thus, the proton flux can be controlled quantitatively by a collimator or by an inverse collimator. In the simulation, the angle cuts of the first and second collimators are 8 and 0.5 mrad, respectively. An imaging magnetic lens group is used before the object to measure the initial proton flux.5 The collimators are made of tungsten and are 60 cm thick. The inverse collimator is also made of tungsten and has a thickness of 70 cm. The aperture parameters of the collimators are listed in Table I. The pixel size is 1 × 1 mm². The total number of protons is 10⁸ in all simulations except those in Sec. IV B. The proton flux is along a single radial line, to make full use of the protons and reduce the running time.

μ and M were obtained for a variety of single-material objects. The characteristic parameters were calculated by Eq. (7), as listed in Table II.

To verify the proposed MMDM, a concentric spherical object (CSO) was designed, as shown in Fig. 1(b). From the inside out, the radii of the spheres are 2, 3.5, and 5 cm. The CSO was made from carbon, sodium, and beryllium. The flux after the protons pass through the CSO while undergoing nuclear reactions can be obtained at the first collimator, which has an angle cut of 8 mrad. The proton flux after the subsequent MCS can be obtained at the second collimator, which has an angle cut of 0.5 mrad.

Figure 2(a) shows the proton flux along the radial direction of the image. The red curve is the flux after the nuclear reactions, and the blue curve is the flux after MCS.

The characteristic parameter was calculated with Eq. (15), as plotted as the red and blue curves in Fig. 2(b). The only difference between the curves is that a regularization term was included for the blue curve. The large oscillations in the red curve in Fig. 2(b) were due to statistical fluctuations of the proton flux. The blue curve shows that the oscillations were significantly reduced by adding the regularization term. According to MMDM, ω was 0.194 in the interval (0, 1.5) cm, 0.101 in the interval (1.9, 3.1) cm, and 0.321 in the interval (3.5, 4.5) cm. Compared with the values in Table II, the materials in the layers of the object are as follows: carbon in (0, 1.5) cm, sodium in (1.9, 3.1) cm, and beryllium in (3.5, 4.5) cm. The green curve for the characteristic parameter in Fig. 2(b) was calculated with Eq. 7. The curve should remain unchanged because the area of object is single material (beryllium) at the outermost layer, but the result is changing because the flux distortion is induced on the edge. The phenomenon can be explained by theoretical calculation, which the result by the SMDM is good fit with the method at the outermost layer in Fig. 5(e). The value tends to 0.2 because the ω average value of the CSO is equal to 0.196 at 0. So it is coincidental that the value by the SMDM is approximately equal to the ω of carbon.

Figure 2(c) shows the ω distribution obtained by the SMDM for the CSO. It is unclear in the inner region. Figure 2(d) shows the result obtained with the MMDM. It is clear that the object has three layers. Thus, although it is difficult to identify the materials in a multi-material object with an SMDM, multi-material diagnosis can be realized by Eq. (15).

Dark-field imaging can improve the contrast resolution for thinner objects. In a dark-field pRad system, the second collimator is an inverse collimator. As before, the proton flux after the nuclear reactions was acquired at the first collimator, which has a large angle cut. The proton flux after MCS was obtained at the inverse collimator, which has an angle cut of 0.5 mrad. The curves of the proton flux for the CSO are plotted in Fig. 3(a), and ω is shown in Fig. 3(b). As can be seen from Fig. 3(b), ω is 0.203 in (0, 1.5) cm, 0.101 in (2.0, 3.1) cm, and 0.331 in (3.5, 4.5) cm. As before, the materials in the CSO can be identified from Table II.

These results demonstrate that information about the material components in an object can be identified by the proposed MMDM. There are errors due to statistical fluctuations and distortion of the proton flux. Therefore, improving the algorithm is pivotal to increasing the accuracy of the results. For example, the errors were greatly reduced by adding the regularization term. The results can also be improved with an iterative method.34

The precision of the calculated characteristic parameter is closely related to the accuracy of measuring the proton flux. The errors are caused by nonideal factors in the pRad system, as well as statistical fluctuations. One common practical problem is that the central axes of the components do not coincide. The collimators are important components in an pRad system, since they control the proton flux. Thus, the errors will be increased if a collimator deviates from the central axis. In particular, the proton flux after MCS depends completely on the second collimator. The size of the aperture at the back end of the collimator is proportional to the angle cut. Protons scattered by more than this angle cannot pass through the collimator. The angular distribution of protons after MCS depends on the object and the proton energy.28 Therefore, the error in the characteristic parameter caused by collimator deviation depends on the specific object, the proton energy, and the pRad system.

We analyzed the impact of collimator deviation on the characteristic parameter through simulation. The object and collimators were the same as the settings in Sec. IV A. We used the following offsets of the second collimator from the central axis: 0 (no deviation), 100, 300, 500, and 800 µm. The total number of protons was 10⁷.

Figure 4(a) shows the proton flux along the radial direction for the different offsets. The resulting statistical fluctuations were worse than those in Fig. 2(a) because the total number of protons was reduced by 9 × 10⁷. However, this test was very useful, since the fluctuations were basically consistent at each position in each simulation. It shows that interference caused by statistical fluctuations can be significantly reduced. The flux error plotted in Fig. 4(b) was calculated asIdeviated−Inot deviatedInot deviated×100%.The error increased with increasing deviation.

The radial values of ω calculated with the optimized equation are plotted in Fig. 4(c). The values at positions less than 0.1 cm with no offset are obviously different from those in Fig. 2(b). This is because the statistical fluctuations were larger, owing to the fact that although the conditions are the same, the total number of protons is lower. This result also demonstrates that the statistical fluctuations have a significant impact on the accuracy of the characteristic parameter.

Figure 4(d) shows the ω error:ωdeviated−ωnot deviatedωnot deviated×100%.The ω error basically increased with increasing offset. It was affected by statistical fluctuation at a few positions. The green curve is for an offset of 100 µm. It has a maximum error of 1.5% at 0.3 cm. The minimum error was 0.0% at 2.7 cm, and the mean error was 0.4%. For a deviation of 300 µm, the mean error was 0.8%. It was 1.4% for a deviation of 500 µm. Thus, the mean ω error for the object, collimator, and pRad system can be less than 1% if the deviation is less than 300 µm. These results demonstrate that the nonideal factors of the system have an important influence on the accuracy of the characteristic parameter. The specific impact needs to be analyzed quantitatively for each situation.

Detonation is an important topic in fluid dynamics. In a detonation, the components mix together, and the boundaries between them usually change. However, high-energy x-ray radiography cannot determine the material compositions of the components. Only the product of the mass absorption coefficient and the density can be obtained by x-ray imaging. Since the density of each material is not constant in a fluid, the product could be for any material. In this subsection, we determine the boundary and material information of each component of a dynamic object at a certain moment using the MMDM.

A CSO was designed. Initially, the materials in the layers from inside to outside are boron, sodium, and beryllium. At some point, the two components in the inner layer mix together, so that the layers change to boron, a boron and sodium mixture, sodium, and beryllium. This process is illustrated in Fig. 5(a). At this moment, the radius of the boron region is 1.5 cm, and its density is 2.370 g/cm³. The radius of the mixed region is 2.3 cm. The density of the boron there is 1.247 g/cm³, and the density of the sodium is 0.460 g/cm³. The radius of the sodium region is 3.5 cm, and its density is 0.971 g/cm³. The radius of the beryllium region is 5 cm, and its density is 1.848 g/cm³ (Table III).

First, we calculated the proton flux of the object using Eqs. (10) and (11), as shown in Fig. 5(b). Figure 5(c) shows the distribution of the characteristic parameter obtained with the SMDM. The value is constant in the outermost layer, but gradually changes in the other regions, and so it is difficult to determine the boundaries of the components. Figure 5(d) shows the ω distribution obtained with the MMDM. It is quite clear that the object has four components.

The ω distribution in the radial direction is plotted in Fig. 5(e). The blue curve is for the MMDM. It can be seen that the characteristic parameter has constant values of 0.245 in the region (0, 1.5) cm, 0.185 in (1.6, 2.3) cm, 0.099 in (2.4, 3.5) cm, and 0.340 in (3.6, 4.9) cm. Using information about the initial object and Table II, it can be deduced that the regions (0, 1.5), (2.4, 3.5), and (3.6, 4.9) cm are boron, sodium, and beryllium, respectively. As the value of the characteristic parameter in the region (1.6, 2.3) cm is not that of any of the original components, but lies between those of boron and sodium, it can be deduced that this region is a mixture of boron and sodium. The green curve was obtained with the SMDM. It is a good fit with the blue curve in the outermost layer, but not in the other layers, and so we cannot distinguish the distribution of materials in the CSO.

The distributions of μ and M can be obtained from the MMDM results. The density distribution of the object can be reconstructed with Eqs. (10) and (11), as shown in Fig. 5(f). Table III shows that the reconstructed densities are in good agreement with the designed densities, because the theoretical proton flux is accurate.

Second, we obtain the proton flux of the object via a numerical simulation with two collimators, as shown in Fig. 5(g). The ω distribution from the SMDM is plotted in Fig. 5(h) and that from the MMDM in Fig. 5(i). We can draw the same conclusions as we did for the theoretical results.

The characteristic parameter in the radial direction is shown in Fig. 5(j). The red curve is for the MMDM, and the green curve is for the SMDM. Compared with the theoretical results, the precision of ω is seriously affected by statistical fluctuations and distortion of the flux at the boundaries. Using the ω distribution obtained with the MMDM, we can deduce that the region (0, 1.3) cm is boron, (1.55, 2.05) cm is a mixture of boron and sodium, (2.25, 3.2) cm is sodium, and (3.5, 4.7) cm is beryllium.

The reconstructed density distribution is shown in Fig. 5(k). Owing to the statistical fluctuations and distortion of the flux at the boundaries, there are differences between the reconstructed and designed densities (Table III). The results do indicate that information about the boundaries and material composition of a dynamic object can be obtained with the MMDM. Moreover, the densities of both components in a mixture can also be reconstructed.

The characteristic parameter ω is completely independent of density. ω provides information about the material composition, regardless of how the density is changing. In this paper, we have proposed a method to obtain the ω distribution of a multi-material object. We have considered its application to a detonation test based on numerical simulations of a CSO. A 20 GeV pRad system has been modeled in Geant4. The proton flux has been obtained with both a collimator–collimator setup and a collimator–inverse collimator setup, both with angle cuts of 8 and 0.5 mrad. The results show that the proposed MMDM can identify information about the materials in a multi-material object. In fluid dynamics, it is pivotal to obtain the material and boundary information about each component of a multi-material object. Our results show that the proposed method can provide this information. In addition, the MMDM can be used to obtain the density distribution of a dynamic object, even if the components are mixing together. Therefore, the proposed MMDM may have significant potential for fluid dynamicaldiagnosis.

Acknowledgment. This work was supported by the National Natural Science Foundation of China (Grant No. 11675021) and the Youth Fund Project of Zhejiang Lab (Grant No. K2023MG0AA08).