Simulation of Deflection Uncertainties on Directional Reconstructions of Muons Using PROPOSAL

Large scale neutrino detectors and muon tomography rely on the muon direction in the detector to infer the muon's or parent neutrino's origin. However, muons accumulate deflections along their propagation path prior to entering the detector, which may need to be accounted for as an additional source of uncertainty. In this paper, the deflection of muons is studied with the simulation tool PROPOSAL, which accounts for multiple scattering and deflection on stochastic interactions. Deflections along individual interactions depend on the muon energy and the interaction type, and can reach up to the order of degrees -- even at TeV to PeV energies. The accumulated deflection angle can be parametrized in dependence of the final muon energy, independent of the initial muon energy. The median accumulated deflection of a propagated muon with a final energy of 500 GeV is $\theta_{\text{acc}} = 0.10{\deg}$ with a 99 % central interval of $[0.01{\deg}, 0.39{\deg}]$. This is on the order of magnitude of the directional resolution of present neutrino detectors. Furthermore, comparisons with the simulation tools MUSIC and Geant4 as well as two different muon deflection measurements are performed.


Introduction
The directional reconstruction of muons is an essential task for muography or large scale neutrino detectors such as IceCube [1] or KM3NeT [2]. In both cases, the muon direction is measured at its crossing through the instrumented volume, which is then utilized to infer its origin or the origin of the parent neutrino. However, muons may propagate many kilometers prior to entering the detector while interacting with the surrounding medium. Along their propagation, muons can undergo many of thousands of interactions, depending on their energy and propagation distance. These interactions can lead to a deflection of the muon that may need to be accounted for as an additional source of uncertainty in these measurements. Current angular resolutions are above 0.1°for TeV to PeV energies in IceCube [3] and below 0.2°for energies greater than 10 TeV in KM3NeT/ARCA (part of KM3NeT dedicated to search for very high-energetic neutrinos) [4].
To study the impact of the muon deflection on the angular resolution of current neutrino detectors, the paper is structured as follows: in section 2, the lepton propagator PROPOSAL is briefly described. In section 3, PROPOSAL [5,6] is used to study the muon deflection per interaction. The accumulated deflection is analyzed and compared to a e-mail: pascal.gutjahr@udo.edu (corresponding author) the propagation codes MUSIC [7,8] and Geant4 [9,10] and data from two experiments in section 4. The findings of this study are summarized in section 5.

Overview of the Simulation Tool PROPOSAL
The tool PROPOSAL [5,6] propagates charged leptons and photons through media and is used in this paper to simulate the deflection of muons. All relevant muon interaction types as bremsstrahlung [11,12], photonuclear interaction [13] with shadowing [14], electron pair production [15] with corrections for the interaction with atomic electrons [16], ionization described by the Bethe-Bloch formula with corrections for muons [17], and the decay are provided by PROPOSAL. The interaction processes are sampled by their cross section. Since energy losses with the massless photon as secondary particle can be arbitrarily small, an energy cut is introduced to avoid an infinite number of bremsstrahlung interactions and furthermore to increase the runtime performance. The cut is applied with a minimum energy loss using two parameters -a total and a relative energy cut denoted as e cut and v cut with the energy E of the par-arXiv:2208.11902v2 [astro-ph.IM] 2 Feb 2023 ticle directly before the interaction. By the introduction of this energy cut, the next significant energy loss with E loss ≥ E loss,min is treated as a stochastic energy loss in the propagation. All energy losses with E loss < E loss,min between two stochastic losses are accumulated and lost continuously, denoted as continuous energy loss. The methodical uncertainties are small for a relative energy cut v cut 1, however, using a small energy cut increases the runtime. Typically, a relative energy cut of v cut ≤ 0.05 is chosen which enables accurate propagations at low runtimes. The total energy cut depends on the minimum visible energy loss in the detector. It is often set to e cut = 500 MeV. The propagation process is defined by an initial muon energy E i and two stopping criteria -a final energy E f, min and a maximum propagation distance d max . If the last interaction of a propagation is sampled by a stochastic interaction, the true final energy E f can become lower. Since muons are unstable, a decay leads to a premature stop, which is negligible for high energies.
In PROPOSAL, the deflections for stochastic interactions are parametrized by Van Ginneken in Ref. [18] with a direct calculation of the deflection in ionization using four-momentum conservation. Furthermore, there are parametrizations for stochastic deflections given in Geant4 [9,10] for bremsstrahlung and photonuclear interaction, which are also available. To estimate the deflection along a continuous energy loss, multiple scattering described by Molière (MSM) [19] and the Gaussian approximation by Highland (MSH) [20] can be chosen. MSM results as a summation of elastic scatterings of one particle at another particle, called single scattering. Thus, the muon is deflected by a single angle for each continuous loss, analogous to a stochastic loss. The orientation of the deflection in the plane perpendicular to the muon direction is sampled uniformly between 0 and 2π. The latest updates with a detailed description of the whole tool can be found in Ref. [21]. The stochastic deflections have been implemented in PROPOSAL recently and they are described and studied in Ref. [22]. A publication describing the updates in PROPOSAL is in preparation [23]. All simulations are done with PROPOSAL 7.3.1.

Muon Deflection per Interaction
First, the stochastic deflections described by Van Ginneken [18] and implemented in PROPOSAL are investigated in combination with the two multiple scattering methods. For this purpose, 1000 muons are propagated from E i = 1 PeV to E f, min = 1 TeV. The deflections per interaction are presented for each interaction type and the sum over all types in Figure 1. The size of individual deflections extend over several orders of magnitude with a median of 3.9 × 10 −6 deg and a 95 % central interval of [

Accumulated Muon Deflection
As shown in Section 3, the deflection per interaction is lower than ∼ 1°in general. Since these deflections accumulate along the propagation path, the angle between the incoming and the outgoing muon direction is analyzed. This angle limits the angular resolution for neutrino source searches utilizing incoming muons, since there is no information about the muon before the detector entry.

Comparison with MUSIC and Geant4
First, the deflections in PROPOSAL are compared to the tools MUSIC and Geant4. MUSIC is a tool to simulate the propagation of muons through media like rock and water considering the same energy losses as in PROPOSAL. Also, the losses are divided into continuous and stochastic energy losses by a relative energy cut. Several cross sections, multiple scattering methods, and parametrizations for stochastic deflections are available. For these studies, the same cross section parametrizations as in PROPOSAL are chosen, except those for photonuclear interaction [24,25,26]. The stochastic deflections are also parametrized by Van Ginneken [18]. The Gaussian approximation [20] is set as multiple scattering. Geant4 is another common toolkit to simulate the passage of particles through matter using the same cross section parametrizations except for photonuclear interaction [27]. The simulation is very precise and especially made for simulations in particle detectors [9,10].
A comparison of all three tools is shown in Figure 2 for the accumulated deflection angle θ acc and the lateral displacement x. The propagation of 1 000 000 muons with an initial energy of E i = 2 TeV is simulated in water with a maximum propagation distance of 3 km. Four different settings are studied in PROPOSAL to compare the results with the two multiple scattering methods and the different stochastic deflection parametrizations. The deflection angles are similar in all cases. The largest displacements are exhibited by Geant4 and PROPOSAL with Molière scattering, which leads to the largest deflections and thus to a larger displacement. PROPOSAL with Highland scattering and MUSIC have less outliers, since large deflections are neglected in the Gaussian approximation [20]. The combination of Highland and Van Ginneken's photonuclear interaction parametrization leads to the smallest displacement. This is due to the fact that the angle is sampled from the root mean squared angle in the exponential distribution in the parameterization for photonuclear interaction by Van Ginneken, which neglects outliers to larger angles. In general, the lateral displacements differ, although the angles are very similar in all simulations. This can be  Table 1. Table 1. The medians of deflections θ per interaction from Figure 1 are presented for each stochastic interaction type, the two multiple scattering methods, and the total distribution including MSM with the upper and lower limits of the 95 % central intervals. The largest median deflection is caused by photonuclear interaction.
explained by the location of the deflection. If larger deflections occur sooner, they lead to further displacements during propagation, although the angle remains the same. Detailed information are given in Table 2. The largest average deflections are obtained in Geant4 with θ = 0.27°a nd x = 3.3 m, while MUSIC provides the lowest ones with θ = 0.22°and x = 2.6 m. The results of PROPOSAL lie between these two tools. Hence, the mean values of all tools are very close to each other and therefore consistent.

Data-MC Agreements
In the following, two comparisons are performed with measured data for different energies and media. A measurement of muon deflections in low-Z materials was done by Attwood et al. [28]. From this it can be seen that for Z < 4 the scattering angle is overestimated by Molière scattering in Geant4. Hence, the lower scattering in PROPOSAL leads to a better agreement especially in the region of outliers. The comparison is done in liquid H 2 with a thickness of 109 mm and an initial particle energy of E i = 199 MeV. This energy is obtained via the energy-momentum relation of a beam momentum of p = 168.9 MeV/c used in Ref. [28]. In PROPOSAL, the simulations are done with two different energy cuts v cut = 10 −3 and v cut = 10 −5 , but there is no significant difference between the resulting deflections. Even though in the logarithmic figure the simulation data agree well with the measured data, it is clear from the data-MC ratio that the deviations are up to 200 % in some cases. Thus, the deflections are described correctly only  Table 2.
The results for MUSIC and Geant4 are taken from Ref. [8]. Table 2. The survival probability ps defined by the ratio of all muons that reach the propagation distance of 3 km and the amount of muons stopping before due to large energy losses and muon decays, the mean survived muon energy E f , the mean scattered angle θ, and the mean displacement x are presented for all cases from Figure 2. For all means, the standard deviation is given. The largest deflection and displacement is observed in the tool Geant4, which has the lowest mean survived energy. The lower the energy, the larger the deflection. in a first approximation. The comparison is presented in Figure 3. The second measurement of muon deflections is done for higher energetic muons of p = 7.3 GeV/c by Akimenko et al. [29]. In total, 31 125 muons are propagated through a 1.44 cm thick copper layer. Again, the two energy cuts mentioned before and the effect of stochastic deflections in comparison with Molière scattering only are checked. Neither between the two energy cuts, nor when using the stochastic deflection a significant difference occurs. PROPOSAL simulates more large deflections than observed in these data. This observation differs from the comparison with Attwood, in which less higher deflections are simulated. In general, the higher muon energy leads to smaller deflections. Also simulations with Geant4 v11.0.3 using the default settings and the PhysicsLists QBBC and FTFP_BERT are performed. At angles larger than 0.2°, Geant4 simulates less large deflections than PROPOSAL and less than expected in the data. Similar to the comparison with Attwood, data-MC mismatches larger than 100 % are observed. There are no differences between the two PhysicsLists in the resulting deflections. The result is presented in Figure 4.

Muon Deflection Impact on Angular Resolutions
For neutrino source searches based on muons entering the detector, it is important to study the impact of the muon deflection on the angular resolution to estimate whether or not this needs to be taken into account as an additional source of uncertainty. For this purpose, four different initial energies from E i = 10 TeV to E i = 10 PeV are used and the final energy is set to E f, min ≥ 1 GeV with E f, min < E i for each simulation. This energy range covers the muon energies typically measured in neutrino experiments. In total 44 simulations are performed. To compare the results of these simulations, the medians of the deflection distributions with a 99 % central interval are presented in Figure 5. The lower the final muon energy, the larger the accumulated deflection. For energies E f = 1 PeV, the median deflection is 10 −4 deg. For energies E f = 100 GeV, angles larger than 1°are possible. For energies E f ≤ 1 TeV, there is a small overlap of the deflection with the angular resolution of KM3NeT/ARCA [4,30]. At low energies of E f = 5 GeV, the upper limit of the deflections affects the resolution of SuperKamiokande [31]. The kinematic scattering angle between the incident neutrino and the produced muon is larger than the deflection in the presented region from 60 GeV to 200 TeV. For energies below 2 TeV, the muon deflections are of the same order of magnitude as the kinematic angle and thus become increasingly relevant. Here it must be noted that the kinematic angle and the resolution of ARCA in Ref. [4] as well as the resolutions of ORCA (part of KM3NeT optimized to study atmospheric neutrinos in the GeV energy range) [32] and ANTARES [33] are presented in dependence of the neutrino energy. Hence, a rescaling to the muon energy is applied using the average energy transfer of the neutrino to the nucleus [34]. This shifts the curves to lower energies. Since all of these simulations are done in ice, the same simulations are done in water to compare the results for water-based experiments. The deviations of the medians are less than 1 % for all energies and therefore not shown. The accumulated deflections and also the propagated distances of Figure 5 are presented in Table 3. Muons are able to propagate various distances for a fixed final muon energy depending on the stochasticity of the energy losses.
Note that the distribution of deflection angles at a given final energy E f in Figure 5 overlap for differing initial energies. This result indicates that the total deflection of a muon primarily depends on the final muon energy. The initial muon energy is nearly irrelevant. Hence, the reconstructed muon energy in a detector can be used to estimate the deflection. For this purpose, a polynomial of degree three as can be used with the parameters a = + 0.0176 ± 0.0018 , c = + 0.0929 ± 0.0527 , b = − 0.2328 ± 0.0185 , d = + 0.0726 ± 0.0404 , in the logarithmic space via In general, the function f (x) in Eq. (2) describes the median deflection of a muon after a propagated distance in ice for a given, respectively measured energy to estimate the deflection before the detector entry. This equation is valid for muon energies between 1 GeV and 50 PeV. Basically, it should be mentioned here that the data-MC comparisons shown earlier are for energies of E i = 199 MeV and E i = 7.3 GeV, which are much lower than the energies in these simulations.
To analyze the impact of the propagation distance on the muon deflection, another simulation is done. This time, the initial energies are not fixed and sampled from an atmospheric muon flux at sea level by [36] with a weighting of E 3.7 for energies from 10 GeV to 10 10 GeV. The final energies are sampled similar. The resulting deflections are presented in Figure 6. From this follows, that the median deflection is not impacted by the propagation distance, if the initial muon energy is unknown. This is a realistic scenario for example for a neutrino telescope, since the only known value is the reconstructed muon energy at the detector entry. There are no information about the initial energy and the propagation distance. Finally, the muon deflection can be estimated only by the reconstructed muon energy.

Relevance for Muography
Muography is a technique to study the inside of structures with a wide field of applications such as the monitoring of volcanic activity and many more. For this, the atmospheric  [28]. The figure presents the normalized counts in dependence of the projected scattering angle θy in degree. In PROPOSAL, 100 simulations each with 10 5 muons are performed for two different settings using the energy cut vcut = 10 −5 . The blue points present the mean of the simulations considering stochastic deflections and Molière scattering (MSM), the orange points present the mean of the simulations taking into account only Molière scattering. The uncertainties on the x-axis result due to the measured bin widths. The y-uncertainties are the standard deviations. The deflections are underestimated in PROPOSAL, except at θy ≈ 0°and at θy ≈ 3°. At deflections 2°< θy < 5°, the result seems to be more accurate than Geant4's. The consideration of the stochastic deflections shows no significant influence. muon flux is measured with a detector located below or even behind the object, which is visualized. The muon flux count rates then depend on the densities of the materials, the higher the density, the stronger the attenuation. Based on this, conclusions can be drawn about materials and cavities in the respective object. Sufficient statistics in reasonable time are obtained for muons in the GeV energy range [37].
Angular resolutions of these detectors are below 0.6°, which is on the order of magnitude of the muon deflection at GeV muon energies. Multiple scattering of muon deflection in several media is also studied in [37,38]. The resulting deflections are about 1°, similar to the deflections expected with PROPOSAL in water and ice. Hence, the scattering of muons can be a limiting factor for the angular resolution in muography at energies of a few GeV.

Conclusion
Stochastic deflection, recently implemented in PROPOSAL 7.3.0, is used to study the muon deflection per interaction. The deflection is dominated by multiple scattering except for a few stochastic outliers by bremsstrahlung. These angles are lower than ∼ 1°.
The results of PROPOSAL are compared with the common tools MUSIC and Geant4 and they are in good agreement. In low-Z materials, the region of outlier deflections fits the measured data better with PROPOSAL, than Geant4. A second data comparison points out that PRO-POSAL simulates more large deflections at higher muon energies and less at lower energies. In the data-MC ratio deviations up to a factor of 3 are observed. This points out that still improvements are required in the deflection parameterizations and in the multiple scattering, respectively. Since the presented measurements of the muon deflection are based on muon energies lower 10 GeV, deflection measurements of muons with energies up to TeV or even higher are required to validate the results at higher energies.
The median accumulated deflection depends primarily on the final muon energy, which can be interpreted as the muon energy at detector entry in neutrino telescopes or other muon detectors. The outcome is fit by a polynomial and can be used for a theoretical estimation of the muon deflection in water and ice. Since the result can be interpreted as the deflection before the detector entry, it defines a lower limit on the directional resolution. At energies lower 1 TeV, there is potentially a small impact of the muon deflection on the angular resolution of KM3NeT.  Each data set includes more than 50 000 events with the requirement that the true final muon energy E f is at most 10 % below the set final energy E f, min , E f > E f, min · 0.9. The energy cuts are ecut = 500 MeV and vcut = 0.05, and Molière scattering is chosen. Simulations are performed in ice, the deviations of the medians in a water-based simulation are smaller than 1 %. Since the medians overlap for different initial energies, there is no strong impact of the initial energy on the median deflection. These medians can be fit by a third degree polynomial in the log-space as shown in Eq. (2). The kinematic angle between the muon and neutrino is taken from Ref. [4]. Since the kinematic angle and the angular resolution of ARCA taken from Ref. [4] are presented in dependence of the neutrino energy as well as the resolutions of ORCA [32] and ANTARES [33], a rescaling to the muon energy is performed using the average energy transfer to the nucleus [34]. For energies E f ≤ 1 TeV, there is a minimal influence of the deflection on the angular resolution of ARCA and at E f = 5 GeV the upper limit of the deflections affects the resolution of SuperKamiokande [31]. The resolutions shown by IceCube [3], ORCA, Baikal [35] and ANTARES are not impacted, but no uncertainty bands are given for these either. The exact values and the propagated distances are presented in Table 3.