Imaging internal density structure of the Laoheishan volcanic cone with cosmic ray muon radiography

Imaging internal density structure of the Laoheishan volcanic cone with cosmic ray muon radiography 1 corresp first-author hanran@ncepu.edu.cn hanran@ncepu.edu.cn hanran@ncepu.edu.cn https://orcid.org/0000-0001-9717-6958 1 https://orcid.org/0000-0002-5566-8193 1 corresp zwli@whigg.ac.cn zwli@whigg.ac.cn zwli@whigg.ac.cn 1 1 2 2 2 3 2 4 Science and Technology on Reliability and Environmental Engineering Laboratory, Beijing Institute of Spacecraft Environment Engineering, Beijing 100094, China Science and Technology on Reliability and Environmental Engineering Laboratory, Beijing Institute of Spacecraft Environment Engineering, Beijing 100094, China School of Nuclear Science and Engineering, North China Electric Power University, 2 BeiNong Road, Beijing 102206, China School of Nuclear Science and Engineering, North China Electric Power University, 2 BeiNong Road, Beijing 102206, China College of nuclear science and technology, Beijing Normal University, 19 Xinjiekou outer St., Beijing 100875, China College of nuclear science and technology, Beijing Normal University, 19 Xinjiekou outer St., Beijing 100875, China Aerospace Information Research Institute, Chinese Academy of Sciences, 9 Dengzhuang South Road, Beijing 100094, China Aerospace Information Research Institute, Chinese Academy of Sciences, 9 Dengzhuang South Road, Beijing 100094, China Muon radiography is a promising technique for imaging the internal density structures of targets such as tunnels, pyramids, and volcanoes up to a scale of a few hundred meters by measuring the flux attenuation of cosmic ray muons after they have travelled through these targets. In this study, we conducted experimental muon radiography of one of the volcanoes in the Wudalianchi area in Northeast China to image its internal density structure. The muon detector used in this study was composed of plastic scintillators and silicon photomultipliers. After approximately one and a half months of observing the crater and conduit of the Laoheishan volcano cone in Wudalianchi from September 23 rd to November 10 th 2019, more than 3 million muon tracks fulfilling the data selection criteria were collected. Based on the muon samples and high-resolution topography obtained through aerial photogrammetry using an unmanned aerial vehicle, a density image of the Laoheishan volcano cone was constructed. The results obtained in this experiment demonstrate the feasibility of using a radiography technique based on plastic scintillator detectors. To obtain the density distribution, we performed a detailed background analysis and found that low-energy charged particles dominated the background noise. Relatively higher densities were found near the surface of the volcanic cone, whereas relatively lower densities were found near the center of the volcanic cone. The experiment in this study is the first volcano muon tomography study performed in China. Our work provides an important reference for future research. Muon radiography Muon transmission imaging Density Introduction 1 Volcanic hazard assessment and risk management are important for both population safety and economic development. On January 15, 2022, a powerful explosive eruption from the underwater volcano Hunga Tonga-Hunga Ha'apai brought powerful tsunami waves and heavy ashfall to islands in Tonga, severely damaging coastal communities. Understanding the internal structures of volcanos is essential for forecasting such volcanic hazards. Conventional geophysical methods have spatial resolutions that typically range from tens of meters to 1 km. Such low resolutions are insufficient for detecting small volumes of magma or magma conduits [ 1 1 ]. Gravimetry is a commonly used geophysical method. However, it suffers from ill-posed inversion matrices, and its sensitivity is limited to relatively shallow depths [ 2 2 ]. As early as 1994, Nagamine et al. proposed the use of cosmic ray muons to probe the internal structures of volcanoes [ 3 3 ]. To overcome the limitation that X-ray radiography cannot be applied to medium-to-large and dense objects, Borozdin et. al. proposed in 2003 that natural background muons with high penetrating power can be used for radiographic imaging [ 4 4 ]. This proposal has attracted broad attention. The use of global density tomography for volcano monitoring was subsequently promoted in [ 5 5 - 7 7 5-7 ]. Muon radiography can provide high-resolution density measurements of shallow structures with a density contrast precision of 1–3%, which is significantly higher than that achievable using conventional geophysical methods such as electromagnetic and seismic techniques [ 5 5 ]. For example, a conduit with a diameter of 102 ± 15 m was successfully imaged using muon radiography, as reported in [ 8 8 ]. In 2010, Jack et al. theoretically investigated the possibility of using muon radiography to determine the internal density distributions of candidate target structures [ 9 9 ]. In this paper, we present our work on the use of muon radiography to perform volcano imaging and to investigate the inner structures of a volcano. The remainder of this paper is structured as follows: First, we provide an overview of muon radiography (Sect. 2). We then provide a detailed description of the muon radiography project at Wudalianchi (Sect. 3). In Sect. 4, we present the data processing method as well as the imaging procedure and the obtained density results. We explain the possible contamination from background events involving charged particles. We finally provide concluding remarks and a discussion in Sect. 5. Muon Radiography 1 Principle 2 Muons are produced through the interaction of primary cosmic rays with the Earth's atmosphere. Muons have drawn much attention because of their utility in an imaging method called muon radiography, which is similar to X-ray radiography. In addition, muon imaging based on multiple muon scattering is also a popular topic in the field of radiation detection imaging and is commonly used for the identification of high- Z materials [ 10 10 ]. As muons pass through the investigated object, they lose energy via interactions with matter such as ionization and radiation. There are different attenuation coefficients for gamma rays and X-rays passing through different materials [ 11 11 ]. The muon flux attenuation due to interaction with matter should be carefully modelled according to the spectrum of the incident cosmic muons. In this study, we used a muon flux model named the modified Gaisser's formula [ 12 12 ]. Compared to Gaisser's formula [ 13 13 ], which is applicable in the high-energy range (100 GeV–100 TeV) and at small zenith angles, the modified Gaisser's formula is extended to low energies and large zenith angles and is given in Eq. 1: <math display="block" id="M1"><mtable><mtr><mtd><mi>Φ</mi><mo stretchy="false">(</mo><mi>θ</mi><mn>,</mn><msub><mi>E</mi><mi>μ</mi></msub><mo stretchy="false">)</mo><mo>=</mo><mn>0.14</mn><msup><mrow><mo>[</mo><mrow><mfrac><mrow><msub><mi>E</mi><mi>μ</mi></msub></mrow><mrow><mtext>GeV</mtext></mrow></mfrac><mrow><mo>(</mo><mrow><mn>1</mn><mo>+</mo><mfrac><mrow><mn>3.64</mn><mtext> </mtext><mtext>GeV</mtext></mrow><mrow><msub><mi>E</mi><mi>μ</mi></msub><msup><mrow><mo stretchy="false">(</mo><mi>cos</mi><msup><mi>θ</mi><mo>*</mo></msup><mo stretchy="false">)</mo></mrow><mrow><mn>1.29</mn></mrow></msup></mrow></mfrac></mrow><mo>)</mo></mrow></mrow><mo>]</mo></mrow><mrow><mo>−</mo><mn>2.7</mn></mrow></msup></mtd></mtr><mtr><mtd><mrow><mo>(</mo><mrow><mfrac><mn>1</mn><mrow><mn>1</mn><mo>+</mo><mfrac><mrow><mn>1.1</mn><mi>E</mi><mi>cos</mi><msup><mi>θ</mi><mo>*</mo></msup></mrow><mrow><mn>115</mn><mtext> </mtext><mtext>GeV</mtext></mrow></mfrac></mrow></mfrac><mo>+</mo><mfrac><mrow><mn>0.054</mn></mrow><mrow><mn>1</mn><mo>+</mo><mfrac><mrow><mn>1.1</mn><mi>E</mi><mi>cos</mi><msup><mi>θ</mi><mo>*</mo></msup></mrow><mrow><mn>850</mn><mtext> </mtext><mtext>GeV</mtext></mrow></mfrac></mrow></mfrac></mrow><mo>)</mo></mrow><mn>,</mn></mtd></mtr></mtable></math> where E μ is the muon energy. The quantity cos θ * is related to the zenith angle θ by: <math display="block" id="M2"><mrow><mi>cos</mi><msup><mi>θ</mi><mo>*</mo></msup><mo>=</mo><mrow><mo>(</mo><mrow><mfrac><mrow><msup><mstyle mathsize="140%" displaystyle="true"><mrow><mi>cos</mi></mrow></mstyle><mn>2</mn></msup><mi>θ</mi><mo>+</mo><msubsup><mi>p</mi><mn>1</mn><mn>2</mn></msubsup><mo>+</mo><msub><mi>p</mi><mn>2</mn></msub><msup><mrow><mo stretchy="false">(</mo><mi>cos</mi><mi>θ</mi><mo stretchy="false">)</mo></mrow><mrow><msub><mi>p</mi><mn>3</mn></msub></mrow></msup><mo>+</mo><msub><mi>p</mi><mn>4</mn></msub><mi>cos</mi><msup><mi>θ</mi><mrow><msub><mi>p</mi><mn>5</mn></msub></mrow></msup></mrow><mrow><mn>1</mn><mo>+</mo><msubsup><mi>p</mi><mn>1</mn><mn>2</mn></msubsup><mo>+</mo><msub><mi>p</mi><mn>2</mn></msub><mo>+</mo><msub><mi>p</mi><mn>4</mn></msub></mrow></mfrac></mrow><mo>)</mo></mrow><mn>,</mn></mrow></math> where p 1 , p 2 , p 3 , p 4 , and p 5 are the parameters given in Table 1 in [ 14 14 ] and have the values of 0.102573, –0.068287, 0.958633, 0.0407253, and 0.817285, respectively. The zenith angle θ =0° corresponds to the vertical direction perpendicular to the ground while θ =90° represents the horizontal direction parallel to the ground. When there are no visible obstructions in the muon path from the sky to the detector, that is, in the open-sky case, the expected number of muons can be predicted as: <math display="block" id="M3"><mrow><msub><mi>N</mi><mrow><mtext>opensky</mtext></mrow></msub><mo stretchy="false">(</mo><mi>θ</mi><mo stretchy="false">)</mo><mo>=</mo><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mrow><mn>0.1</mn></mrow><mi>∞</mi></msubsup></mrow></mstyle><mi>Φ</mi><mo stretchy="false">(</mo><mi>θ</mi><mn>,</mn><mi>E</mi><mo stretchy="false">)</mo><msub><mi>S</mi><mrow><mtext>eff</mtext></mrow></msub><mo stretchy="false">(</mo><mi>θ</mi><mo stretchy="false">)</mo><mi>ϵ</mi><mi>Δ</mi><mi>T</mi><mtext>d</mtext><mi>E</mi><mn>,</mn></mrow></math> where S eff ( θ ) denotes the effective area of the detector, ϵ is the total efficiency in the experiment, Δ T is the data acquisition time, and Φ ( θ, E ) is the differential muon flux as a function of the zenith angle θ and muon energy E , as shown in Eq. 1. The integration starts from 0.1 GeV, which was determined based on experience. When there is a volcano, the expected number of muons is <math display="block" id="M4"><mrow><msub><mi>N</mi><mrow><mtext>volcano</mtext></mrow></msub><mo stretchy="false">(</mo><mi>θ</mi><mo stretchy="false">)</mo><mo>=</mo><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mrow><msub><mi>E</mi><mrow><mtext>min</mtext></mrow></msub></mrow><mi>∞</mi></msubsup></mrow></mstyle><mi>Φ</mi><mo stretchy="false">(</mo><mi>θ</mi><mn>,</mn><mi>E</mi><mo stretchy="false">)</mo><msub><mi>S</mi><mrow><mtext>eff</mtext></mrow></msub><mo stretchy="false">(</mo><mi>θ</mi><mo stretchy="false">)</mo><mi>ϵ</mi><mi>Δ</mi><mi>T</mi><mtext>d</mtext><mi>E</mi><mn>,</mn></mrow></math> where E min denotes the minimum energy required for a muon to cross the investigated object without being absorbed. E min is dependent on the density-length, which is the product of the average material density ρ and path length l inside the investigated object traversed by the muon. In this study, the density-length is calculated as the product of the volcano rock density and the thickness encountered by the muon along a given line of the muon track before it enters the detector. In this method, the path is divided into many small line segments called steps. Whether the step is inside the volcano or not is determined by the z -coordinate of the line segment. If the z -coordinate is smaller than the elevation of the volcano, the step is accumulated. Otherwise, the track is located outside the volcano. The thickness of the volcano encountered along the muon track can thus be derived. For a given density-length, there is a minimum muon energy ( E min ) required, which can be estimated using either the empirical Bethe-Bloch formula [ 15 15 ] or the numerical values provided by the Particle Data Group [ 16 16 ]. The continuous slowing down ability (CSDA) table by [ 17 17 ] gives the relationship between the initial muon kinetic energy and the average penetration depth in various materials. We show the penetration depth for water (red solid line) in Fig. 1 Fig. 1 ). Based on the modified Gaisser's formula (Eq. 1), the percentage of muons with energies larger than the minimum muon energy (muon survival ratio) can be calculated. By measuring the muon flux attenuation, E min and the density-length of the investigated object can be deduced. Together with the path length obtained from geodetic measurements, the average density of the research target can then be determined. This is the principle behind muon radiography. Muon Detector 2 With the development of particle detection techniques, particle detectors have become sufficiently precise to perform reliable measurements of muon trajectories. Various tracking detectors such as gaseous detectors, nuclear emulsions, and plastic scintillator detectors have been used for volcano density investigations. Gaseous detectors, such as resistive plate chambers (RPC) [ 18 18 ], are widely used to detect high-energy charged particles, particularly muons. Typical RPCs have excellent det ection efficiency as well as timing and spatial resolutions that are on the order of nanoseconds and centimeters, respectively. In addition, they can provide a large detection area at a low cost per unit area. Emulsion detectors ([ 8 8 ]) have high angular and spatial resolutions. In addition, they are portable and do not require electric power. The plastic scintillator detector ([ 19 19 , 20 20 , 7 7 ]) is another widely used particle detection device in nuclear and particle physics. The detection principle of a plastic scintillator detector is as follows: When the detector is hit by charged particles or radiation, the plastic scintillator bar emits flashes of scintillation light. Through coupling to photomultipliers or silicon photomultipliers, the scintillation light can be converted into electric signals. The fired scintillation bar is tagged in the same manner. Plastic scintillator detectors offer good time resolution on the order of nanoseconds and relatively high light output. They are easily produced commercially and relatively inexpensive. In addition, their flexibility and portable assembly offer great advantages for investigating volcanic density. The muon telescope in our experiments was a three-layer plastic scintillator detector. Each detection layer has N x = N y = 16 scintillation bars with the dimensions of 80 cm× 5 cm ×1 cm. There are two subplanes in each layer placed orthogonal to each other to provide the x and y coordinates of the fired spot caused by energy deposition from a muon. A total of 256 pixels with the detection area size of 25 cm 2 are formed. The separation between the top and bottom panels is 100 cm. The readout electronics and details of the detector can be found in [ 21 21 , 22 22 ]. Muon radiography at Wudalianchi 1 Detector placement 2 Our research target Laoheishan is located in Heihe, Heilongjiang Province, China. Laoheishan is one of the volcanoes of Wudalianchi. The volcano was measured to be approximately 1400 m in diameter and 200 m in height. Wudalianchi volcanoes (referring to Laoheishan unless otherwise noted) are very interesting because the main stages in the geological effects of ongoing geomorphological evolution can be studied in these volcanoes, and they are the youngest volcanoes in China. The last eruption occurred approximately 300 years ago [ 23 23 ]. Before the measurements, we conducted a geological survey using an UAV. A high-precision 3D image of the volcano was obtained ( Fig. 2a Fig. 2a ). The horizontal resolution of the drone image is 1 cm and the vertical resolution 10 cm. A plastic scintillator detector was deployed at a distance of about 250 m from the Laoheishan volcano. From the 3D image of the volcano, we calculated the path length of the muon trajectory. The calculated path lengths are presented in Fig. 2b Fig. 2b . The path length was obtained from the detector point of view and is given in the geological frame in which the horizontal axis denotes the azimuth angle and the vertical axis denotes the elevation angle. The elevation angle is 0 at ground level and 90 directly overhead. The azimuth angle is 0 at true north in the horizontal plane, and ranges from 0 to 180 clockwise, and from 0 to -180 anticlockwise. The path lengths are represented by different colors and the 100 m contour lines are superimposed on them. The path length of the volcano is a crucial input for volcano density devolution. Based on the path length image, the detector was elevated to 20 to improve the scanning of the volcano. The acceptance was thus increased and the time for data taking shortened. The detector was located at an altitude of 400 m. The highest part of the volcano is approximately 600 m; therefore, the relative height of the volcano with respect to the detector was approximately 200 m. The detector angle of the volcano was a cone with an opening angle of approximately 43. The angular resolution of the plastic scintillator detector is 50 mrad and its spatial resolution is a few tens of meters at a distance of 250 m. Data Taking 2 The measurements were performed between September 23 rd and November 10 th 2019 for approximately 68 days at the Laoheishan site. More than 3 million muons were collected. The data presented in this study were obtained over three experimental stages: First, the detector was oriented 10 west of north to perform volcano imaging ( Fig. 3a Fig. 3a ). The detector was then rotated to the true west to record the muon flux from the open sky ( Fig. 3b Fig. 3b ). The open-sky measurements were used for comparison with the theoretical predictions from the modified Gaisser formula. The measured fluxes should agree well with the modelled flux. Finally, the detector was oriented 30 east of north to cover the missing part of the volcano ( Figs. 3c, d) Figs. 3c, d) corresponding to the part outside the field of view of the detector in the first stage. The fields of view of the detector facing the volcano in the two stages are shown in Fig. 2b Fig. 2b . Data were collected over a period of approximately one day. All the muon events were recorded for later analysis. Data Analysis 1 Data Processing Method 2 The following data selection method was used to determine the muon tracks in our measurements: Muon tracks were selected from the event samples in which with all three detection panels fired. The muon trajectories are supposed to be straight lines and were determined using the following line selection criteria: <math display="block" id="M5"><mtable columnalign="left"><mtr><mtd><mrow><mo>|</mo><mrow><mo stretchy="false">(</mo><msub><mi>x</mi><mn>1</mn></msub><mo>−</mo><msub><mi>x</mi><mn>2</mn></msub><mo stretchy="false">)</mo><mo>−</mo><mo stretchy="false">(</mo><msub><mi>x</mi><mn>2</mn></msub><mo>−</mo><msub><mi>x</mi><mn>3</mn></msub><mo stretchy="false">)</mo></mrow><mo>|</mo></mrow><mo><</mo><mn>1</mn></mtd></mtr><mtr><mtd><mrow><mo>|</mo><mrow><mo stretchy="false">(</mo><msub><mi>y</mi><mn>1</mn></msub><mo>−</mo><msub><mi>y</mi><mn>2</mn></msub><mo stretchy="false">)</mo><mo>−</mo><mo stretchy="false">(</mo><msub><mi>y</mi><mn>2</mn></msub><mo>−</mo><msub><mi>y</mi><mn>3</mn></msub><mo stretchy="false">)</mo></mrow><mo>|</mo></mrow><mo><</mo><mn>1</mn></mtd></mtr></mtable></math> where 1 on the right-most side denotes the ID difference between the scintillation bars, and ( x 1 , y 1 ), ( x 2 , y 2 ), and ( x 3 , y 3 ) are the coordinates of the fired spots in each detection panel. Euler rotation was performed to obtain the zenith and azimuth coordinates ( θ g , ϕ g ) (subscript g denotes geological) of the muon tracks in the geological coordinate frame. The muon track in the geological coordinate frame is thus derived. To evaluate the muon flux, various corrections were made to remove the impact of the detector response. Efficiency correction was performed on the plastic scintillator bars. The overall efficiency correction factor for each scintillator bar was determined from the vertical muon samples. Because there are three detection panels, the efficiency of each scintillation bar can be derived by dividing the number of events registered by each bar by the number of coincident events from the remaining two bars. The overall efficiency includes the contributions from the light yield of the plastic scintillator, the photon detection efficiency of the silicon photomultiplier (SiPM), and the threshold of the electronics. The overall efficiency value extracted for each scintillator bar are presented in Fig. 4a Fig. 4a . The average efficiency of the 96 scintillator bars is approximately 85%. There were two scintillator bars with extremely low efficiencies of less than 60%, which might have been caused by poor connection to the SiPM. Because the measured number of muon events is proportional to the solid angle and detection area, it is necessary to perform solid angle and acceptance correction to obtain the absolute muon flux. The muon trajectory is restricted by the detection pixel pair in the top and bottom panels. The solid angle subtended by the muon trajectory from the detection pixel on the top panel to the detection pixel on the bottom panel was calculated following the method introduced in [ 24 24 ] and shown in Fig. 4b Fig. 4b . In Fig. 4b Fig. 4b , θ x ( θ y ) was defined as atan(5×Δ x/d ) (atan(5×Δ y/d )). Δ x (Δ y ) is the ID difference between the detection pixels along the x ( y ) direction and d is the distance between the top and bottom detection panels. The numerical value of 5 is the width of the scintillation bars. All the values are given in centimeters. The solid angle decreases rapidly with increasing θ x or θ y . For vertical muons (i.e., θ x , θ y = 0), the solid angle spanned by each detection pixel is approximately 2.5×10 -3 sr. The detector acceptance is introduced to account for differences in the effective detection area of the muon detector for muons with different injection directions. For example, for muons injected perpendicularly into the detection panel, there are 16×16 detection pixels with a detection area of 25 cm 2 for each detection pixel, forming an acceptance of approximately 16 cm 2 · sr. The acceptance function (cm 2 · sr) corresponding to the detector configuration is shown in Fig. 4c Fig. 4c . The acceptance also decreases quickly with increasing θ x or θ y . The muon flux was measured by applying the above-mentioned corrections, and the results are presented in Fig. 5 Fig. 5 . The measured results as a function of the zenith angle are shown in Fig. 5a Fig. 5a . The open-sky muon flux (red dots) and its comparison with the prediction from the modified Gaisser formula (solid blue line) are also presented. The discrepancy with the latter is probably due to contamination from muon-like particles, such as electrons and protons, and will be discussed later. The pixel-by-pixel muon flux measurement results in the open-sky and volcano cases for each detection pixel are also shown. Figure 5b Figure 5b shows the results in the open-sky case. The muon flux decreased rapidly with the elevation angle. For elevation angles of less than 20, the flux was smaller than 0.001. Figure 5c Figure 5c shows the same variable for the volcano case. A comparison between the open-sky flux measurement and the modified Gassier prediction was also performed. There is good overall agreement with ratio values varying from 0.9 to 1.2, as shown in Fig. 5d Fig. 5d . This proves that the muon detector, including its solid angle and acceptance, was understood correctly. Another part of the data analysis is volcano density analysis. For this purpose, the muon transmission power is defined as the ratio of the muon event rates in the volcano and open-sky cases, as shown in the first line of Eq. 5. By substituting Eq. 2 and Eq. 3 into the first line of Eq. 5, the muon transmission power can be further rewritten as the second line of Eq. 5. <math display="block" id="M6"><mtable><mtr><mtd><msup><mi>K</mi><mi>α</mi></msup><mo>=</mo><mfrac><mrow><msub><mi>N</mi><mtext>V</mtext></msub><mo stretchy="false">(</mo><mi>θ</mi><mo stretchy="false">)</mo><mo>/</mo><mi>Δ</mi><msub><mi>T</mi><mtext>V</mtext></msub></mrow><mrow><msub><mi>N</mi><mtext>O</mtext></msub><mo stretchy="false">(</mo><mi>θ</mi><mo stretchy="false">)</mo><mo>/</mo><mi>Δ</mi><msub><mi>T</mi><mtext>O</mtext></msub></mrow></mfrac></mtd></mtr><mtr><mtd><mo>=</mo><mfrac><mrow><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mrow><msub><mi>E</mi><mrow><mtext>min</mtext></mrow></msub></mrow><mi>∞</mi></msubsup></mrow></mstyle><mi>Φ</mi><mo stretchy="false">(</mo><mi>θ</mi><mn>,</mn><mi>E</mi><mo stretchy="false">)</mo><msub><mi>S</mi><mrow><mtext>eff</mtext></mrow></msub><mo stretchy="false">(</mo><mi>θ</mi><mo stretchy="false">)</mo><mi>ϵ</mi><mtext>d</mtext><mi>E</mi></mrow><mrow><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mrow><mn>0.1</mn></mrow><mi>∞</mi></msubsup></mrow></mstyle><mi>Φ</mi><mo stretchy="false">(</mo><mi>θ</mi><mn>,</mn><mi>E</mi><mo stretchy="false">)</mo><msub><mi>S</mi><mrow><mtext>eff</mtext></mrow></msub><mo stretchy="false">(</mo><mi>θ</mi><mo stretchy="false">)</mo><mi>ϵ</mi><mtext>d</mtext><mi>E</mi></mrow></mfrac></mtd></mtr></mtable></math> Here, the subscript V(O) denotes the volcano (open-sky) case. By integrating the muon energy spectrum given by the modified Gaisser formula in Eq. 1, E min , the minimum energy the muon must possess to cross the volcano without being absorbed, could be obtained. The CSDA table from [ 17 17 ], which lists the CSDA values for various materials, was used to evaluate the energy deposited by the muon when it crossed the volcano. Here, we used the numerical values for standard rock. The functional relationship between the density-length ρ.l and muon transmission power K α , that is, the muon survival ratio, can then be deduced. The results are presented in Fig. 1 Fig. 1 and used for the volcano density analysis. Muon radiography results 2 Muon radiography can be used to obtain the profile of the volcano through the muon flux ratio plot. The flux ratio plot was obtained from the pixel-by-pixel ratios of the flux measurement results for the volcano and open-sky cases. In addition, detector-related corrections, such as the efficiency correction for the scintillation bars and the solid angle and acceptance corrections, can be cancelled out naturally in the pixel-by-pixel ratio calculation. The muon flux ratios are presented in Fig. 6a Fig. 6a . The black curves represent the contour lines of the path length in the volcano. The volcano boundary is denoted by the outermost black curve. The muon flux ratio is represented by different colors in which blue represents ratio values above 0.6. There is good agreement between the superposition of the flux ratio plot and the volcano path length. The impact of the volcano can be observed from the ratio plot. Muons with energies less than E min were stopped inside the volcano, resulting in a small ratio value, which reflects the average density of the volcano along the specific direction of the muon track. The large ratio values in the regions without the volcano indicate that the muons passed through only air, as shown by the blue part in Fig. 6a Fig. 6a . In this way, the profile image of the Laoheishan volcano was obtained using the muon telescope. The inner density structure of the volcano is the research focus of this study. Based on the muon flux ratio analysis, the density was unfolded through the following procedures: The flux ratio value for a given azimuth and elevation angle was first obtained. Using the modified Gaisser formula, the minimum energy required to travel through the volcano at the given zenith angle was then calculated. The density-length could subsequently be deduced from the stopping power curve. In this study, we used the CSDA curve for standard rock. Because we have already obtained the path length from the drone-scanning data, the density could be deduced directly by dividing the density-length by the path length. The density map could then be obtained using the above procedure. The density at most regions was approximately 0.5 g/cm 3 . It was assumed that the average density of the volcano is 1.7 g/cm 3 . The measured ratio was approximately one order of magnitude larger than the expected value. Contamination from background events could be the source of the discrepancies between the observed and expected ratios. A background analysis will be performed in the next section. We therefore show only the relative density ( Fig. 6b Fig. 6b ). Background analysis 2 The background could be the source of the discrepancies between the observed and expected muon fluxes. The measured flux was much larger than the value expected from the modified Gaisser formula. To investigate the causes of the higher flux measurement results, a Geant4 [ 25 25 ] simulation package was developed to evaluate the data taking process in the open-sky case. The simulation chain consisted of particle generation, particle transportation, detector response, and digitization. Particle generation was performed using the CRY software [ 26 26 ], which can provide muons with realistic angular and momentum distributions. In addition to muons, background particles, such as electrons and hadrons, can also be generated according to a realistic model. The plastic scintillator detector was constructed and lifted up to 20 as in the actual experiment. The particle interactions with matter, such as ionization and scattering effects, were accounted f or using the EM and hadronic physics list models in Geant4. Monte Carlo truth information, such as the precise positions of the hits, energy deposition in the scintillator bars, and precise track lengths, could be obtained. The simulated hits were digitized according to the pixel size of the detector. The number of hits collected for each detection pixel was obtained from the simulation framework. We did not implement the volcano at present because realistic volcano modelling is complex and time consuming. There were also some backward muons, that is, muons coming from the side without the volcano. Because the volcano is located in the north (azimuth angle 0), the forward muon azimuth angle is assumed to be within the angular range of (-90, 90), and the remaining half of the angular range is defined as the backward direction. Considering the conditions of the practical observation field, the imaging performance can be degraded by backward muons. Based on the Monte Carlo truth, the ratio of muons from the forward direction to muons from the backward direction after selecting the muon samples that passed through all three panels was 100:1.494 under our detector configuration ( Fig. 7 Fig. 7 ). Therefore, the effect of the backward muons was trivial and had little effect on the final imaging results. Using the simulation framework, we studied the charged particle background. The same analysis procedures were applied to the simulated events. The muon-only simulation results are shown as black crosses in Fig. 5a Fig. 5a . The simulated muon-only flux agrees with the prediction results of the modified Gaisser formula. The calculated flux when muon-like particles such as electrons and protons were included in the simulation is represented by the green triangles in Fig. 5a Fig. 5a . The results deviate from the modified Gaisser prediction and overlap with the measured flux data points (red points in Fig. 5a Fig. 5a ). Therefore, our detector received non-negligible amounts of charged particles, and the excess muon flux measured can be explained by the charged particle background. Results and Discussions 1 In this study, we performed the first survey of muon transmission imaging of a volcano in China, and proved its capability to probe the internal density of large geological objects. The first muon radiography system based on plastic scintillator bars was constructed and demonstrated. A muon flux deficit in the elevation angle range covered by the volcano, i.e., 0–20, could be clearly seen for the volcano case. We performed a relative ratio measurement and successfully obtained a profile image of the volcano. The boundary of the volcano is in good agreement with the drone-scanned image ( Fig. 6a Fig. 6a ). We also used the relative density to inspect the inner structures of the volcano ( Fig. 6b Fig. 6b ). The relative density map suggests the existence of significant spatial variation within the volcanic cone. Relatively higher densities were found near the surface of the volcanic cone, whereas relatively lower densities were found near the center of the volcanic cone. This is largely consistent with previous muon radiography results for the Showa-Shinzan lava dome in Japan [ 27 27 ]. Another finding is that the absolute internal density was contaminated by an overwhelming background. One of the major contributions to the background was the mimicking of muon events by charged particles such as electrons and protons ( Fig. 5a Fig. 5a ). Systems with precise timing will be highly useful for suppressing the low-energy electron and proton backgrounds through particle identification based on time-of-flight differences. Radiography detectors with a higher momentum threshold, for example, 0.5 mm-thick aluminum sandwiched between detector layers, can stop the majority of electrons and protons below 1 GeV. In addition, the introduction of a lead plate between the scintillator layers is suggested. The lead plate can stop the passage of low-energy electrons or protons as well as deflect high-energy electrons and protons. Electrons and muons can then be discriminated by comparing their directions before and after the lead plate. These improvements will be implemented in future volcano imaging studies using muon radiography. In summary, muon radiography is a promising technique for imaging internal density structures. Good background control and evaluation strategies are necessary to exploit the potential of muon radiography. References G. Macedonio , M. Martini , Motivations for muon radiography of active volcanoes . Earth, Planets and Space 62 , 139 - 143 ( 2010 ). doi: 10.5047/eps.2009.03.005 http://doi.org/10.5047/eps.2009.03.005 Motivations for muon radiography of active volcanoes A. Barnoud , V. Cayol , V. Niess , et al. , Bayesian joint muographic and gravimetric inversion applied to volcanoes . Geophysical Journal International 218 , 2179 - 2194 ( 2019 ). doi: 10.1093/gji/ggz300 http://doi.org/10.1093/gji/ggz300 Bayesian joint muographic and gravimetric inversion applied to volcanoes 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 . Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 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 K. Borozdin , G. 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 H.K. Tanaka , T. Nakano , S. Takahashi , et al. , High resolution imaging in the inhomogeneous crust with cosmic-ray muon radiography: The density structure below the volcanic crater floor of mt. asama, japan . Earth and Planetary Science Letters 263 , 104 - 113 ( 2007 ). doi: 10.1016/j.epsl.2007.09.001 http://doi.org/10.1016/j.epsl.2007.09.001 High resolution imaging in the inhomogeneous crust with cosmic-ray muon radiography: The density structure below the volcanic crater floor of mt. asama, japan H.K.M. Tanaka , T. Nakano , S. Takahashi , et al. , Imaging the conduit size of the dome with cosmic-ray muons: The structure beneath showa-shinzan lava dome, japan . Geophysical Research Letters 34 , L22311 ( 2007 ). doi: 10.1029/2007GL031389 http://doi.org/10.1029/2007GL031389 Imaging the conduit size of the dome with cosmic-ray muons: The structure beneath showa-shinzan lava dome, japan D. Gibert , F. Beauducel , Y. Déclais , et al. , Muon tomography: Plans for observations in the lesser antilles . Earth, Planets and Space 62 , 153 - 165 ( 2010 ). doi: 10.5047/eps.2009.07.003 http://doi.org/10.5047/eps.2009.07.003 Muon tomography: Plans for observations in the lesser antilles H. Tanaka , T. Nakano , S. Takahashi , et al. , Development of an emulsion imaging system for cosmic-ray muon radiography to explore the internal structure of a volcano, mt. asama . Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 575 , 489 - 497 ( 2007 ). doi: 10.1016/j.nima.2007.02.104 http://doi.org/10.1016/j.nima.2007.02.104 Development of an emulsion imaging system for cosmic-ray muon radiography to explore the internal structure of a volcano, mt. asama N. Lesparre , D. Gibert , J. Marteau , et al. , Geophysical muon imaging: feasibility and limits . Geophysical Journal International 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 X.Y. Pan , Y.F. Zheng , Z. Zeng , et al. , Experimental validation of material discrimination ability of muon scattering tomography at the tumuty facility . Nuclear Science and Techniques 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 X.D. Su , G.L. Zhang , S.P. Xu , et al. , Attenuation coefficients of gamma and x-rays passing through six materials . Nuclear Science and Techniques 31 , 3 ( 2020 ). doi: 10.1007/s41365-019-0717-9 http://doi.org/10.1007/s41365-019-0717-9 Attenuation coefficients of gamma and x-rays passing through six materials M. Guan , M.C. Chu , J. Cao , et al. , A parametrization of the cosmic-ray muon flux at sea-level . https://doi.org/10.48550/arXiv.1509.06176 https://doi.org/10.48550/arXiv.1509.06176 A parametrization of the cosmic-ray muon flux at sea-level T.K. Gaisser , T. Stanev , High-energy cosmic rays . Nuclear Physics A 777 , 98 - 110 ( 2006 ). Special Isseu on Nuclear Astrophysics . doi: 10.1016/j.nuclphysa.2005.01.024 http://doi.org/10.1016/j.nuclphysa.2005.01.024 High-energy cosmic rays D. Chirkin , Fluxes of atmospheric leptons at 600 gev - 60 tev . ( 2004 ). https://doi.org/10.48550/ARXIV.HEP-PH/0407078 https://doi.org/10.48550/ARXIV.HEP-PH/0407078 Fluxes of atmospheric leptons at 600 gev - 60 tev H.A. Bethe , Passage of radiations through matter . Experimental nuclear physics . Experimental nuclear physics P.D. Group , P. Zyla , R. Barnett , et al. , Review of particle physics . Progress of Theoretical and Experimental Physics 2020 , 083C01 ( 2020 ). doi: 10.1093/ptep/ptaa104 http://doi.org/10.1093/ptep/ptaa104 Review of particle physics D.E. GROOM , N.V. MOKHOV , S.I. STRIGANOV , Muon stopping power and range tables 10 mev-100 tev . Atomic Data and Nuclear Data Tables 78 , 183 - 356 ( 2001 ). doi: 10.1006/adnd.2001.0861 http://doi.org/10.1006/adnd.2001.0861 Muon stopping power and range tables 10 mev-100 tev F. Ambrosino , A. Anastasio , A. Bross , et al. , Joint measurement of the atmospheric muon flux through the puy de dôme volcano with plastic scintillators and resistive plate chambers detectors . Journal of Geophysical Research: Solid Earth 120 , 7290 - 7307 ( 2015 ). doi: 10.1002/2015JB011969 http://doi.org/10.1002/2015JB011969 Joint measurement of the atmospheric muon flux through the puy de dôme volcano with plastic scintillators and resistive plate chambers detectors H.K.M. Tanaka , T. Uchida , M. Tanaka , et al. , Cosmic-ray muon imaging of magma in a conduit: Degassing process of satsuma-iwojima volcano, japan . Geophysical Research Letters 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, japan H. Tanaka , T. Nakano , S. Takahashi , et al. , Radiographic imaging below a volcanic crater floor with cosmic-ray muons . American Journal of Science 308 , 843 - 850 ( 2008 ). Radiographic imaging below a volcanic crater floor with cosmic-ray muons R. Han , Q. Yu , Z. Li , et al. , Cosmic muon flux measurement and tunnel overburden structure imaging . Journal of Instrumentation 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 J. Li , Z. Li , R. Han , et al. , Investigation of structures in tunnel overburdens by means of muon radiography . Journal of Instrumentation 17 , P05029 ( 2022 ). doi: 10.1088/1748-0221/17/05/P05029 http://doi.org/10.1088/1748-0221/17/05/P05029 Investigation of structures in tunnel overburdens by means of muon radiography Z. Li , S. Ni , B. Zhang , et al. , Shallow magma chamber under the wudalianchi volcanic field unveiled by seismic imaging with dense array . Geophysical Research Letters 43 , 4954 - 4961 ( 2016 ). doi: 10.1002/2016GL068895 http://doi.org/10.1002/2016GL068895 Shallow magma chamber under the wudalianchi volcanic field unveiled by seismic imaging with dense array H. Gotoh , H. Yagi , Solid angle subtended by a rectangular slit . Nuclear Instruments and Methods 96 , 485 - 486 ( 1971 ). doi: 10.1016/0029-554X(71)90624-0 http://doi.org/10.1016/0029-554X(71)90624-0 Solid angle subtended by a rectangular slit S. Agostinelli , J. Allison , K.a. Amako , et al. , Geant4–a simulation toolkit . Nuclear instruments and methods in physics research section A: Accelerators, Spectrometers, Detectors and Associated Equipment 506 , 250 - 303 ( 2003 ). doi: 10.1016/S0168-9002(03)01368-8 http://doi.org/10.1016/S0168-9002(03)01368-8 Geant4–a simulation toolkit C. Hagmann , D. Lange , D. Wright , in 2007 IEEE Nuclear Science Symposium Conference Record , Cosmic-ray shower generator (cry) for monte carlo transport codes . Vol. 2 , IEEE , 2007 , pp. 1143 - 1146 . doi: 10.1109/NSSMIC.2007.4437209 http://doi.org/10.1109/NSSMIC.2007.4437209 Cosmic-ray shower generator (cry) for monte carlo transport codes R. Nishiyama , Y. Tanaka , S. Okubo , et al. , Integrated processing of muon radiography and gravity anomaly data toward the realization of highresolution 3-d density structural analysis of volcanoes: Case study of showa-shinzan lava dome, usu, japan . Journal of Geophysical Research 119 , 699 - 710 ( 2014 ). doi: 10.1002/2013JB010234 http://doi.org/10.1002/2013JB010234 Integrated processing of muon radiography and gravity anomaly data toward the realization of highresolution 3-d density structural analysis of volcanoes: Case study of showa-shinzan lava dome, usu, japan All authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by Ya-Ping Cheng, Jing-Tai Li and Zhi-Wei Li. The first draft of the manuscript was written by Ya-Ping Cheng and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript. 10.1007/s41365-022-01072-4 This work was supported by the National Natural Science Foundation of China (Nos. 41974064, 42174076 and U1865206) and young Elite Scientists Sponsorship Program by CAST (No. 2019QNRC001) . 2022-01-10 2022-05-22 2022-05-22 2022-07-16 2022-08-27 © 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-01072-4 http://dx.doi.org/10.1007/s41365-022-01072-4 . 2022 Nuclear Science and Techniques 07 33