{'fraction_non_alphanumeric': 0.03432257073751296, 'fraction_numerical': 0.018044304526432988, 'mean_word_length': 4.520983467571004, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 2, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Three-Dimensional Code for Muon Propagation through the Rock: MUSIC

We present a new three-dimensional Monte-Carlo code MUSIC (MUon SImulation Code) for muon propagation through the rock. All processes of muon interaction with matter with high energy loss (including the knock-on electron production) are treated as stochastic processes. The angular deviation and lateral displacement of muons due to multiple scattering, as well as bremsstrahlung, pair production and inelastic scattering are taken into account. The code has been applied to obtain the energy distribution and angular and lateral deviations of single muons at different depths underground. The muon multiplicity distributions obtained with MUSIC and CORSIKA (Extensive Air Shower simulation code) are also presented. We discuss the systematic uncertainties of the results due to different muon bremsstrahlung cross-sections.


Introduction
Muon transport through the rock plays an important role in 'underground physics', in particular in the underground cosmic-ray experiments. The muon 'depth-intensity curve' related to the muon energy spectrum at sea level and to the primary cosmicray specrum, muon multiplicity distribution and 'decoherence curve' related to the mass composition of primaries and to the characteristics of high-energy hadronnucleus interactions are influenced by the muon propagation through the rock. The large underground detectors designed to search for rare events (such as proton decay, neutrino interactions etc.) are not free from the background of cosmic-ray muons and muon-produced secondary particles. The propagation of muons produced by high-energy astrophysical and atmosheric neutrinos in the rock should be also taken into account. These examples show the importance of an accurate three-dimensional simulations of the muon transport in matter.
There are several three-dimensional codes for the simulation of muon propagation through the rock (see, for example, [1,2,3]). They have been used to obtain the characteristics of muon flux and muon bundles deep underground and to describe the experimental data. However, the increase of the accuracy of the recent and future measurements requests an adequate increase of the accuracy of the calculations. The accuracy of the calculations is restricted by the uncertainties in the cross-sections and the assumptions made to simplify the computational procedure (due to the restricted CPU time). Recently, a new corrections to the muon bremsstrahlung and knock-on electron production cross-sections were calculated [4]. A new generation of powerful computers allows to avoid many simplifications in the calculation procedure.
In our code, named MUSIC (MUon SImulation Code) we have used the most recent and accurate cross-sections of the muon interactions with matter and we tried to avoid many simplifications used in the previous simulations. This threedimensional code is based on an algorithm of a one-dimensional muon propagation code described in [5] and used to fit the 'depth-intensity relation' measured by LVD [6,7]. We have treated all processes of muon interaction with high energy transfer (bremsstrahlung, inelastic scattering, pair production and knock-on electron production) as the stochastic ones if the fraction of energy lost by muon is more than 10 −3 . The mean muon path between two interactions with the energy loss more than 10 −3 is about 20 hg/cm 2 . The angular deviation of muons was taken into account not only in the process of multiple scattering (as it was usually done) but also in the processes of bremsstrahlung, inelastic scattering and pair production. We investigated the effect due to the last three processes.
We have performed three-dimensional simulations of the muon transport through Gran Sasso and standard rocks. The results of such simulations for Gran Sasso rock, as well as the complete description of the code and related topics can be found in [8]. Here we present the results of the calculations for standard rock using two different bremsstrahlung cross-sections [4,9]. In this way we investigated the effect due to the uncertainties in the cross-sections. We have also applied the code to calculate muon multiplicity distribution using CORSIKA [10] code for the simulations of Extensive Air Showers.
The main features of the MUSIC code are described in Section 2. The results of the simulation in the standard rock are presented in Section 3. We also discuss in Section 3 the dependence of our results on the cross-sections and assumptions used in the code. The results of the full Monte-Carlo simulation of muon component using CORSIKA and MUSIC can be found in Section 4. In Section 5 we present our conclusions.

The code
In the code we have treated energy losses due to ionization, bremsstrahlung, pair production and muon -nucleus inelastic scattering.
The mean energy loss is usually expressed by where α is the energy loss due to ionization and β is the relative energy loss due to other processes. Ionization can be considered as a continuous process and it is well described by Bethe-Bloch formula. On the contrary, in other processes a particle can lose a big fraction of its energy in a single interaction and it is necessary to treat these kinds of energy losses in two separate parts, one continuous and the other stochastic. To improve the simulation procedure and to study the importance of the stochasticity of the knock-on electron production we have divided the energy loss due to this process into two parts: continuous loss with v < v cut = 10 −3 which we calculate using the Bethe-Bloch formula, and stochastic loss which we treat in the same way as other stochastic processes. Thus, the energy losses due to knock-on electron production, bremsstrahlung, pair production and inelastic scattering are expressed by: where A is the mass number of the material, N is the Avogadro number, v is the fraction of the transferred muon energy and dσ dv is the cross-section of the process. The differential probability that a muon loses a fraction v of its energy per g/cm 2 in a stochastic way is given by It is quite difficult to choose a suitable value for v cut , because a little value of v cut , that means an accurate simulation, results in an increase of CPU's time. We have chosen v cut = 10 −3 that seems meet all requirements. On a DEC AlphaStation 250 4/266 the simulation of 10 TeV muons propagated through 3 km.w.e. takes 0.02 s/muon.
We note that most of the other codes use higher value of v < v cut and/or are considering the knock-on electron production as a continuous process.
The aim of the code is to propagate a muon through a layer X of standard rock (we used A=22, Z=11, ρ= 2.65 g/cm 3 ). If we know the initial energy of the particle E 0 , its starting point x 0 , y 0 , z 0 in the rock and its direction θ 0 , φ 0 , we can simulate with an iterative procedure, step by step, muon interaction point according where − → r is the vector which defines the interaction point in the space and λ is the mean free path. Then the continuous and stochastic energy losses, as well as lateral and angular displacement due to multiple Coulomb scattering are evaluated. Moreover the code simulates the angular deviations caused by the stochastic processes.
Finally we check if the muon has crossed the layer X, and in this case we store its final energy E f , its direction θ f , φ f and the coordinates of its arrival point To describe pair production we have used the cross-section given by Kokoulin-Petrukhin [11]. The code samples separately muon energy loss and angular deflection. The second one is sampled following the parametrization proposed in [12] (see Appendix A).
Two different muon bremsstrahlung cross-sections, one given by Bezrukov and Bugaev [9] and the other recently proposed by Kelner, Kokoulin and Petrukhin [4] are available in MUSIC. These two calculations differ mainly by the account of bremsstrahlung on atomic electrons, when the photon is emitted by the target electron.
The code samples the fraction of muon energy transferred to the photon and separately the muon scattering angle, following also in this case the parametrisation proposed by [12] (see Appendix A).
We carried out additional test of the simulation of muon deflection angle due to pair production and bremsstrahlung process using a simplified procedure based on the GEANT parametrizations [13,14]. We concluded that this approach leads to a little underestimation of the muon deflection due to these processes with respect to Van Ginneken's parametrizations.
Among stochastic processes, the angular deflection due to inelastic scattering on nuclei is expected to be dominant [15]. In this case to calculate the muon energy loss and angular deflection we have used double-differential cross section d 2 σ dvdQ 2 obtained by Bezrukov and Bugaev [16], where Q is the four-momentum transferred. We have sampled the muon scattering angle θ µ according to the aforementioned cross-section using the relation Q 2 = 2(E µ0 E µf − P µ0 P µf cosθ µ − m 2 µ ). Since the muon inelastic scattering is expected to be the most important for angular deviation of muons (with the exception of multiple scattering), this process has been treated in a more correct way (using the double-differential cross-section) in comparison with other stochastic processes. However, the cross-section of this process proposed in [16] is valid only at low Q 2 (the accuracy is not worse than 11% for Q 2 ≤ 3 GeV 2 and energy transfers from 1 to 10 6 GeV [16]).
Some authors [12,17,18], pointed out the complexity of the problem to sample simultaneously the emission angle of the photon, the deflection angle of the muon and the fractional energy loss v. The accurate and simultaneous simulation of the angles and energies can be important for muons exceeding 1 TeV, where the deflection due to bremsstrahlung could be dominant [12]. As it will be discussed in the next section, however, even if we used the same approach as in [12], we didn't find such behaviour, and the multiple scattering appears to be the dominant mechanism of deflection of muons passing through the large thickness of matter. A complete treatment of these aspects is currently beyond the aim of our code.
We have treated multiple Coulomb scattering in the Gaussian approximation [20]. We have considered the projections of muon angular deviation, θ µ , and lateral displacement, l, on two orthogonal planes which include the initial direction of muon. The projections of angular deviation and lateral displacement of muon with respect to its initial direction at each step between two interactions with v > v cut are distributed according to correlated Gaussians with mean values and dispersions calculated using the formulae: where ǫ ms = 21.2 MeV, λ = 26.48 g/cm 2 is the radiation length in standard rock and ρ = 2.65 g/cm 3 is the standard rock density. As discussed in the next section, we have checked also if the use of a more accurate treatment of multiple scattering (e.g. Molière theory) leads to different results.

Results
We have propagated 100000 muons through 3 km.w.e. and 1000000 through 10 km. w.e. with MUSIC using Kelner-Kokoulin-Petrukhin's cross-section for bremsstrahlung. Initial muon energy has been sampled according to a primary energy spectrum with index γ = 3.7. In the case of 3 km.w.e. minimal muon energy E 0 was fixed to 900 GeV, while for 10 km.w.e. E 0 was 9 TeV, corresponding in each case to a survival probability of less or about 0.3%. Figures 1, 2 and 3 show energy, lateral and angular distributions at 3 and 10 km.w.e. The mean muon energy is 250 GeV at 3 km.w.e. and 367 GeV at 10 km.w.e. This result is in agreement with the simple consideration that the mean muon energy at a depth h increases with h.
Mean angular deviation < α > is 0.56 • at 3 km.w.e. and 0.45 • at 10 km.w.e. Angular deviation is mainly caused by multiple Coulomb scattering. This process produces a deviation that increases when energy decreases. This explains why at 3 km.w.e., when < E h > is lower, < α > is higher than at 10 km.w.e. Figures 4 and 5 show lateral and angular distributions obtained with MUSIC and PROPMU -a code for muon propagation in rock written by P. Lipari and T. Stanev [2]. The codes are completely independent, but their results are practically the same. This fact confirms the correctness of both simulations and shows that the effects of angular deviations due to stochastic processes, with the approximations discussed in previous section, are really very little, 1-2 orders of magnitude less than the effects due to multiple scattering. In fact, as Lipari and Stanev treat only angular deviation due to multiple scattering, and in MUSIC all possible angular deviations are taken into account, the distributions obtained with these two codes do not show any appreciable difference.
In Figure 6 we show the contribution of each process to the angular deviation of muons which pass 3 km.w.e. To do this plot we have repeated five times the same simulation leaving active one mechanism of deflection each time. For multiple scattering we have used either Gaussian treatment described in previous section or Molière theory subroutines provided by GEANT [19]. A more detailed treatment (Molière theory) results in the similar distribution compared with that of Gaussian treatment, with a small increase of mean deflection angle. However, the original Molière theory (see [19] and references therein) does not provide the lateral displacement, more important from the experimental point of view. As can be clearly seen from Figure 6 the multiple scattering dominates over other processes. This result remains true for 10 km.w.e. and for a monochromatic beam of 10 TeV muons.
To calculate survival probabilities, P (E 0 , h), and muon energy distributions underground, P (E 0 , h, E), we have propagated muons with different initial energy (100000 muons for each initial energy) through the standard rock. P (E 0 , h, E) means the probability that muon with initial energy E 0 will reach the depth h with energy E. Figure 7 shows the survival propabilities obtained with MUSIC for E 0 from 1 to 10 4 TeV. Figure 8 shows P (E 0 , h) for h = 4 km.w.e. and h = 10 km.w.e., calculated with MUSIC and with PROPMU. The curves are very similar also in this case. Our values are on average (1-2)% lower: this difference can be due to the new cross-section for bremsstrahlung used in MUSIC.
We have studied the effects of the cross-section difference and of the stochasticity of different processes on the characteristics of the muon flux underground. Using the sea-level vertical muon energy spectrum, proposed by T.Gaisser [21] with the spectral index of primary spectrum γ = 2.70, and the distributions P (E 0 , h) and P (E 0 , h, E) we have calculated the vertical muon intensities I(h) and mean muon energies < E(h) > at different depths underground: where dI(E 0 ) dE 0 is the vertical muon energy spectrum at sea level and dI(E,h) dE is the vertical muon spectrum at the depth h. The latter one is obtained using the formula: We have calculated the muon intensities and mean muon energies for two different muon bremsstrahlung cross-sections (from [4] and [9]) and for two values of v cut (10 −3 and 10 −2 ). We have carried out additional tests considering the knock-on electron production and pair production as continuous processes. The results of the calculations are presented in Table 1 for I(h) and Table 2 for < E(h) >. The first column in both tables shows the depth in km.w.e. of standard rock. Other columns show the intensity in units cm −2 s −1 sr −1 (Table 1) and the mean muon energy in GeV (Table 2). 2nd column present our basic values calculated with the muon bremsstrahlung cross-section from [4] and v cut = 10 −3 . The results obtained with bremsstrahlung cross-section from [9] and v cut = 10 −3 are presented in column 3. The values in column 4 were calculated with the cross-section from [4], v cut = 10 −3 and with continuous energy loss due to knock-on electron production. The simulation related to column 5 has been carried out under the same assumptions as that for column 4 but both knock-on electron and pair production have been treated as continuous processes. Finally, column 6 shows the results of the simulation similar to that for column 4, but with v cut = 10 −2 . The accuracy of the values presented in Tables 1 and 2 is of the order of 1%.
Let us first compare the intensities obtained with different assumptions used in the code. The intensities calculated with the muon bremsstrahlung cross-section from [9] are systematically higher (by (4-8)%) than those with the cross-section from [4]. The relative difference slightly increases with depth. If the depth-intensity relation measured in some experiment is used to obtain the power index of the muon energy spectrum at sea level (or the primary spectrum), then, the use of the crosssection from [9] will result in a higher power index (softer muon spectrum) than that with cross-section [4]. The difference can be of the order of 0.01.
The treatment of the knock-on electron production and pair production as continuous processes results in the decrease of muon intensities (up to 10% at high depths). The relative difference also increases with depth. The power index of the primary spectrum which can be evaluated from the underground measurements will be lower in this case. The difference in the power index is negligible if only knock-on electron production is treated as a continuous process but can be of the order of 0.02 if both knock-on electron and pair production are not treated stochastically.
The increase of v cut from 10 −3 up to 10 −2 results in a small decrease of the muon intensities ((2 − 3)%). However, the increase of the relative difference with depth is small and may not influence much the power index which can be evaluated from the measured depth -intensity relation. The increase of power index in this case should not exceed 0.01.
The mean muon energies calculated with our basic assumptions (2nd column of Table 2) are higher than those obtained with another cross-section and/or other assumptions. The difference is not negligible in the cases of the cross-section from [9] and/or the continuity of the pair production process.
We have checked also the influence of the cross-section difference and v cut on the angular and lateral distributions. The initial sea-level muon spectrum in the power-law form with power index 3.7 has been used. With the cross-section from [4] we have obtained the mean angular deviation 0.56 • and 0.45 • for 3 and 10 km.w.e., respectively, while with the cross-section from [9] the values are 0.56 • and 0.46 • . Statistical uncertainty is of the order of 0.01 • . The mean lateral displacement is 1.56 m and 1.87 m with the cross-section from [4] for 3 and 10 km.w.e., respectively, and 1.59 m and 1.93 m with the cross-section from [9]. Statistical uncertainty is about 0.01 m. We can conclude that the difference in the cross-sections does not influence much the angular and lateral distributions at any depth. A small but statistically significant difference of the mean lateral displacements at large depth can be attributed to the difference of the mean muon energies at this depth (see Table 2 and Eq. (5)).
We have carried out also the propagation of muons with v cut = 10 −2 . Actually, the increase of v cut from 10 −3 to 10 −2 results in an increase of the average step between two consecutive interactions by one order of magnitude. Thus, this procedure allows us to test the dependence of the resulting angular deviation and lateral displacement on the step between two consecutive interactions. We didn't find any significant difference between the results with v cut = 10 −3 and v cut = 10 −2 . This proves the correctness of the three-dimensional simulation procedure.

Muon multiplicity
We have roughly investigated also the effect of bremsstrahlung cross-section on muon multiplicity distribution. We have used the code CORSIKA [10] to simulate the development of a shower produced by a primary cosmic-ray nucleus in atmosphere. We have sampled two sets of 10000 showers produced by protons and iron nuclei of 10 4 TeV to have a good statistics of muons.
To transport shower muon component through 3000 m.w.e. of standard rock with MUSIC, we have used two cross sections for bremsstrahlung from [9] and [4]. For the case of infinite underground detector, we have produced two muon multiplicity distributions obtained with different cross-sections. Figures 9 and 10 show these distributions for showers produced by protons and iron nuclei. In Table 3 the mean values < N µ > are presented, their percentage difference Cross − section < N µ > p < N µ > F e [4] 7.89 ± 0.04 20.96 ± 0.07 [9] 8.20 ± 0.04 21.71 ± 0.07 In order to evaluate the possibility that the differences between distributions could be caused by statistical fluctuactions, we have done a χ 2 test: the probability that muon multiplicity distributions obtained with different cross sections come from a single distribution is ∼ 10%.
The effect of different bremsstrahlung cross sections on muon multiplicity can be underlined using a fit with a negative binomial [22]. It is evident, as shown in Fig.  11, that fit to muon multiplicity distribution, obtained with cross-section from [9], is practically a translation of the other fit due to different values of < N µ >.
These simple considerations proves that the muon interaction cross-sections can influence the muon multiplicity distributions underground and, hence, the primary cosmic-ray composition which can be evaluated from the underground measurements.
In real experiments, however, the observed multiplicity distribution is driven by the primary energy spectrum and the assumed primary composition and a such complete simulation to investigate the differences due to the muon cross sections is beyond the aim of this paper. We can state only that the effect is quite small (but not negligible) and dominated by other systematic uncertainties.

Conclusions
We have presented the new three-dimensional Monte-Carlo code MUSIC for muon propagation through the rock. The main advantages of this code are: i) the treat-ment of all main processes of muon interaction with matter (bremsstrahlung, inelastic scattering, pair production and knock-on electron production) with the fractional energy loss v > 10 −3 as stochastic processes; ii) the use of the most accurate cross-sections of muon interactions with matter; iii) the accounting of the angular deviation of muons due to stochastic processes. We can conclude that the angular and lateral distributions of muons at any depth are fully determined by multiple Coulomb scattering, while other processes give a contribution which is 1-2 orders of magnitudes less. However, the very rare processes of muon scattering to large angles (for example, backward scattering) which can be due to inelastic scattering and/or bremsstrahlung need more accurate treatment of these processes at large Q 2 . We have studied the dependence of the characteristics of muon flux underground on the bremsstrahlung cross-section and assumptions used in the code. We can conclude that the treatment of pair production as a stochastic process is quite important for the characteristics of muon flux underground. The separation of stochastic and continuous energy losses at v cut = 10 −3 is optimal from the point of view of CPU time and the accuracy of the results. The accurate knowledge of the cross-sections is needed to calculate the intensity of single muons and muon bundles underground. Even a small (1-2)% increase of the muon bremsstrahlung cross-section [4] comparing with that from [9] results in a significant (6%) decrease of the single muon intensity at 3 km.w.e. and in a shift (however, quite small) of muon multiplicity distributions toward low multiplicities. This could produce some bias in the interpretation of the muon 'depth-intensity relation' and multiplicity distribution measured in an underground experiment when they are analysed in terms of the primary spectrum and composition.
, and k 5 is related to k 4 by continuity of the function.
The mean angle between incident and scattered muon due to pair production process has been parametrised as [12]: where m e is the electron mass.
For a given value of v sampled according to the cross-sections from [4,11] without any additional approximation the mean angle < θ 2 > 1/2 between incident and scattered muon has been calculated using the Eqs. (9,10). The spatial angle θ has been chosen from an exponential in the variable θ 2 with the mean value < θ 2 >.
We note again that we did not use any additional parametrisation for the muon inelastic scattering cross-section which has been taken from [16]. For a given value of v sampled according to the cross-section integrated over Q 2 , the muon scattering angle θ has been simulated using the double-differential cross-section.
For all stochastic processes we have sampled randomly the angle φ in the plane perpendicular to the muon direction.
The code is available upon the request by e-mail to: antonioli@bo.infn.it or kudryavtsev@vaxmw.tower.ras.ru        : Comparison of muon multiplicity distributions at 3 km.w.e. obtained using bremsstrahlung cross sections from [9] (open circles) and [4] (triangles). Muons are produced in 10000 proton showers of energy E o = 10 4 TeV Figure 10: Comparison of muon multiplicity distributions at 3 km.w.e obtained using bremsstrahlung cross sections from [9] (open circles) and [4] (triangles). Muons are produced in 10000 showers initiated by iron nuclei of energy E o = 10 4 TeV Figure 11: Comparison of the fits to the multiplicity distributions with negative binomial functions. Multiplicity distributions have been obtained using muon bremsstrahlung cross sections from [9] (dashed line) and [4] (solid line) for proton and iron primaries.