Atmospheric and Geodesic Controls of Muon Rates: A Numerical Study for Muography Applications

Abstract

Muon tomography or muography is an innovative imaging technique using atmospheric muons. The technique is based on the detection of muons that have crossed a target and the measurement of their attenuation or deviation induced by the medium. Muon flux models are key ingredients to convert tomographic and calibration data into the 2D or 3D density maps of the target. Ideally, they should take into account all possible types of local effects, from geomagnetism to atmospheric conditions. Two approaches are commonly used: semi-empirical models or Monte Carlo simulations. The latter offers the advantage to tackle down many environmental and experimental parameters and also allows the optimization of the nearly horizontal muons flux, which remains a long-standing problem for many muography applications. The goal of this paper is to identify through a detailed simulation what kind of environmental and experimental effects may affect the muography imaging sensitivity and its monitoring performance. The results have been obtained within the CORSIKA simulation framework, which offers the possibility to tune various parameters. The paper presents the simulation’s configuration and the results obtained for the muon fluxes computed in various conditions.

1. Introduction

The Earth’s atmosphere is bombarded by charged particles called primary cosmic rays. The least energetic of them come from the Sun, while for energies between 0.1 GeV and 10 GeV, the sources are of galactic origin (mainly supernovae). Beyond that, the particles are so energetic that the sources can only be of extragalactic origin and associated with more important events in the Universe (black holes, active galactic nuclei, gamma-ray bursts, pulsars, etc.). When entering the atmosphere, these primary particles will interact with the molecules of the atmosphere, generating atmospheric showers of different particles, which are called secondary cosmic rays. Among all charged particles reaching ground level, muons are the most numerous. Their energy loss or scattering when crossing a given length of matter is exploited to reconstruct the density of the medium. The rate of high-energy cosmic ray muons as measured underground is shown to be strongly correlated with upper-air temperatures during short-term atmospheric [1,2] phenomena and with pressure [3].

Applications of atmospheric cosmic rays (CR) have grown in numbers in the last decade in the field of muography. Measurements of the muons flux attenuation (absorption muography) or deviation (scattering muography) have been successfully applied to the imaging or monitoring of large geophysical [4,5,6,7,8], archaeological [9,10,11,12,13] or industrial structures [14,15,16,17,18]; however, the measurement of the muon energy is usually not possible in small and compact field detectors, which consist of standard trackers using particle physics technologies (scintillators, resistive plate chambers, micro-mesh gaseous structures, etc). The data post-processing depends on the knowledge of the incident muon flux.

Moreover, the simulation of the muon flux follows in general two different approaches: analytical models which are derived from experimental data and cosmic rays shower generators which use Monte Carlo (MC) techniques to simulate Extended Air Showers (EAS) from the primary cosmic ray to the particles at ground level. Among the analytical models, Tang et al. [19], Shukla et al. [20], Honda et al. [21], Guan et al. [22] and Gaisser et al. [23] are often used for first estimates of the muon flux at sea level or for few elevations only. These muon generators do not require computational resources but are not customizable. Since the aim of this paper is to present results obtained with a CR generator, we propose a more extended review on the most popular ones:

• MUPAGE [24] is a fast Monte Carlo generator of bundles of atmospheric muons for underwater/ice neutrino telescopes. It is based on parametric formulas obtained from Monte Carlo simulations of cosmic ray showers generating muons in bundle, with constraints from measurements of the muon flux with underground experiments. The range of validity extends from 1.5 km to 5.0 km of water vertical depth, and from 0° up to 85° for the zenith angle.
• Matrix Cascade equation MCeq [25,26,27,28,29] uses numerical equation cascades to study fluxes. It is a complete Monte Carlo calculation scheme, capable of calculating neutrino, electron and muon fluxes up to 100 TeV, with a statistical accuracy of about a few percent. All particles have their own cascades of equations that represent the evolution of the energy spectrum as a function of atmospheric depth.
• PARMA (PHITS-based Analytical Radiation Model in the Atmosphere) allows to estimate the terrestrial cosmic ray fluxes of neutrons, protons and ions, muons, electrons, positrons and photons almost anywhere on Earth and in the Earth’s atmosphere [30]. The model is based on analytical numerical functions whose parameter values are adjusted to reproduce the EAS results. The accuracy of the EAS simulation has been well verified with various experimental data.
• CRY (Cosmic-ray Shower Library) generates showers distributions for three observation levels (sea level, 2100 m and 11,300 m) and for primary particles from 1 GeV to ${10}^5$ GeV according to Hagmann et al. [31], and secondary particles from ${10}^{-3}$ to ${10}^5$ GeV. The showers are generated in a specific area (maximum size $300\times 300$ ${\mathrm{m}}^2$) from pre-computed tables as explained in Hagmann et al. [32] and primary protons are produced at an altitude of 31 km in the 1976 US atmosphere [33]. In this generator, the east–west effect is not taken into account but the latitude dependence with the geomagnetic cutoff and the CR spectrum modulation are provided. It is possible to set the type of secondary particles to be studied, the altitude, the latitude, the date, the number of particles and the size of the surface of interest. The date allows to take into account the solar modulation described by Papini et al. [34]. It is possible to use pre-calculated tables from GEANT4 to take into account the configuration of the detector [35]. CRY has limitations when you investigate multi-track events in a cosmic ray experiment or identifying background events in muography.
• CORSIKA (COsmic Ray SImulations for KAscade) [36,37] is a Monte Carlo code for simulating atmospheric particle showers initiated by high-energy cosmic ray particles. Primary particles (protons or light nuclei) are tracked in the atmosphere until they interact, decay or are absorbed. All secondary particles are explicitly followed along their trajectories and their parameters are stored when they reach an observational level.

All the models presented above have their own particularity and the choice of the fluxes generator is essential. The sensitivity of the absorption muography technique relies on the model’s accuracy, which should take into account the experimental conditions in the most realistic possible way. For instance, it has been shown in Jourde et al. [3] and Tramontini et al. [2], that atmospheric conditions strongly impact the muon flux. It is also clear that the geomagnetic field may play a significant role by deflecting charged particles towards the poles [38], which leads to a decrease in the flux at the equator and an increase at high latitudes. Furthermore the altitude of the experiment plays a role since it requires to measure the particle flux at a different stage of the cosmic shower development.

In this article, we detail the methodology and the results obtained with the CORSIKA simulation framework. In Section 2, we define the simulation main configuration parameters concerning the primary flux, the detection device, the hadronic models, the atmosphere and the geomagnetic field properties. In Section 3, we contrast the results of the simulation with analytical models, experimental data and numerical results on the effects of elevation effect, geomagnetic field and atmospheric conditions (pressure, temperature, etc.) on the muon flux are presented. Finally, in Section 4, we draw a global conclusion on the obtained tool performance.

2. Materials and Methods

2.1. Standard Use Cases

Several developments use CORSIKA for the muon fluxes simulation and for various applications. In Biallass et al. [39], they based their parametrization of the differential muon flux at ground level on a flux obtained from the air shower simulation program CORSIKA. According to Wentz et al. [40], the number of simulated particles should be normalized by taking into account the primary flux. In Mitrica et al. [41], they used this generator to determine accurately the overburden in mwe (meter water equivalent) of an underground laboratory. Pethuraj et al. [42] used CORSIKA to generate secondary particles at an experimental site where atmospheric neutrinos oscillations are studied in order to assess the long-term stability of their RPCs detectors. They studied the azimuthal dependence of the muon flux at different zenith angles and they compared they results using CORSIKA with Honda et al. [21]. For a different purpose, Kovylyaeva et al. [43] studied variations of the extra-atmospheric origin CR by introducing atmospheric effects corrections. They used CORSIKA to estimate barometric and temperature coefficients for various components of CR. Finally, Useche et al. [44] chose CORSIKA to gauge the expected cosmic ray muons flux attenuation by a structure. We use those users’ feedback to tune our CORSIKA’s options that are believed to be essential to simulate the muon flux in a specific way and we detail some of them below.

2.2. Simulation’s Configuration

We follow the muon production path from the incident primary particles to the detection level. Detailed knowledge on their properties is an important issue for the secondary particles generation.

2.2.1. Primaries at the Top of the Atmosphere

The primary cosmic ray flux is composed of several types of particles (H, He, C, O, Fe, etc). The associated fluxes are different by several orders of magnitude and the most influential flux is that of the single proton. Nevertheless, we considered different components to obtain an accurate primary flux. Several analytical models are used to describe the behavior of primary cosmic rays where the total flux $J\left(E\right)$ is the sum of its various constituents and it is expressed in term of energy per nucleon (Papini et al. [34]) or per nucleus (Hörandel et al. [45], Wiebel-Sooth et al. [46]). They are our possible choices when we normalize our simulated fluxes of secondary particles with primaries spectrum. It is possible to run simulations for each primary particle type or to apply what we call the superposition principle explained by Spurio [47]. We choose the second way and we only simulated the arrival of primary protons in the atmosphere to save computation time.

Primary CR models are compared on Figure 1, their amplitudes are different depending on the energy of the primary particle. Wiebel et al. [46] seems to overestimate the flux across the full range of energy compared to Papini et al. [34]. Hörandel et al. [45] has a different behavior for energy ranges; it overestimates as well as underestimates. Figures in this paper are made using the Papini’s model, which produces an “intermediate” primary flux amplitude when compared to the others and allows to take into account the solar modulation on the expected primary CR.

• Number of simulated showers and energy ranges

It is possible to choose the number of simulated showers per run. It is preferable (knowing the shape of the primary CR particle spectrum) to have a larger statistic for low energy particles. We choose the energy range of the primary particle to match the muon energy measurable in our tomography experiments: between 1 and 10${}^7$ GeV. This total range is split into several parts with a defined step. The number of simulated primary particles is weighted in each bin with an empiric law for the number of showers: $N={10}^{9-E}$ with $E=[1,7]$. We made this choice in order to obtain the most accurate flux possible by launching enough showers. The duration and the memory requirement of the simulations is also a factor to take into account, so we must limit ourselves to a reasonable number of showers. To be consistent, our simulated muon fluxes will be corrected by known analytical primary fluxes to cover the whole chosen spectrum.

• Slope of the primary spectra

The energy spectrum of primary CR follows a roughly exponential law: $E^{-\gamma }$. Usually $\gamma$ is fixed to 2.7 for an “intermediate” energy range of primary CR spectrum according to Nesterenok et al. [48] and Pethuraj et al. [49]. In our case, $\gamma$ is calculated for each energy range and run with a fit defined by an analytical model describing the primary spectrum. We selected a cosmic ray primary model based on real data fits between Papini at al. [34], Hörandel et al. [45] and Wiebel et al. [46]. Indeed, the spectrum of primaries is different from one model to another and the slope of the curve depends on the primary CR model (see Figure 2).

• Particles angles trajectories

On top of energy and momentum parameters, the trajectories of particles are defined by their zenith and azimuth angles, ranging from 0° to 90° and from −180° to 180°, respectively. The coordinate system used in CORSIKA is shown in Figure 3. Furthermore, for zenith angles higher than 60° the curvature of the atmosphere is taken into account, it cannot be neglected [50].

2.2.2. Hadronic Interaction Models

When a primary particle reaches the top of the atmosphere, it undergoes “hadronic” interactions leading to the production of secondary CR and among those particles, pions and kaons decay into muons. Different types of primary particle interaction models are available on CORSIKA. We chose FLUKA for the low-energy interactions and QGSJET-II-4 for high energies, the best candidates in their energy domains for our expected accuracy and computational time.

To be more specific, FLUKA is a hadronic interaction model of low energy particles used by Usokin et al. [51], Apel et al. [52] and Heck et al. [53]. Its main strength is to allow for rapid simulation. QGSJET is a high energy hadronic interaction model. It has a high level of sophistication, approximates the data very well but is a bit slow especially because it allows few free parameters and deals with the effect of saturation by partons [54]. It is notably chosen by Tapia et al. [55] and Fedynitch et al. [26] for their simulations. In Atri et al. [56], they combine it with FLUKA.

2.2.3. Particles Energy Thresholds

To decrease the computation time and to stay compatible with the detection efficiency of our detector, we set energy thresholds values under which particles (hadrons, muons, electrons and photons) are no longer taken into account in order to simulate only the muons (or their parents) detectable by our telescopes. The choices are different in the bibliography according to the

Table 1. Energy thresholds values used by different authors and for this study.

E Threshold	Hadrons	Muons	Electrons	Photons
1. [57]	0.1 GeV	0.1 GeV	0.1 MeV	0.1 MeV
2. [58]	0.2 GeV	0.2 GeV	0.1 MeV	0.1 MeV
3. [59]	0.05 GeV	0.05 GeV	0.003 GeV	0.003 GeV
4. [52]	300 MeV	100 MeV	250 keV	250 keV
5. [60]	0.05 GeV	50 keV	50 keV	50 keV
6. This study	0.05 GeV	0.01 GeV	0.001 GeV	0.001 GeV

.

2.2.4. Detector Geometry and Observation Levels

CORSIKA offers three types of detector (flat, volume and vertical) associated with different angular distribution of the atmospheric shower and flux intensity [37]. For our muography detectors, the “volume” and “flat” types are preferred.

The altitude of the detector is defined by the observation levels (in cm); up to 10 levels for a simulation.

2.3. External Parameters

2.3.1. Earth’s Magnetic Field

The Earth magnetic field is created by its magnetosphere, which reduces the intensity of the high-energy flux reaching the ground. The geomagnetic field (B) modifies the spectrum of particles bombarding our atmosphere with a low-energy cutoff. The Earth’s magnetic field is able to deflect primary CR below 10 GeV near the equator and close to 1 GeV at higher latitudes. The primary CR intensity also varies with longitude because of the asymmetry of the geomagnetic axis with respect to the Earth’s rotation axis [40,49]. Those “east–west” fluxes show differences in energy intensity up to 100 GeV. This difference is more marked at high altitude than on the ground. Finally, there are significant local variations of the geomagnetic field, which affect the intensity of CR, the most famous being the South Atlantic Anomaly (SAA).

For each place, we declare the horizontal $B_x$ and vertical $B_z$ components of the Earth’s magnetic field Z (in µT), $B_x$ being the magnetic north H (see Figure 4). They are generated by the NOAA geomagnetic calculator [55]. CORSIKA computes the total magnetic field and its inclination from these two components.

2.3.2. Atmosphere Parameters

An input to CORSIKA is the state of the atmosphere in which the CR ray showers are generated. The CORSIKA atmosphere model is composed of N₂, O₂ and Ar with volume fractions of 78.1%, 21.0% and 0.9%, respectively. The density variation of the atmosphere with altitude is modeled by 5 layers. The state of the atmosphere is described by the density of the air at each altitude level. This parameter is calculated by converting the relative humidity into saturation vapor pressure with the Magnus formula [61]. We computed the parameters and altitudes of the layer boundaries from ERA5 data, the latest climate reanalysis produced by the ECMWF (European Centre for Medium-Range Weather Forecasts), which combines large amounts of meteorological observations with estimates made from advanced modeling and data assimilation systems [62]. The choice of ERA5 is guided by the abundance of pressure levels and the availability of data at higher altitudes.

2.3.3. Physical Validity Range

The muon rate at sea level is a logarithmic function spamming multiples order of magnitude in the energy range, and a $cos{\left(\theta \right)}^2$ function in the angular range. The simulations encounter difficulty at the end of each of these spectrum’s range. We want to stress the fact that the estimates have well known limitations in the 70–90° angular range. They are quite limited above ${10}^5$ GeV, where the rarity of muons is not negligible. Since the end of range in energy is not a particular issue for the muography perspectives, we consider this to be satisfactory and therefore, the chosen energy range is between 1 and ${10}^6$ GeV.

2.4. Normalization Issues

From the CORSIKA output files we calculate the muon rate, according to the zenith angle and the muon energy, based on ${10}^{9-E}$ showers (see Section 2.2.1) that correspond to the validity range. Muon rates are then corrected by different factors as primary flux, geometric factor and normalization linked to the number of showers generated by the energy range.

The differential muon energy flux is written as:

$$\Phi (E_{\mu },\theta \mu )=F\left({\theta }_{\mu }\right)\times \sum _{i=1}^{N_{prim}}\sum _{j=1}^{N_{sim}}{w}_{i,j}\times \frac{H_j(E_{\mu },{\theta }_{\mu })}{N_{sh,j}\times \Delta E_{\mu }}$$

(1)

where $N_{prim}$ is the number of primary particles we want to consider in the flux calculation. $N_{sim}$ is the number of simulation j. H is the histograms created as a function of the energy E and the zenith angle $\theta$ during each simulation j. $\Delta E_{\mu }$ is the width of the energy bin of the histogram. The width of the zenith angle is taking into account in $F\left({\theta }_{\mu }\right)$. Variations due to the azimuth angle are not taken into account here. In the interval j, we simulate $N_{sh,j}$ defined in Section 2.2.1.

$w_{i,j}$ is the normalization factor. It takes into account the contribution of the primary particle flux i [40] in the energy interval defined by the run j. We perform all our simulations with hydrogen and apply the principle of superposition: a nucleus with atomic mass number A and energy E is equivalent to a single nucleon with energy $E/A$. In other words, a shower initiated by a particle i containing $n_i$ nucleon, of energy $n_i\times$ E (GeV) is equivalent to $n_i$ showers of protons (of hydrogen) of energy E. The method is validated by summing the contribution of each primary particle on the muon flux and comparing it to the muon flux obtained with the superposition model. A difference is observed at low energy but since CORSIKA requires an energy higher than 80 GeV/nucleon for the primary particles other than hydrogen, the results are correct.

In order to consider the solar modulation of 4.5% (up to 100 GeV) it is recommended to use the parametrization of Papini [34] to modulate its primary particle spectrum; however, we have seen that Hörandel [45] and Wiebel [46] were two possible alternatives.

In the Papini model [34], the primary flux J is expressed in terms of energy per nucleon (in GeV/A):

$$w_{i,j}=A_i{\int }_{E_j^{inf}}^{E_j^{sup}}{J}_i\left(E\right)dE$$

(2)

$A_i$ is the atomic mass of the primary particle i.

F(${\theta }_{\mu }$) is the geometric correction factor related to the solid angle $d\Omega$ adapted to the type of detector (see [47]). For example, in the case of a flat detector: $d\Omega =sin\left(\theta \right)cos\left(\theta \right)d\theta d\psi$ so:

$$F\left({\theta }_{\mu }\right)=\frac{1}{2{\int }_{{\theta }_{inf}}^{{\theta }_{sup}}sin\left(\theta \right)cos\left(\theta \right)d\theta }$$

(3)

To summarize, the normalization takes into account the sampling effects of the histograms, the number of showers $N_{sh,j}$ generated in each simulation j, the geometric effects related to the zenith angle and the modulation of the primary i contribution over the energy range [$E_{inf,j}$, $E_{sup,j}$].

The integrated energy flux, for a minimum energy $E_{min}$ is obtained from the differential flux thus normalized:

$$\Psi (E_{min},\theta )={\int }_{E_{min}}^{\infty }\Phi (E,\theta )\times dE$$

(4)

The discrepancy and uncertainties of muon rates are estimated using these simulations, using mean and standard deviation estimates.

3. Results

3.1. Comparison with Analytical Models

Our fluxes have been cross-checked with the analytical models listed in the Introduction. We present the results of the comparison with Tang et al. [19] regarding their integrated flux ratio $I_{\mu }/I_{\mu ,2}$ as a function of the muon zenith angle and for several muon energy ranges (different colors) in Figure 5. This figure presents the normalized intensity distribution when considering different energy ranges: 10–100 GeV, 100–${10}^3$ GeV, ${10}^3$–${10}^4$ GeV and the distribution for the whole energy range. The histogram of varying flux ratio values over zenith angle is meant to highlight the existence of systematic bias over the entire range of zenith values. It shows that when considering a larger energy range, the intensity ratios distribution tends to a narrow peak centered at 1.0 ± 0.2; on the whole range, the fluxes have the same intensity. Right panel shows the flux ratio for different zenith angles $\theta$ in the range 0 to 90°, when considering the same energy ranges. The models fit well on intermediate energy ranges (10–${10}^2$ GeV in dark blue) except at high zenith angles; however, their intensity is different at high energy (${10}^2$–${10}^4$ GeV, in light blue and green). The 1–10 GeV is range not shown: the ratio starts around 4 at 0° and increases with the zenith angle following a power law because the Tang flux is unstable at low energy and especially at high zenith angles. Our CORSIKA relative fluxes are plotted in Figure 6 only up to 2000 GeV because beyond this point our fluxes are too unstable. The part of the most energetic muons is more important than the low energetic ones at high zenith angles $\theta$ and the opposite occurs at low angles. Low-energy muons are important for calibration and high-energy muons for conducting tomography experiments; however, analytical models are known to be poorly adapted to small and large energies, because few measurements are available for their fitting equations. The CORSIKA model instead is probably more reliable over the whole energy range. Furthermore, analytical models are not extrapolated for all zenith angles and they do not take into account geodesics parameters, a limitation overcome by the CORSIKA approach.

3.2. Comparison with Real Data

The best cross-checks for any simulation are the comparison to real data. Data presented here were taken in Lyon (France, latitude 45° and close to the sea level), in open sky conditions, with a 3-planes muon tracker (so-called muon telescope). We tilted the telescope progressively by step of 15° from the zenith ($\theta =0$°) to the horizontal ($\theta =90$°) directions. Multi-tracks are suppressed, electrons are partly suppressed by the multiplicity cuts on the detection planes, but there is no PID in the system. The muon flux was simulated in Lyon (France) respecting the geodesic constraints. Figure 7 displays the data/simulation comparison. Experimental points are represented with crosses in blue and the flux from CORSIKA’s simulations is represented in black. A fit was made on the real dataset with a simple $co{s}^2\left(\theta \right)$ at first order (in red). Despite the still pending disagreement at large zenith angles, we observed a real improvement in the data/model comparison with respect to the usual $co{s}^2\left(\theta \right)$ (or analytical) fit. Horizontal muons are best described by CORSIKA.

3.3. Geodesic and Atmospheric Factors

In this section, we used our CORSIKA model, validated on a large dataset to simulate muon fluxes for various particle detectionaltitude, Earth magnetic field B and atmosphere conditions. We wanted to quantify those effects on the measurable muon flux at open sky conditions. To achieve this, we varied one parameter at the time while keeping the others constant. All the results presented in this section are subject to significant uncertainty. It is statistical and increases with the energy. Indeed, we simulate much less high-energy muons and extreme-angles muons. We have limits imposed by computational time and linked to the randomness of CORSIKA. Our errors are all set to 1$\sigma$.

For our preliminary tests we chose the city of Lyon, France (lat = 45.75° N, long = 4.75° E, $B_x$ = 22.71 µT, $B_z$ = 40.96 µT.) where we performed the open sky measurements shown in Figure 7. It is located at a medium latitude with a temperate climate. Then, we have used the parameters of Dallol, Ethiopia (lat = 14.14° N, long = 40.18° E, $B_x$ = 36.2752 µT, $B_z$ = 10.1924 µT where the hottest temperatures on Earth have been recorded. Oimiakon, Russia (lat = 69.18° N, long = 141.6° E, $B_x$ = 9.9943 µT, $B_z$ = 58.5077 µT is a plateau where the coldest temperatures in the northern hemisphere have been recorded. These two places were chosen to study extreme weather effects.

3.3.1. Altitude Effects

The altitude of the experiment plays an important role as it involves measuring the particle flux at a different stage of the cosmic shower development. We compared the muon flux in Lyon, France at 1000 m with the muon flux at sea level. Figure 8 presents the normalized intensity distribution when considering different energy ranges: 1–10 GeV, 10–100 GeV, 100–${10}^3$ GeV, ${10}^3$–${10}^4$ GeV and the distribution for the whole energy range. The histogram of varying flux ratio values over zenith angle is meant to highlight the existence of systematic bias over the entire range of zenith values. It shows that when considering a larger energy range, the intensity ratios distribution tends to a narrow peak centered at 1.17 ± 0.01, which means that there are 17% more muons at 1000 m than at sea level in overall. Right panel shows the flux ratio for different zenith angles $\theta$ in the range 0 to 90°, when considering the same energy ranges. It shows that at an elevation of 1000 m the flux ratio is much higher for low energies (1 to 10 and 10 to 100 GeV, in dark and light blue), and less important at high zenith angle (80 to 90°). For higher energies, the flux ratio remains more constant at each angle (in green and yellow).

At high altitude, low energy muons fluxes are more important than at sea level. Low energy muons present at high altitude will probably decay before reaching the ground. The muons of high zenith angle travel longer distances in the atmosphere and lose more energy. This is why we find more low energy and high zenith angle muons on the ground than at 1000 m and why the part of high zenith angle and energy muons decreases. Finally, depending on where you are on the globe, the ratio of fluxes at a given altitude to sea level varies. By looking at two different places on Earth, Dallol (Ethiopia) and Lyon (France), we have observed a different ratio between “flux at 1000 m/flux at sea level”. This means that other factors than altitude influence the muon flux. Magnetic fields and atmospheric conditions are different between these two places and probably affect the integrated muon flux on the column in a distinct way. It seems that the difference increases as the latitude increases.

3.3.2. Geomagnetic Field Effects

In this part, the effects of the magnetic field are tested. For this purpose, the atmospheric and altitude parameters have been kept constant between the compared fluxes. We have set the vertical component as a constant ($B_z$ = 20 µT) and the horizontal component ($B_x$ = 15 µT) for a simulation of reference. We then compute another simulation by only changing $B_x$ = 45 µT. We simulated the muon flux with those two configurations and computed their ratio for four different energy ranges and different zenith angles.

Figure 9, on left panel, presents the normalized intensity distribution of the flux ratio over the zenith angle, for each individual energy range. The histogram of varying flux ratio values over zenith angle is meant to highlight the existence of systematic bias over the entire range of zenith values. It shows that when considering a higher energy range, the intensity ratios distribution (1–${10}^4$ GeV) tends to a narrow peak centered at 0.95 ± 0.01; the median value, corresponding to a small effect of about 5%. As expected the ratio tends to one for energy ranges greater than 10 GeV but affects the low-energy particles that are more deflected towards the poles. The right panel shows the flux ratio with respect to the zenith angle $\theta$ for fixed energy ranges. At low energies (1 to 10 GeV in dark blue) the flux is higher for $B_x$ = 45 µT until 60° and between 85° and 90°, and lower between 60 and 85°. This probably arises from the fact that high-angles particles have to cross a larger section of the atmosphere and therefore may start their travel with higher energies, making them less sensitive to the geomagnetic field’s effect. Note that this sample dominates over the total integrated flux ratio (in black), which follows more or less the same behavior. As expected, at high energies, typically above 10 GeV, the flux ratio remains more constant with the zenith angle (in light blue, green and orange). Note that this high energy sample is usually not used in scattering mode; therefore, one has to pick corrections for the geomagnetic effects to perform absolute measurement with scattering muography.

The same procedure is followed this time with a fixed horizontal component and various vertical components. We simulate the muon flux with these two configurations and we compute the ratio as a function of the energy and of the zenith angle. Only a slight difference at high zenith angles is observed. The $B_z$ component does not seem to affect the muon flux.

3.3.3. Atmospheric Thermodynamics Effects

We have tested the effects of altitude and magnetic field and their components on the muon flux. Finally we studied the impact of atmospheric density on the flux, which is ultimately controlled by the temperature and pressure state of the atmosphere. Hence, at the same location, seasonal variations affect the muon flux over time. For these tests we chose Dallol (Ethiopia) and Oimiakon (Russia) to test the effects of extreme temperature. These atmospheres are presented in the next subsection. The atmosphere density parameters used to simulate muon fluxes with CORSIKA are determined by using ERA5 datasets.

• Atmosphere simple parametrization

Temperature and density profiles of Dallol, Oimiakon and Lyon are displayed in Figure 10. In Figure 11, relative humidity profiles are plotted for December 30th. We observe that the density of the atmosphere decreases with altitude, and that colder atmospheres are denser especially at lower altitudes. For the relative humidity profiles shown in Figure 11, the water vapor content in the air varies with altitude. Similar to temperature, relative humidity also depends on the pressure of the given system. It is low at ground level but important at high altitude. In very hot and/or dry atmospheres, it never reaches 100%.

• Atmosphere properties

The atmosphere can be divided into five layers: the troposphere, stratosphere, mesosphere, thermosphere and ionosphere. The atmospheric profiles presented in Figure 10 stop at the end of the stratosphere (∼50 km). The troposphere is the part of the Earth’s atmosphere located between the surface and an altitude of about 8 to 15 km, depending on the latitude and the season. It is thicker at the equator than at the poles. This layer concentrates three quarters of the atmospheric mass and the temperature decreases rapidly with altitude. The stratosphere extends, on average, between 12 to 50 km. It is characterized by an increase in temperature with altitude. The stratosphere begins at low altitude near the poles, because the temperature is lower there. The distribution of atmospheric density is therefore different at opposite latitudes.

Despite the fact that the temperatures are very different at ground level, all temperatures decrease with increasing altitude (up to 10–15 km altitude) regardless of their location. Then, the temperatures increase in the stratosphere and decrease in the mesosphere. We notice that the change of behavior of the curve at Oimiakon is made at lower altitude. Indeed, the stratosphere starts lower at the pole. The relative humidity (RH) increases globally with altitude and is maximum at Oimiakon at high altitude, at 2 km Dallol have the highest RH and on the ground the values are almost the same everywhere on the globe in Figure 11.

Muons are typically produced at an altitude of 10–15 km (troposphere/stratosphere boundary). Their abundance is affected by the density differences in the atmosphere either by direct re-interaction or by modification of their parent mesons survival probabilities before decay [63,64]. The effect is more important for high-energy muons, which result from high-energy mesons with larger lifetime due to time dilation and therefore with longer paths in the atmosphere; thus, high-energy muons rates are more sensitive to temperature changes.

Figure 12 presents the normalized intensity distribution when considering different energy ranges: 1–10 GeV, 10–100 GeV, 100–10${}^3$ GeV, 10${}^3$–10${}^4$ GeV and the distribution for the whole energy range. The histogram of varying flux ratio values over zenith angle is meant to highlight the existence of systematic bias over the entire range of zenith values. It shows that when considering a larger energy range, the intensity distribution tends to a narrow peak centered at 0.90 ± 0.26, the median value, which means that the flux is 10% higher in Oimiakon (Russia) than in Dallol (Ethiopia); however, this statement is only valid for the integration over the total energy range. Right panel shows the flux ratio for different zenith angles $\theta$ in the range 0 to 90°, when considering the same energy ranges. It shows that in Dallol the flux is higher for high energies (100 to 10${}^4$ GeV, in green and yellow), and lower in Oimiakon. For lower energies, the flux ratio is higher in Dallol (1 to 100 GeV, in light and dark blue), and lower in Oimiakon. These effects increase with the zenith angle $\theta$. Statistical errors of fluxes simulations are fixed to 1 $\sigma$.

Knowing that Dallol is hotter than Oimiakon, we understand that the flux is higher in Oimiakon for high energy muons. At low energy, the integrated flux is higher in Dallol.

We estimated and summarized the maximum observed variations for the studied effects: magnetic field, altitude and state of the atmosphere (extreme and seasonal effect: summer vs. winter) for several muon energy ranges in

Table 2. Maximum observed variations comparing two CORSIKA fluxes for the studied effects: magnetic field, altitude and state of the atmosphere (extreme and seasonal effect) for several muon energy ranges.

	1–10,000 GeV	1–10 GeV	1000–10,000 GeV
Magnetic field	5 (±4)%	6 (±1)%	0.2 (±1)%
Altitude: - 1000 m/0 m - 5000 m/0 m	15 (±1)% 115 (±10)%	17 (±2.5)% 106 (±25)%	7 (±2.7)% 2 (±1)%
Atmosphere: - Extreme - Seasonal	10 (±7)% 8 (±1)%	13 (±10)% 9 (±3)%	5 (±10)% 2 (±10)%

. Altitude effect is the most impacting and the observed flux fluctuations are all dependent on the muon energy ranges for all effects. Seasonal effects have less impact than extreme atmospheres.

4. Conclusions

In this paper, we presented a muon flux simulation workflow accounting for muon-atmosphere interactions, based on the CORSIKA framework. We detailed our simulation strategy and the various relevant inputs from the hadronic interactions models to the atmosphere conditions. In particular, we used meteorological ERA5 pressure and temperature datasets to compute the required atmospheric density profiles. The workflow has been cross-validated against experimental evidences and standard semi-empiric models found in the literature. Simulations prove themselves to be a powerful tools to study and make predictions on tiny effects induced by the geomagnetic field, detection altitude or the atmospheric seasonal variations. The altitude effect seems to be the most important. In Lyon for example, at 1000 m altitude, the flux is 17% higher than at sea level. At medium latitude, seasonal effects affect the muon flux by modifying these proportions at high and low energy for a higher overall integrated flux at low temperatures (in winter). By analyzing the role of the temperature (and the atmospheric density) in this process with extreme atmospheres, we find the same results (∼10%). In contrast, the short term effects, such as an atmospheric depression, are not visible. An increase in the geomagnetic field exclusively increases the flux of low energy muons (∼5%) and only the horizontal component $B_x$ variation seems to play a role. Atmospheric impact appears to be dominant over the geomagnetic field effect for opposites latitudes (poles/equator). Those effects are of increasing importance when one wants to produce muon imagery and/or long-term internal state surveys, on both ends of the opacity spectrum. Low-opacity targets imagery is controlled by low-energy muons filtered out by the density of the atmosphere. On the other side, high opacity targets imagery is largely affected by the process at stake for high-energy muon production. This study opens the gate to develop semi-empiric formulas predicting the evolution of the muon energy spectrum for each zenith angle, in relation with the atmospheric state. These formulas will be useful to correct recorded muon fluxes on the fly, when direct open-sky measurements are not available or not sufficiently refined in terms of energy description.