A novel 4D resolution imaging method for low and medium atomic number objects at the centimeter scale by coincidence detection technique of cosmic-ray muon and its secondary particles

A novel 4D resolution imaging method for low and medium atomic number objects at the centimeter scale by coincidence detection technique of cosmic-ray muon and its secondary particles 1 corresp first-author wangxd@usc.edu.cn wangxd@usc.edu.cn wangxd@usc.edu.cn https://orcid.org/0000-0003-4006-5282 1 1 1 1 1 1 1 School of Nuclear Science and Technology, University of South China, Hengyang 421001, China School of Nuclear Science and Technology, University of South China, Hengyang 421001, China The muon radiography imaging technique for high-atomic-number objects ( Z ) and large-volume objects via muon transmission imaging and muon multiple scattering imaging remains a popular topic in the field of radiation detection imaging. However, few imaging studies have been reported on low and medium Z objects at the centimeter scale. This paper presents an imaging system that consists of three layers of a position-sensitive detector and four plastic scintillation detectors. It acquires data by coincidence detection technique of cosmic-ray muon and its secondary particles. A 3D imaging algorithm based on the density of the coinciding muon trajectory was developed, and 4D imaging that takes the atomic number dimension into account by considering the secondary particle ratio information was achieved. The resultant reconstructed 3D images could distinguish between a series of cubes with 5-mm side lengths and 2-mm intervals. If the imaging time is more than 20 days, this method can distinguish intervals with a width of 1 mm. The 4D images can specify target objects with low, medium, and high Z values. Image reconstruction Monte Carlo simulation Non-destructive detection Introduction 1 In 1936, a cosmic-ray muon was discovered in an electromagnetic field [ 1 1 ] . A muon is a particle produced by high-energy cosmic rays interacting with atoms in the Earth’s atmosphere through cascade clustering, with a mass of 105.7 MeV/C 2 , which is 207 times the mass of an electron, and a lifetime of approximately 2.2 microseconds. It reaches the Earth’s surface with an average energy of approximately 3–4 GeV (maximum energy can reach TeV) and a flux scale of 1 cm -2 min -1 . The angular distribution of the muon intensity is approximately proportional to cos 2 θ , where θ is the zenith angle [ 2 2 ] . Cosmic-ray muons lose their kinetic energy by Coulomb scattering, which is described by the Bethe-Block formula [ 3 3 ] . Apart from the cosmic-ray muon source, the determination of an artificial muon source is a popular research topic [ 4 4 - 7 7 4-7 ] . <math display="block" id="M1"><mrow><mrow><mo>〈</mo><mrow><mfrac><mrow><mtext>d</mtext><mi>E</mi></mrow><mrow><mtext>d</mtext><mi>x</mi></mrow></mfrac></mrow><mo>〉</mo></mrow><mo>=</mo><mi>K</mi><msup><mi>z</mi><mn>2</mn></msup><mfrac><mi>Z</mi><mi>A</mi></mfrac><mfrac><mn>1</mn><mi>β</mi></mfrac><mrow><mo>[</mo><mrow><mfrac><mn>1</mn><mn>2</mn></mfrac><mi>ln</mi><mfrac><mrow><mn>2</mn><msub><mi>m</mi><mtext>e</mtext></msub><msup><mi>c</mi><mn>2</mn></msup><msup><mi>β</mi><mn>2</mn></msup><mi>γ</mi><msubsup><mi>W</mi><mrow><mi>max</mi></mrow><mn>2</mn></msubsup></mrow><mrow><msup><mi>I</mi><mn>2</mn></msup></mrow></mfrac><mo>−</mo><msup><mi>β</mi><mn>2</mn></msup><mfrac><mrow><mi>δ</mi><mo stretchy="false">(</mo><mi>β</mi><mo stretchy="false">)</mo><mi>γ</mi></mrow><mn>2</mn></mfrac></mrow><mo>]</mo></mrow></mrow></math> Cosmic-ray muons have the unique advantage of being highly penetrating and non-destructive. Therefore, they can be used in imaging techniques. Current muon imaging technologies can be divided into muon multiple scattering imaging and muon transmission imaging. Muon multiple scattering imaging was first applied to detect high-atomic-number ( Z ) objects in a short period of time by analyzing the scattering angle of muons together with the point of closest approach (POCA) algorithm in 2003 [ 8 8 ] . Subsequently, this technology is primarily applied for imaging high- Z objects. A large-area muon detector system named CRIPT was constructed, which is based on scintillators and wavelength-shifting fibers to detect potential uranium or plutonium material in a container [ 9 9 , 10 10 ] . After the Fukushima nuclear accident, nuclear material debris inside the reactor were imaged using the muon scattering imaging technique [ 11 11 ] . It can also be used to identify different materials [ 12 12 - 14 14 12-14 ] . For muon transmission imaging, this technology reconstructs the image of the target objects by analyzing the change of muon flux after muons penetrate objects and is generally aimed at a large volume and natural terrain. In 1970, hidden chambers in pyramids were discovered by applying the muon transmission imaging technique [ 15 15 ] . Various research groups have investigated the application of similar techniques to volcano monitoring [ 16 16 - 21 21 16-21 ] , bulk ship imaging [ 22 22 ] , and lunar shadow observations [ 23 23 ] . In 2017, the study on the imaging of a large void in Khufu’s pyramids using muon transmission imaging was published [ 24 24 ] . Furthermore, significant progress was made in China in 2020, where a research team completed the muon transmission imaging of the Changshu Underground Trench [ 25 25 ] . In addition, many types of position-sensitive detectors to detect muons have been designed to adapt to different target objects. For the imaging of targets with large volumes, such as volcanoes and pyramids, many research teams have applied scintillation detectors as position-sensitive detectors [ 26 26 - 28 28 26-28 ] . These detectors have very large sensitive volumes (more than 1000 mm <math id="M2"><mo>×</mo></math> 1000 mm), but their position resolution is relatively poor (less than 10 mm). For imaging materials with a high Z , a series of gas detectors with a high position resolution (better than 0.5 mm) has been constructed [ 29 29 - 31 31 29-31 ] , and Micromegas detection with a position resolution of 0.075 mm was introduced [ 32 32 ] . Both muon transmission imaging and muon multiple scattering imaging essentially reconstruct images based on the change in muon trajectory or flux before and after the muon penetrates the target objects. They are more adaptable for reconstructing high- Z objects and imaging large-volume objects. The secondary particles produced by the muon also carry information about the target object. If the information carried by these secondary particles can be fully utilized, more characteristics of the target object can be observed in a single experiment. Based on secondary particle information, the Muon Camera (MUCA) was designed and used for experiments on imaging bovine bones [ 33 33 - 35 35 33-35 ] , compensating for the insufficient imaging of low and medium  Z  objects. However, the information implied by secondary particles should be further explored for muon imaging. A new imaging system based on muons and its secondary particle detection technique is proposed and designed in this work. The system consists of three layers of a position-sensitive detector and four fast time response plastic scintillation detectors with a total volume of 0.125 m 3 . It fully combines the advantages of a high position resolution and fast response time. Incident muon trajectories have been obtained using position-sensitive detectors, and secondary particles have been measured using plastic scintillators. In terms of the imaging principle, an algorithm was developed for 3D imaging based on the coinciding muon trajectory density. The resolution of the atomic number dimension is accomplished by secondary particle information, that is the reconstruction of 4D images. Secondary particle production and performance 1 The physical processes, energy spectra, and polar angle distributions of the secondary particles are studied in this section. It should be noted that secondary particles include secondary gamma rays, electrons, and a few positrons. As the number of positrons is several orders of magnitude lower than the those of the other particles, only secondary gamma rays and electrons were considered here. Geant4 simulation 2 Geometry of the detector system 3 Simulations were performed using GEANT4.10.6p-03. The generic model was visualized using the Geant4 OpenGL interface, as described in Figure 1 Figure 1 . The detector system consisted of two parts. The first part was the muon trajectory detector module, which was made up of three layers of position-sensitive detectors and used to obtain the precise position of muon hits and reconstruct muon trajectories. Each detector was 400 mm <math id="M3"><mo>×</mo></math> 400 mm <math id="M4"><mo>×</mo></math> 10 mm in size, with a distance of 50 mm between the bottommost top detector and the topmost bottom detector. This part was located at the top of the detector system. They consisted of micro-filter gaseous detectors. The second part was defined as the secondary particle detector module, which was used to detect secondary particles. It consists of four EJ200 scintillation detectors and has a box-like shape surrounding the target object. This module was designed to be very large to detect as many secondary particles as possible. The dimensions of each scintillation detector were 500 mm <math id="M5"><mo>×</mo></math> 500 mm <math id="M6"><mo>×</mo></math> 50 mm. Compared with the muon multiple scattering imaging system, the new imaging system only added four scintillation detectors; thus, the two can be readily adapted to each other. The white cube at the center was a target object with dimensions of 20 mm <math id="M7"><mo>×</mo></math> 20 mm <math id="M8"><mo>×</mo></math> 20 mm. When an incident muon penetrated the cube, there was a certain probability that secondary particles would be produced and recorded by the secondary particle detector module. A series of Geant4 simulations were developed to study secondary particle generation processes and to carry information about the target objects. Primary particles 3 The G4VUserPrimaryGeneratorAction class in Geant4 defines the primary particles. The particle name was defined as the mu-, and the energy spectrum obeyed the natural distribution; more detailed information can be found in M. Tanabashi, P.D. Grp, K. Hagiwara et al. [ 3 3 ] (The average energy is 3–4 GeV, and the maximum energy can reach the TeV level). In order to simulate the zenith angle of the cosmic-ray muon, the momentum direction was selected from the cos 2 θ distribution. Therefore, the larger the polar angle and the smaller the amount of incident muon, the more the amount of the initial muon that is inserted in the vertical angle will be. Each initial position of the muon is chosen as follows: The z -coordinate was fixed above the muon trajectory detector module. The x -coordinate and y-coordinate were selected from the uniform range, from x min to x max and y min to y max, respectively. Physics lists 3 Physical processes can be defined by the G4VModularPhysics class in Geant4. There were four modules of physics in the simulation: G4DecayPhysics, G4RadioactiveDecayPhysics, G4EmstandardPhysics, and G4OpticalPhysics. EventAction 3 The EventAction class is defined based on G4UserEventAction. At the beginning of each event, the muon trajectory detector module and secondary particle detector module were assigned a value of 0. When the secondary particles passed through the scintillation detectors, the value of the secondary particle detector module was assigned as 1, and the information of the secondary particles was recorded. Secondary particles produced by a muon 2 Ionization, bremsstrahlung, and pair production are the three principal types of interactions between muons and the target objects. When the target object’s material is Fe, the energy loss of the muon by different interactions with limited physical processes is shown in Fig. 2 Fig. 2 . It can be seen that the ionization energy loss varies little with the change in its energy, and this is the primary means through which muons lose their energy below 570 GeV. In addition, the energy loss caused by the bremsstrahlung and pair effect increases rapidly with the increase in muon energy and finally takes the dominant position after the muon energy goes beyond 570 GeV. Thus, for cosmic ray muons with an average energy of approximately 3–4 GeV, the secondary electrons are primary ionized δ electrons and secondary ionized electrons. Secondary gamma is produced by the bremsstrahlung of the secondary electrons. Secondary particle energy spectrum 2 Figure 3(a) and 3(b) Figure 3(a) and 3(b) shows the secondary electron and gamma energy spectra derived from the Geant4 simulation results, respectively. In terms of the secondary electrons, the energy spectrum did not vary significantly for different materials. The counts first increase and then decrease with Z , owing to ionization and self-absorption effects that vary with Z . Therefore, a target objects with a medium Z has the highest secondary electron count. As for the secondary gamma, because the δ electrons produced by muon ionization can continue to ionize with the target object, and then bremsstrahlung, the peak of the secondary gamma energy spectrum is primarily concentrated at below 1 MeV. Furthermore, with an increase in Z , the peak shifts to the right and the counts increase owing to the growth of ionization energy loss and the weak self-absorption of the gamma rays. It can be seen that the counts of secondary particles vary considerably in the energy range of 0.1 MeV to 100 MeV. Therefore, the counts of secondary electrons and secondary gamma in this range can potentially distinguish Z of the target object.  The polar angle distribution of secondary particles 2 Figure 3(c) and 3(d) Figure 3(c) and 3(d) show the polar angle distribution of the secondary particles. To obtain the angle information of all secondary particles, the angle leaving the target object was recorded in place of the information recorded by the scintillators. The horizontal coordinate represents the angle between the secondary particle and the muon trajectory. It can be seen that both the secondary electrons and gamma rays are most likely to exit along the muon trajectory. For Pb and U, the probability of the secondary electrons being emitted in the direction opposite to the incident muon is higher than that of other materials. This is because the secondary electrons are more likely to be deflected at a large angle for target objects with a higher Z . In addition, Z of the target objects has a limited effect on the polar angle distribution. Principle 1 Coinciding muon trajectory 2 In the Geant4 simulation, the incident muon trajectory is only recorded as a coincident muon trajectory if both the muon trajectory detector module and the secondary particle detector module respond in a short time window. Because the speed of muons and their secondary particles is very close to the speed of light, this short time is approximately 1 ns. Furthermore, the EJ200 scintillator for secondary particle detection has fast response time (rise time 0.9 ns, decay time 2.1 ns) and short pulse width (FWHM=2.5 ns) to distinguish secondary particles from background noise. Because secondary particles are generated by the process in which muons pass through the material and then ionize, the coinciding muon trajectory is a line that penetrates the target material in an ideal environment. Considering that an actual experimental environment exists in air, the coinciding muon trajectory has a high probability of passing through the target object. Because the muon trajectory changes very little when penetrating the low and medium Z objects, it is assumed that the incident muon trajectory does not change; that is, the trajectory is a straight line. If sufficient coincident muon trajectories have been obtained, the target object can be reconstructed. Imaging algorithm 2 The section where a large number of coinciding muon trajectories cross each other reflects the area where the target object is located. A method to confirm this is shown in Fig. 4 Fig. 4 . The first step is to set a series of mathematical planes and then calculate each coordinate of the intersection points of all coinciding muon trajectories with all given mathematical planes, as shown in Fig. 4(a) Fig. 4(a) . The second step is to select the highest density area of the intersection points, as shown in Fig. 4(b) Fig. 4(b) . The intersection point is represented by a dot. The hollow white dot indicates that the intersection point is outside the object, and the solid green dot indicates that it is inside the object or it indicates a good reconstruction event. The third step is to convert the selected point into the assignment of each voxel, as shown in Fig. 4(c) Fig. 4(c) . As a result, by showing the assignment of each voxel, the 2D and 3D images of the target object can be reconstructed. When each given mathematical plane is perpendicular to the Z -axis, the method of calculating the intersection point coordinates is equivalent to calculating the intersection of a three-dimensional line with a plane perpendicular to the Z -axis. The distribution of the intersection points along the Z -axis is different, and the highest density of intersection points reflects the Z -value of the target object’s height. The production of secondary particles was dominant during this interval. When the given mathematical plane is higher or lower than the height of the target object, the distribution of the intersection points shows a radial and mesh pattern. Thus, if we select a high-density intersection area by a certain threshold, the height of the target material will be determined. This algorithm can be modified one step further. The given mathematical planes are not perpendicular to the Z axis, as shown in Fig. 4(d) Fig. 4(d) . Similarly, the target objects distribution along the direction of the normal vector to the mathematical plane can be obtained. Furthermore, it is possible to obtain the distribution of the object along different normal vectors by continuously changing the given mathematical plane, which can increase the accuracy of the tomography. However, it is more complex to calculate the intersection point coordinates when the given mathematical plane is not perpendicular to the Z -axis. A less complicated calculation method is displayed as follows: First, we establish a new coordinate system X <math id="M12"><mtext>’</mtext></math> Y <math id="M13"><mtext>’</mtext></math> Z <math id="M14"><mtext>’</mtext></math> , whose Z <math id="M15"><mtext>’</mtext></math> -axis is perpendicular to the given mathematical plane, as shown in Fig. 4(e) Fig. 4(e) . The basis vectors ( i , j , k ) of the original coordinate system XYZ can be converted to the basis vectors ( i <math id="M16"><mtext>’</mtext></math> , j <math id="M17"><mtext>’</mtext></math> , k <math id="M18"><mtext>’</mtext></math> ) of the new coordinate system X <math id="M19"><mtext>’</mtext></math> Y <math id="M20"><mtext>’</mtext></math> Z <math id="M21"><mtext>’</mtext></math> using the following equation: <math display="block" id="M22"><mrow><mrow><mo>(</mo><mrow><msup><mi>i</mi><mo>'</mo></msup><mo>,</mo><msup><mi>j</mi><mo>'</mo></msup><mo>,</mo><msup><mi>k</mi><mo>'</mo></msup></mrow><mo>)</mo></mrow><mo>=</mo><mrow><mo>(</mo><mrow><mi>i</mi><mo>,</mo><mi>j</mi><mo>,</mo><mi>k</mi></mrow><mo>)</mo></mrow><mrow><mo>(</mo><mrow><mtable><mtr><mtd><mrow><mi>cos</mi><msub><mi>α</mi><mn>1</mn></msub></mrow></mtd><mtd><mrow><mi>cos</mi><msub><mi>β</mi><mn>1</mn></msub></mrow></mtd><mtd><mrow><mi>cos</mi><msub><mi>γ</mi><mn>1</mn></msub></mrow></mtd></mtr><mtr><mtd><mrow><mi>cos</mi><msub><mi>α</mi><mn>2</mn></msub></mrow></mtd><mtd><mrow><mi>cos</mi><msub><mi>β</mi><mn>2</mn></msub></mrow></mtd><mtd><mrow><mi>cos</mi><msub><mi>γ</mi><mn>2</mn></msub></mrow></mtd></mtr><mtr><mtd><mrow><mi>cos</mi><msub><mi>α</mi><mn>3</mn></msub></mrow></mtd><mtd><mrow><mi>cos</mi><msub><mi>β</mi><mn>3</mn></msub></mrow></mtd><mtd><mrow><mi>cos</mi><msub><mi>γ</mi><mn>3</mn></msub></mrow></mtd></mtr></mtable></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> where α 1 , α 2 , and α 3 denote the angle between the X -axis and X <math id="M23"><mtext>’</mtext></math> -axis, Y <math id="M24"><mtext>’</mtext></math> -axis as well as Z <math id="M25"><mtext>’</mtext></math> -axis, respectively. Similarly, β 1 , β 2 , β 3 , γ 1 , γ 2 , and γ 3 represent the angles between the Y -axis, Z -axis, and the axes of the new coordinate system. Thus, in the new coordinate system X <math id="M26"><mtext>’</mtext></math> Y <math id="M27"><mtext>’</mtext></math> Z <math id="M28"><mtext>’</mtext></math> , the coordinates of the intersection point system can be calculated by imitating the method in which the given mathematical plane is perpendicular to the Z -axis. Second, as shown in Fig. 4(f) Fig. 4(f) , the coordinates of the intersection points in the new coordinate system X <math id="M29"><mtext>’</mtext></math> Y <math id="M30"><mtext>’</mtext></math> Z <math id="M31"><mtext>’</mtext></math> are converted to the coordinates in the original coordinate system XYZ using the following equation: <math display="block" id="M32"><mrow><mrow><mo>(</mo><mrow><msub><mi>x</mi><mn>1</mn></msub><mo>,</mo><msub><mi>y</mi><mn>1</mn></msub><mo>,</mo><msub><mi>z</mi><mn>1</mn></msub></mrow><mo>)</mo></mrow><mo>=</mo><mrow><mo>(</mo><mrow><msub><mi>x</mi><mn>1</mn></msub><msup><mrow/><mo>'</mo></msup><mo>,</mo><msub><mi>y</mi><mn>1</mn></msub><msup><mrow/><mo>'</mo></msup><mo>,</mo><msub><mi>z</mi><mn>1</mn></msub><msup><mrow/><mo>'</mo></msup></mrow><mo>)</mo></mrow><mrow><mo>(</mo><mrow><mtable><mtr><mtd><mrow><mi>cos</mi><msub><mi>α</mi><mn>1</mn></msub></mrow></mtd><mtd><mrow><mi>cos</mi><msub><mi>α</mi><mn>2</mn></msub></mrow></mtd><mtd><mrow><mi>cos</mi><msub><mi>α</mi><mn>3</mn></msub></mrow></mtd></mtr><mtr><mtd><mrow><mi>cos</mi><msub><mi>β</mi><mn>1</mn></msub></mrow></mtd><mtd><mrow><mi>cos</mi><msub><mi>β</mi><mn>2</mn></msub></mrow></mtd><mtd><mrow><mi>cos</mi><msub><mi>β</mi><mn>3</mn></msub></mrow></mtd></mtr><mtr><mtd><mrow><mi>cos</mi><msub><mi>γ</mi><mn>1</mn></msub></mrow></mtd><mtd><mrow><mi>cos</mi><msub><mi>γ</mi><mn>2</mn></msub></mrow></mtd><mtd><mrow><mi>cos</mi><msub><mi>γ</mi><mn>3</mn></msub></mrow></mtd></mtr></mtable></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> where x 1 , y 1 , and z 1 are the intersection coordinates when the given mathematical plane is not perpendicular to the Z -axis. Reconstruction of 3D image 1 The result of 3D imaging 2 As ionization loss and self-absorption effects are influenced by the Z of the target object, the imaging effects of various materials will be different. Figure 5(a) and 5(b) Figure 5(a) and 5(b) show the 3D image and top view image of five cubes with a volume of 20 mm×20 mm×20 mm. The materials of the cubes were H 2 O, Al, Fe, Pb, and U, respectively. It is important to note that because the amount of data of U is much more considerable than that of H 2 O, direct imaging will not show water. Hence, simulations were performed for each of the five materials under the same conditions and then converged to compare the imaging effect. As for the images of Pb and U, due to the strong self-absorption effect of the high- Z material, secondary particles generated inside the cube are difficult to penetrate, and there are large voids in the Pb and U images at the center of the cubes. Therefore, while imaging a high- Z target object, there is a chance that the central part’s information of the object will be lost. In addition, when the muon passes through the high- Z target object, it cannot be assumed that the muon trajectory is a straight line. Therefore, when applying this algorithm to objects with a higher Z , the error will be larger than that of other materials. As for Al and H 2 O, since the ionization loss increases with Z , meaning that low- Z materials have a smaller amount of data of the coincident muon trajectories, there are also some voids in the images of Al and H 2 O. However, because this problem is caused by insufficient data, it can be solved by extending the imaging time. Moreover, Fe has the best imaging result compared to the above materials. In summary, this method can be used for imaging objects of medium and low Z , with medium Z imaging being the best, but not perfectly applicable to high Z . 3D imaging accuracy 2 To better study the imaging accuracy of this method, Fe was used for the eight cubes because it provides the best imaging effect. Figure 5(c) and 5(d) Figure 5(c) and 5(d) show a 3D image and a top-view image of eight iron cubes with a height of 20 mm. For one column, the intervals of the cubes gradually decreased to 4 mm, 2 mm, and 1 mm. For the other column, the volume of cubes progressively decreased in the order of 20 mm, 15 mm, 10 mm, and 5 mm. All eight cubes were imaged simultaneously. Because the best position resolution achievable in the experiment is approximately 0.1 mm, in order to simulate the actual experimental condition, the position information of the muon trajectory detector module randomly added 0.1 mm position resolution with a Gaussian distribution. There were approximately 9×10 6 muons generated in the muon trajectory detector module, which was equal to 4 days muon flux. As a result, the size and position of all the other restricted cubes match the modelling of Geant4. A cube with a side length of 5 mm and an interval of 2 mm were clearly visible. However, a 1-mm interval of cubes in the first line cannot be observed. In addition, because the volume of the cube with a side length of 5 mm is too small, after determining the height of the target object using the threshold value, the height of the reconstructed image is approximately 4 mm shorter than that of the other cubes. As a result, after setting the muon flux to a 4 days imaging time and adding a 0.1-mm position resolution to the imaging system, this method can accurately reconstruct the size and position of the target object, and its accuracy reaches a level that targets objects with a side length of 5 mm, and an interval of 2 mm can be easily distinguished. Influence of position resolution to 3D imaging 2 Position resolution is a very important property for detectors with a good sensitivity to position. This means that the minimum interval between two incident particles can be distinguished by a position-sensitive detector. Because of this, position resolution directly influences the quality of the reconstructed image. Considering that the position resolution is no better than 0.1 mm on most occasions, the imaging quality of poorer position resolution has been studied. In the Geant4 simulation, maintaining the muon flux to a 4 days level, the reconstructed images with position resolutions of 0 mm, 0.2 mm and 0.4 mm are shown in Figs. 5(e), 5(f), and 5(g) Figs. 5(e), 5(f), and 5(g) . As the position resolution increases, the edge of the reconstructed image becomes increasingly dim, and its volume increases slightly. The perfect position resolution was 0 µm. Based on this condition, the 2-mm interval of cubes in the first line was clear, but the 1-mm interval remained invisible. When the position resolution grew to 0.2 mm, the 2-mm interval was mostly filled, making them difficult to distinguish. Finally, at 0.4 mm, the widest 4-mm interval can hardly be distinguished. In addition, among the four cubes of the second line, the last cube with a side length of 5 mm was reconstructed to an image that is close to a circle, which is not equal to the modeling. However, the same distortion did not occur in the three bigger cubes in front of it. A conclusion can be drawn by referring to the imaging effect with a position resolution of 0.1 mm shown in Fig. 5(d) Fig. 5(d) . That is, when the position resolution is better than 0.1 mm, the 2-mm interval can be distinguished, but as it falls to 0.2 mm, the minimum distance between two distinguishable objects has to grow to 4 mm. If the position resolution is only 0.4 mm, the 4-mm intervals and 5-mm side length will all be indistinguishable. However, despite the poor position resolution, cubes with side lengths of 20 mm, 15 mm, and 10 mm can still be correctly reconstructed. Influence of imaging time to 3D imaging 2 The imaging time directly influences the number of acquired data. After adding a 0.1-mm position resolution to the muon trajectory detector, the imaging results with imaging times of 20 days, 1 day, and 0.25 days were simulated, as shown in Figs. 5(h), 5(i), and 5(j) Figs. 5(h), 5(i), and 5(j) . With 20 days imaging time, the minimum distinguishable interval is 1 mm, because there are approximately 5×10 7 muons injected in the detector. When the imaging time is 4 days or 1 day, the 2- and 4-mm interval can be distinguished. However, if the imaging time to 0.25 days, a significant shower distortion will occur on the reconstructed image and the 4-mm interval together with the cubes in the first line will be difficult to discern. Overall, if the imaging time is 20 days, 4 days, and 1 day, respectively, this method will accurately distinguish cubes with a side length of 5 mm and respective intervals of 1 mm, 2 mm, and 4 mm. If the imaging time is no longer than 0.25 days, reconstructing objects at the centimeter scale in three dimensions would be challenging. Reconstruction of 4D image 1 The ratio of secondary electron to secondary gamma 2 According to Fig. 3(a) and 3(b) Fig. 3(a) and 3(b) , the counts of the secondary electrons and gamma vary with the Z of the target object. Thus, it is possible to study the likeliness of secondary particle generation to distinguish the Z of the target objects in the 0.1 MeV to 100 MeV range, as shown in Fig. 6(a) Fig. 6(a) . For Pb, 100 incident muons produce approximately 17 secondary electrons and 74 secondary gamma rays. It can be seen that the ratio of these two secondary particles varies with the material. Thus, the ratio of secondary electron counts to secondary gamma counts with various values of Z was calculated, as shown in Fig. 6(b) Fig. 6(b) . As Z increased, the ratio consistently decreased. Therefore, the Z of the target object can be determined by the section in which the ratio lies. Result of 4D imaging 2 4D imaging adds atomic number dimension resolution to 3D imaging. The 4D images were reconstructed as follows: First, secondary electrons and secondary gamma are used to reconstruct the coinciding muon trajectories and reconstruct the 3D images, respectively. In this way, two three- -dimensional matrices can be obtained. Then, we divide the corresponding voxels in the two three-dimensional matrices and obtain a new three-dimensional matrix. The assignment of voxels in this new matrix is the ratio of secondary electrons to secondary gamma. If these voxels are converted into colors, a 4D image is achieved. The 4D imaging results of an imaging time of 20 days and a 0.1-mm position resolution are shown in Figs. 7(a) and 7(b) Figs. 7(a) and 7(b) . The target objects were five cubes with dimensions of 20 mm×20 mm×20 mm, whose materials were H 2 O, Al, Fe, Pb, and U, respectively. All five cubes were imaged at the same time. From the top view, it is apparent that the five materials can be distinguished by the ratio of secondary electrons to secondary gamma. However, loops with different colors surround the cubes because the secondary electrons generated on the surface of the target objects are not easily absorbed, leading to a higher voxel ratio compared to that on the interior. It should be noted that the ratio also varies with volume. Therefore, a database that includes different volume ratios and different material ratios can better specify the target objects. Influence of imaging time to 4D imaging 2 The 4D imaging results for different imaging times are shown in Fig. 7(c), (d) and (e) Fig. 7(c), (d) and (e) . When the imaging time is 5 days, the imaging system can distinguish objects with low, medium, and high Z . If the imaging time decreases to 12 h, the target objects of medium and high Z will be difficult to discern, but the target objects with low Z can be clearly visualized. As the imaging time decreases to only 3 h, the reconstructed 4D image fails to distinguish the shape of target objects, but is still available to specify low- Z target objects from those with medium Z . Therefore, target objects with low and medium Z can be easily distinguished, but those with medium and high Z are less distinguishable. Furthermore, it is widely known that muon multiple scattering imaging is more applicable to target objects with medium and high Z . Considering that these two imaging systems can adapt to each other, a much better imaging result will be achieved by combining them. Conclusion 1 A muon imaging method based on the coincident muon trajectory was studied, which includes the physical processes of secondary particle generation, 3D imaging, and 4D imaging. Secondary particles are primarily produced by muon ionization and δ electron bremsstrahlung, including electrons and gamma rays. The secondary electron energy spectrum does not vary much with Z , while the peak of the secondary gamma energy spectrum gradually moves to the higher energy part as Z increases. For the imaging result, for the condition of 0.1-mm position resolution on the muon trajectory detector module, 3D images with an acceptable level of accuracy of 5-mm side length and 2-mm intervals can be obtained. If the imaging time is longer than 20 days, this method can characterize a 1-mm interval. In 4D imaging, low-, medium-, and high- Z target objects can be distinguished. However, when the imaging time is less than 0.5 days, it can only distinguish target objects with low Z between those with medium or high Z . Moreover, a database of different volumes and different materials would perform better. References S.H. Neddermeyer , C.D. Anderson , Cosmic-ray particles of intermediate mass . Phys. Rev. 54 , 88 ( 1938 ). doi: 10.1103/PhysRev.54.88.2 http://doi.org/10.1103/PhysRev.54.88.2 Cosmic-ray particles of intermediate mass J.W. Lin , Y.F. Chen , R.J. Sheu et al. , Measurement of angular distribution of cosmic-ray muon fluence rate . Nucl. Instru. and Meth. A 619 , 24 - 27 ( 2010 ). doi: 10.1016/j.nima.2009.12.017 http://doi.org/10.1016/j.nima.2009.12.017 Measurement of angular distribution of cosmic-ray muon fluence rate M. Tanabashi , P.D. Grp , K. Hagiwara et al. , Review of particle physics: Particle data group . Phys. Rev. D. 98 , 030001 ( 2018 ). doi: 10.1103/PhysRevD.98.030001 http://doi.org/10.1103/PhysRevD.98.030001 Review of particle physics: Particle data group Z.W. Pan , J.Y. Dong , X.J. Ni et al. , Conceptual design and update of the 128-channel μSR prototype spectrometer based on musrSim . Nucl. Sci. Tech. 30 , 123 ( 2019 ). doi: 10.1007/s41365-019-0648-5 http://doi.org/10.1007/s41365-019-0648-5 Conceptual design and update of the 128-channel μSR prototype spectrometer based on musrSim L.P. Zhou , Q.L. Mu , H.T. Jing et al. , A possible scheme for the surface muon beamline at CSNS . Nucl. Sci. Tech. 30 , 169 ( 2019 ). doi: 10.1007/s41365-019-0684-1 http://doi.org/10.1007/s41365-019-0684-1 A possible scheme for the surface muon beamline at CSNS R.H. Bernstein , P.S. Cooper , Charged lepton flavor violation: An experimenter's guide . Phys. Rep. 532 , 27 - 64 ( 2013 ). doi: 10.1016/j.physrep.2013.07.002 http://doi.org/10.1016/j.physrep.2013.07.002 Charged lepton flavor violation: An experimenter's guide Y. Kuno , A search for muon-to-electron conversion at J-PARC: the COMET experiment . Prog. Theor. Phys. 2013 , 2 , 022C01 ( 2013 ). doi: 10.1093/ptep/pts089 http://doi.org/10.1093/ptep/pts089 A search for muon-to-electron conversion at J-PARC: the COMET experiment K.N. Borozdin , G.E. Hogan , C. Morris et al. , Surveillance: Radiographic imaging with cosmic-ray muons . Nature 422 , 277 ( 2003 ). doi: 10.1038/422277a http://doi.org/10.1038/422277a Surveillance: Radiographic imaging with cosmic-ray muons P.L. Rocca , V. Antonuccio , M. Bandieramonte et al. , Search for hidden high- Z materials inside containers with the Muon Portal Project . J. Instrum. 9 , C01056 ( 2014 ). doi: 10.1088/1748-0221/9/01/C01056 http://doi.org/10.1088/1748-0221/9/01/C01056 Search for hidden high-Z materials inside containers with the Muon Portal Project J. Armitage , J. Botte , K. Boudjemline et al. , First images from the cript muon tomography system . Int. J. Mod. Phys. Conf. Ser 27 , 1460129 ( 2014 ). doi: 10.1142/S201019451460129X http://doi.org/10.1142/S201019451460129X First images from the cript muon tomography system K. Borozdin , S. Greene , Z. Lukic et al. , Cosmic ray radiography of the damaged cores of the Fukushima reactors . Phys. Rev. Lett. 109 , 152501 ( 2012 ). doi: 10.1103/physrevlett.109.152501 http://doi.org/10.1103/physrevlett.109.152501 Cosmic ray radiography of the damaged cores of the Fukushima reactors X.Y. Pan , Y.F. Zheng , Z. Zeng et al. , Experimental validation of material discrimination ability of muon scattering tomography at the TUMUTY facility . Nucl. Sci. Tech. 30 , 120 ( 2019 ). doi: 10.1007/s41365-019-0649-4 http://doi.org/10.1007/s41365-019-0649-4 Experimental validation of material discrimination ability of muon scattering tomography at the TUMUTY facility S. Xiao , W.B. He , M.C. Lan et al. , A modified multi-group model of angular and momentum distribution of cosmic ray muons for thickness measurement and material discrimination of slabs . Nucl. Sci. Tech. 29 , 28 ( 2018 ). doi: 10.1007/s41365-018-0363-7 http://doi.org/10.1007/s41365-018-0363-7 A modified multi-group model of angular and momentum distribution of cosmic ray muons for thickness measurement and material discrimination of slabs J.N. Dong , Y.L. Zhang , Z.Y. Zhang et al. , Position-sensitive plastic scintillator detector with WLS-fiber readout . Nucl. Sci. Tech. 29 , 117 ( 2018 ). doi: 10.1007/s41365-018-0449-2 http://doi.org/10.1007/s41365-018-0449-2 Position-sensitive plastic scintillator detector with WLS-fiber readout L.W. Alvarez , J.A. Anderson , F.E. Bedwei et al . Search for hidden chambers in the pyramids . Science. 167 , 3919 ( 1970 ). doi: 10.1126/science.167.3919.832 http://doi.org/10.1126/science.167.3919.832 Search for hidden chambers in the pyramids HKM. Tanaka , T. Uchida , M. Tanaka et al . Cosmic-ray muon imaging of magma in a conduit: Degassing process of Satsuma-Iwojima Volcano . GEOPHYS. RES. LETT. 36 , L01304 ( 2009 ). doi: 10.1029/2008GL036451 http://doi.org/10.1029/2008GL036451 Cosmic-ray muon imaging of magma in a conduit: Degassing process of Satsuma-Iwojima Volcano K. Jourde , D. Gibert , J. Marteau et al. , Muon dynamic radiography of density changes induced by hydrothermal activity at the La Soufrière of Guadeloupe volcano . Sci. Rep. 6 , 33406 ( 2016 ). doi: 10.1038/srep33406 http://doi.org/10.1038/srep33406 Muon dynamic radiography of density changes induced by hydrothermal activity at the La Soufrière of Guadeloupe volcano V. Tioukov , A. Alexandrov , C. Bozza et al. , First muography of Stromboli volcano . Sci. Rep. 9 , 6695 ( 2019 ). doi: 10.1038/s41598-019-43131-8 http://doi.org/10.1038/s41598-019-43131-8 First muography of Stromboli volcano N. Lesparre , J. Marteau , Y. Declais et al. , Design and operation of a field telescope for cosmic ray geophysical tomography . Geosci, Instrum. Method. Data Syst. 1 , 33 ( 2012 ). doi: 10.5194/gi-1-33-2012 http://doi.org/10.5194/gi-1-33-2012 Design and operation of a field telescope for cosmic ray geophysical tomography C. Carloganu , V. Niess , S. Bene et al. , Towards a muon radiography of the Puy de Dme . Geosci. Instrum, Method. Data Syst. 2 , 55 ( 2013 ). doi: 10.5194/gi-2-55-2013 http://doi.org/10.5194/gi-2-55-2013 Towards a muon radiography of the Puy de Dme A. Anastasio , F. Ambrosino , D. Basta et al. , The MU-RAY experiment. An application of SiPM technology to the understanding of volcanic phenomena . Nucl. Instru. and Meth. A 718 , 134 - 137 ( 2013 ). doi: 10.1016/j.nima.2012.08.065 http://doi.org/10.1016/j.nima.2012.08.065 The MU-RAY experiment. An application of SiPM technology to the understanding of volcanic phenomena VV. Kobylyansky , VV. Moiseichenko , AI. Myagkikh et al. , The deep-sea Cerenkov muon detector with spatial structure . Sov. J. Phys. Oceanogr. 3 , 149 - 153 ( 1992 ). doi: 10.1007/BF02197621 http://doi.org/10.1007/BF02197621 The deep-sea Cerenkov muon detector with spatial structure N. Lesparre , D. Gibert , J. Marteau et al. , Geophysical muon imaging: feasibility and limits . Geophys. J. Int. 183 , 1348 - 1361 ( 2010 ). doi: 10.1111/j.1365-246X.2010.04790.x http://doi.org/10.1111/j.1365-246X.2010.04790.x Geophysical muon imaging: feasibility and limits K. Morishima , M. Kuno , A. Nishio et al. , Discovery of a big void in Khufu's Pyramid by observation of cosmic-ray muons . Nature 552 , 386 - 390 ( 2017 ). doi: 10.1038/nature24647 http://doi.org/10.1038/nature24647 Discovery of a big void in Khufu's Pyramid by observation of cosmic-ray muons R. Han , Q. Yu , Z. Li et al . Cosmic muon flux measurement and tunnel overburden structure imaging . J. Instrum. 15 , P06019 ( 2020 ). doi: 10.1088/1748-0221/15/06/P06019 http://doi.org/10.1088/1748-0221/15/06/P06019 Cosmic muon flux measurement and tunnel overburden structure imaging K. Nagamine , M. Iwasaki , K. Shimomura et al. , Method of probing inner-structure of geophysical substance with the horizontal cosmic-ray muons and possible application to volcanic eruption prediction . Nucl. Instru. Meth. A 356 , 585 - 595 ( 1995 ). doi: 10.1016/0168-9002(94)01169-9 http://doi.org/10.1016/0168-9002(94)01169-9 Method of probing inner-structure of geophysical substance with the horizontal cosmic-ray muons and possible application to volcanic eruption prediction A. Anastasio , F. Ambrosino , D. Basta et al. , The MU-RAY detector for muon radiography of volcanoes . Nucl. Instru. Meth. A 732 , 423 - 426 ( 2013 ). doi: 10.1016/j.nima.2013.05.159 http://doi.org/10.1016/j.nima.2013.05.159 The MU-RAY detector for muon radiography of volcanoes S. Pesente , S. Vanini , M. Benettoni et al. , First results on material identification and imaging with a large-volume muon tomography prototype . Nucl. Instru. and Meth. A 604 , 3 ( 2009 ). doi: 10.1016/j.nima.2009.03.017 http://doi.org/10.1016/j.nima.2009.03.017 First results on material identification and imaging with a large-volume muon tomography prototype K. Gnanvo , L.V. Grasso , M. Hohlmann et al. , Imaging of high-Z material for nuclear contraband detection with a minimal prototype of a Muon Tomography station based on GEM detectors . Nucl. Instru. Meth. A 652 , 16 - 20 ( 2011 ). doi: 10.1016/j.nima.2011.01.163 http://doi.org/10.1016/j.nima.2011.01.163 Imaging of high-Z material for nuclear contraband detection with a minimal prototype of a Muon Tomography station based on GEM detectors S. Bouteille , D. Attie , P. Baron et al. , A Micromegas-based telescope for muon tomography: The WatTo experiment . Nucl. Instrum. Meth. A 834 , 223 - 228 ( 2016 ). doi: 10.1016/j.nima.2016.08.002 http://doi.org/10.1016/j.nima.2016.08.002 A Micromegas-based telescope for muon tomography: The WatTo experiment J. Feng , Z. Zhang , J. Liu et al. , A thermal bonding method for manufacturing Micromegas detectors . Nucl. Instru. Meth. A 989 , 164958 ( 2020 ). doi: 10.1016/j.nima.2020.164958 http://doi.org/10.1016/j.nima.2020.164958 A thermal bonding method for manufacturing Micromegas detectors I. Bikit , D. Mrdja , K. Bikit et al. , Novel approach to imaging by cosmic-ray muons . EPL 113 , 58001 ( 2016 ). doi: 10.1209/0295-5075/113/58001 http://doi.org/10.1209/0295-5075/113/58001 Novel approach to imaging by cosmic-ray muons D. Mrdja , I. Bikit , K. Bikit et al. , First cosmic-ray images of bone and soft tissue . EPL 116 , 48003 ( 2016 ). doi: 10.1209/0295-5075/116/48003 http://doi.org/10.1209/0295-5075/116/48003 First cosmic-ray images of bone and soft tissue G. Galgoczi , D. Mrdja , I. Bikit et al. , Imaging by muons and their induced secondary particles—a novel technique . J. Instrum. 15 , 6 ( 2020 ). doi: 10.1088/1748-0221/15/06/C06014 http://doi.org/10.1088/1748-0221/15/06/C06014 Imaging by muons and their induced secondary particles—a novel technique A.G. Bogdanov , H. Burkhardt , V.N. Ivanchenko et al . Geant4 simulation of production and interaction of Muons . IEEE Trans. Nucl. Sci. 53 , 2 ( 2006 ). doi: 10.1109/TNS.2006.872633 http://doi.org/10.1109/TNS.2006.872633 Geant4 simulation of production and interaction of Muons All authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by Xuan-Tao Ji, Si-Yuan Luo, Yu-He Huang and Xiao-Dong Wang. The first draft of the manuscript was written by Xuan-Tao Ji and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript. 10.1007/s41365-022-00989-0 This work was supported by the Ministry of Science and Technology of China Foundation (No. 2020YFE0202001), the National Natural Science Foundation of China (No. 11875163), and the Natural Science Foundation of Hunan Province (No. 2021JJ20006). 2021-08-12 2021-09-08 2021-09-26 2022-01-15 2022-02-16 © China Science Publishing & Media Ltd. (Science Press), Shanghai Institute of Applied Physics, the Chinese Academy of Sciences, ChineseNuclear Society and Springer Nature Singapore Pte Ltd. 2022 This is a post-peer-review, pre-copyedit version of an article published in Nuclear Science and Techniques. The final authenticated version is available online at: http://dx.doi.org/10.1007/s41365-022-00989-0 http://dx.doi.org/10.1007/s41365-022-00989-0 . 2022 Nuclear Science and Techniques 01 33