Principle

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].

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]. 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]. Compared to Gaisser's formula [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: $\begin{array}{c}\Phi (\theta \mathrm{,}{E}_{\mu })=0.14{\left[\frac{E_{\mu }}{\text{GeV}}\left(1+\frac{3.64\text{GeV}}{E_{\mu }{(\cos {\theta }^{*})}^{1.29}}\right)\right]}^{-2.7}\\ \left(\frac{1}{1+\frac{1.1E\cos {\theta }^{*}}{115\text{GeV}}}+\frac{0.054}{1+\frac{1.1E\cos {\theta }^{*}}{850\text{GeV}}}\right)\mathrm{,}\end{array}$ (1) where Eμ is the muon energy. The quantity cosθ^* is related to the zenith angle θ by: $\cos {\theta }^{*}=\left(\frac{{{\displaystyle \cos }}^2\theta +p_1^2+p_2{(\cos \theta )}^{p_3}+p_4\cos {\theta }^{p_5}}{1+p_1^2+p_2+p_4}\right)\mathrm{,}$ where p₁, p₂, p₃, p₄, and p₅ are the parameters given in Table 1 in [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: $N_{\text{opensky}}(\theta )={\displaystyle {\int }_{0.1}^{\infty }}\Phi (\theta \mathrm{,}E)S_{\text{eff}}(\theta )ϵ\Delta T\text{d}E\mathrm{,}$ (2) 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 $N_{\text{volcano}}(\theta )={\displaystyle {\int }_{E_{\text{min}}}^{\infty }}\Phi (\theta \mathrm{,}E)S_{\text{eff}}(\theta )ϵ\Delta T\text{d}E\mathrm{,}$ (3) 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] or the numerical values provided by the Particle Data Group [16]. The continuous slowing down ability (CSDA) table by [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). 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.

Fig. 1Full-screen View

(Color online) The muon survival ratios (percentage of muons with energy larger than the energy along the x axis) at different zenith angles. The calculation was based on the modified Gaisser formula and the CSDA table from [17]. The energy spectrum for horizontal muons (purple line) is much harder than that for vertical muons (green line). The density-length and minimum muon energy required are also shown (red solid line, left y axis).

Muon Detector

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], are widely used to detect high-energy charged particles, particularly muons. Typical RPCs have excellent detection 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]) have high angular and spatial resolutions. In addition, they are portable and do not require electric power.

The plastic scintillator detector ([19, 20, 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 Nx = Ny = 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² 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, 22].

Detector placement

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].

Before the measurements, we conducted a geological survey using an UAV. A high-precision 3D image of the volcano was obtained (Fig. 2a). The horizontal resolution of the drone image is 1 cm and the vertical resolution 10 cm.

Fig. 2Full-screen View

(Color online) (a): 3D image of volcano collected by the unmanned aerial vehicle system. (b): Muon path length inside the volcano from the detector point of view. The black solid lines are 100 m contour lines superimposed on the path length map. The black and red dashed regions denote the fields of view of the detector facing the volcano during the two stages of data taking (10 west of north and 30 east of north).

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. 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

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). The detector was then rotated to the true west to record the muon flux from the open sky (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) 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. Data were collected over a period of approximately one day. All the muon events were recorded for later analysis.

Fig. 3Full-screen View

(Color online) Detector arrangement in the 3 stages of data taking.

Data Processing Method

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: $\begin{array}{l}\left|(x_1-x_2)-(x_2-x_3)\right|<1\\ \left|(y_1-y_2)-(y_2-y_3)\right|<1\end{array}$ (4) where 1 on the right-most side denotes the ID difference between the scintillation bars, and (x₁, y₁), (x₂, y₂), and (x₃, y₃) 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. 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.

Fig. 4Full-screen View

(Color online) Detector response-related corrections. (a): Overall efficiency of the 96 scintillator bars. There are 6 panels with 16 scintillation bars on each panel. (b) and (c): Solid angle and acceptance as functions of θx and θy.

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] and shown in Fig. 4b.

In 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² for each detection pixel, forming an acceptance of approximately 16 cm²· sr. The acceptance function (cm²· sr) corresponding to the detector configuration is shown in 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. The measured results as a function of the zenith angle are shown in 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 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 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. This proves that the muon detector, including its solid angle and acceptance, was understood correctly.

Fig. 5Full-screen View

(Color online) Flux measurement results. (a): Muon flux measurement results in the open-sky case (red dots), prediction by modified Gaisser formula (solid blue line), and simulated flux results under CRY input w/ and w/o e^-, p backgrounds (green and black dots, respectively).(b) and (c): Pixel-by-pixel muon flux measurement results in the open-sky and volcano cases. (d): Pixel-by-pixel comparison of the open-sky flux measurement results and the modified Gassier formula prediction.

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. $\begin{array}{c}{K}^{\alpha }=\frac{N_{\text{V}}(\theta )/\Delta T_{\text{V}}}{N_{\text{O}}(\theta )/\Delta T_{\text{O}}}\\ =\frac{{\displaystyle {\int }_{E_{\text{min}}}^{\infty }}\Phi (\theta \mathrm{,}E)S_{\text{eff}}(\theta )ϵ\text{d}E}{{\displaystyle {\int }_{0.1}^{\infty }}\Phi (\theta \mathrm{,}E)S_{\text{eff}}(\theta )ϵ\text{d}E}\end{array}$ (5) 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], 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 and used for the volcano density analysis.

Muon radiography results

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. 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. In this way, the profile image of the Laoheishan volcano was obtained using the muon telescope.

Fig. 6Full-screen View

(Color online) Muon radiography results. (a): Muon flux ratio plot of the volcano and open-sky cases. There is an obvious separation between air, which is represented by the blue part, and the volcano, which is represented by the colorful part. (b) Unfolded density plot of the volcano. Because of noise, the absolute density does not have any physical meaning. Here, we show the relative density map.

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³. It was assumed that the average density of the volcano is 1.7 g/cm³. 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).

Background analysis

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] 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], 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 for 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). Therefore, the effect of the backward muons was trivial and had little effect on the final imaging results.

Fig. 7Full-screen View

(Color online) θ, ϕ distribution of muons that can pass through all three detection panels.

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. 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. The results deviate from the modified Gaisser prediction and overlap with the measured flux data points (red points in 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.