On the mystery of the multi-muon ﬂux at the TeV cosmic-ray energy range

. Current Monte Carlo simulations do not provide a good description of the muon component of extensive air showers. Many air shower experiments report discrepancies between their data and Monte Carlo predictions, ranging from the TeV scale up to the highest energies. In these proceedings, we address the seasonal variation of the multi-muon events observed by the NOvA Near Detector (ND). For our studies, we use the general-purpose Monte Carlo code FLUKA to treat the transport and interaction of the air-shower particles in the atmosphere and other media. Our design considers a multilayered atmosphere and a layered underground approximated to match the NOvA ND location and detector geometry. Our atmospheric model uses air densities for winter and summer calculated from the temperature and geopotential information for the pressure levels given by the European Center for Medium-Range Weather Forecasts (ECMWF) datasets in situ. Understanding the multi-muon ﬂux at the cosmic ray high-energy range may lead to a better description of the muon production mechanisms in ultra-high-energy extensive air showers. In addition, it can help to improve future Monte Carlo codes or hint at new physics processes or interactions.


Introduction
The seasonal variation of the underground muon flux, originating from the interaction of high-energy cosmic rays in the Earth's atmosphere, has been studied by the MINOS and NOvA near detectors (NDs) located 100 m underground [1,2]. MINOS cosmic data [3] show that the total (single)-muon flux increases in summer and decreases in winter by ≈ 1%. In contrast, both experiments find the multi-muon event rates to reach a maximum during winter and a minimum in the summer [1,2]. The magnitude of the seasonal multi-muon flux oscillation is five times larger than the single-muon flux oscillation and has an opposite phase: the maximum is shifted by approximately six months. In [2], the MINOS Collaboration proposed several mechanisms to explain the winter maximum for the multi-muon flux. However, none of them fully describe the magnitude of the oscillation. A recent work [5] using CORSIKA simulations [6] has also addressed this issue. Still, a quantitative agreement was not obtained: the predicted seasonal multi-muon oscillation magnitude (see Fig. 9, in [5]) is a factor four times smaller than the MINOS and NOvA data, showing that even in the wellknown TeV range, discrepancies between data and Monte Carlo simulations exist. A proper understanding of this problem could help to improve future Monte Carlo codes or even hint at new physics processes in hadronic interactions. It could also contribute to the work of the WHISP [7] group reporting on the muon puzzle arising at √ s ∼ 10 TeV.
In this work, we use the general-purpose Monte Carlo code FLUKA [9] to treat the transport and interaction of the shower particles in the atmosphere and other media. In particular, we have based on the FLUKA's 100-layered atmospheric model, where we have implemented the air density for the summer and winter profiles as provided by the ECMWF database [10] for the NOvA/MINOS NDs location. The muon component was analyzed at the surface (226 m above sea level) and the detector plane (99 m underground) for both atmospheric profiles. The atmospheric model and geometry are described in section 2. Preliminary results are discussed in section 3. A summary of our results is given in section 4.

Atmospheric Model and Detector Geometry
Our atmospheric model is based on the FLUKA's 100layered model [8], adapted for the NOvA ND geometry and location, considering different atmospheric density profiles for winter and summer. Since our detector area is ≈ 20 km × 20 km at the surface and ≈ 1 km × 1 km at the detector plane, we can neglect the effects of the Earth's curvature and have adopted a cylindrical geometry instead of the one from the original model. Our layout includes four further underground layers, allowing for more proper treatment of the interaction and propagation of the muons underground.
The description of our atmosphere starts at 72 km altitude, with a density of 8.78 × 10 −8 g cm −3 for both atmospheric profiles. The air density is gradually increasing over the one hundred layers, reaching 1.21 × 10 −3 and 1.28 × 10 −3 g cm −3 at the surface (226 m a.s.l.) for the sum- mer and winter profiles, respectively. The air densities for winter and summer have been calculated from the temperatures and geopotential information for 37 pressure levels given by the ERA5 ECMWF dataset in situ (averaged data from January and July 2017 periods were used). We have implemented a simple uniform magnetic field permeating the whole volume with components B The first underground layer has a thickness of 99 m, ranging from the surface level down to an altitude of 127 m a.s.l., where the NOvA ND experiment is installed. It has a molasse composition with a density of 2.35 g cm −3 . At this depth, we implement the second underground layer, which consists of a cavity of 2 m height filled with air, with a density of 1.205 × 10 −3 g cm −3 . The end of the second layer coincides with the altitude of the top of the NOvA detector, at 125 m a.s.l., followed by a four meters layer of mineral oil (1 g cm −3 ). Finally, we add a molasse layer that reaches the sea level, which further slows down and absorbs the particles, as illustrated in Figure 1.

FLUKA cosmic ray shower library
For the work presented in this paper, we have used the Monte Carlo code FLUKA [9], version 4-2.2, with the DP-MJET package [4] as the hadronic interaction model. We simulate a flux of high-energy cosmic rays arriving from two zenith angles, θ = 30 • and 50 • , at an altitude of z ≈ 72 km, which coincides with the top of our atmosphere. We adjust the XY initial coordinates of cosmic rays so that the shower axis intersects the center of the surface plane, located at 125 m a.s.l., at x = y = 0. For each zenith angle, we simulate six energies: 10, 30, 50, 100, 500, and 1000 TeV. In order to mitigate the statistical uncertainties at the lowest energies, where few muons reach the NOvA detector, we have produced a shower library in which the number of simulations decreases with energy as follows: 10 4 showers at 10 TeV, 4 × 10 3 at 30 and 50 TeV, and 10 3 at 100, 500, and 1000 TeV. We are assuming a pure proton cosmic ray mass composition since protons are 5 -6 times more abundant than helium nuclei [11] at this energy range (with reference to energy per nucleon)

Preliminary Results
In this section, we report on our findings regarding the differences in the muon fluxes for the winter and summer atmospheres. In subsection 3.1, we present our results for the total flux of muons (single-muon flux), and in subsection 3.2, we focus on the multi-muons case, which is a sub-sample of the total muon flux. We use the following notation for all our figures and tables: winter atmosphere results are denoted by a W and plotted in blue, while for the summer, we will use an S and red color. For our analysis, we used the FLUKA shower library described in subsection 2.1, from which we stored the information regarding the position and momenta of all muons produced along the shower development.

Single-muon flux
We observe that the total muon flux at the surface is always higher in the winter than in the summer. On average, 12 − 22% (10 − 18%) more muons are detected at the surface in winter than in summer for θ = 30 • (50 • ) showers. This difference was found to increase with the cosmic ray energy. We find that our results are in agreement with the observed seasonal muon flux variation reported in [12].
In a second step, we calculated the muon flux summerto-winter ratio at the detector level, located at 99 m depth. We found that, at the detector plane, the ratio becomes compatible with one within statistical errors, meaning that a seasonal variation is not significant. A summary of the results discussed above for the surface and detector planes is given in Tables 1 and 2.  We have calculated the fraction of muons reaching the detector plane from the surface to be 3.2 − 6.1% (2.8 − 5.4%) in winter (summer). Assuming that this ratio is mainly related to the decay of low-energy muons propagating in the rock, the energy spectrum of muons should be softer in summer than in winter. One possible explanation is that, in winter, the atmosphere is more compressed due to lower temperatures leading to cosmic rays interacting deeper and producing air showers developing in a denser medium. That leads to an increased probability of pions and kaons interacting more times in the atmosphere before decaying into (lower energy) muons. In Figure 2  We have also calculated the mean distance of muons from the shower axis as σ r = (x 2 i + y 2 i ) 1/2 at the detector plane for both profiles. Our results are shown in Table 3. From Figure 2 and Table 3, we confirm that the lateral distribution of muons is slightly denser in winter, suggesting that it could play a key role in the seasonal variation detected by the NOvA ND experiment [1].

Multi-muons
To reduce the statistical uncertainties in the measurement of multi-muon events in each shower, we used a grid of identical NOvA-ND detectors (16 × 4 m large) located in |x| < 400 m; |y| < 400 m area, thus effectively simulating a total of 10 4 detectors. For each shower, we check for multi-muon events in each one of the detectors. This way, we simulate the randomly shifted position of NOvA-ND relative to the distance from the shower axis.
We have evaluated the average multi-muon flux for the whole simulated shower library. The number of multimuon events is shown in Tables 4 and 5, where s = n S and w = n W are the number of multi-muons detected in summer and winter, s/w is their ratio and D SW quantity is the S − W difference in percentage of S + W (which is two times the average flux), as given by: with statistical error δ evaluated as: A negative value of D SW means a deficit of multi-muons in summer with respect to winter.  to 50 TeV and decrease afterward. For θ = 50 • , the maximum difference is obtained at 30 TeV, decreasing above 100 TeV, and reaching minimal values at the highest energies. The multiplicity dependence of the simulated winter and summer multi-muon fluxes as a function of energy for θ = 30 • is shown in Figure 3, from which observes a higher flux of single-muon events in summer for E = 50 TeV (0.47 ± 0.26%) and E = 100 TeV (1.33 ± 0.40%). We find our results compatible with the observed sea-  sonal oscillation of a summer maximum of the singlemuon flux observed by the MINOS Near Detector [3] located at the same depth and location as NOvA-ND. For the multi-muon (N µ > 2) events, an excess of the flux in winter, relative to summer, is clearly visible. Accumulated distributions for the whole simulation library, shown in Figure 4, suggest a multi-muon flux excess in winter is increasing with N µ muon event multiplicity, as observed by the NOvA-ND [1].

Conclusions
We reported on our attempt to explain the winter multimuon flux excess observed by the NOvA and MINOS Near Detectors [1,2] as all previous works based on Monte Carlo simulations proved unsuccessful. In particular, in [5], a study using CORSIKA simulations has predicted an amplitude for the observed seasonal variation by a factor of four smaller than the one observed. That motivated us to use the Monte Carlo FLUKA package [8] since it allows for a more realistic modelling of the atmosphere while adequately treating the muon propagation in various media. We show that even considering the simple case of a pure proton nuclear mass composition for a set of fixed zenith angles and energies, we were able, for the first time, to describe the multi-muon flux excess in winter, as reported in [1,2]. Our results show that the winter excess seems to increase with the muon number N µ of the multimuon events. It rises from the lowest simulated shower energy of 10 TeV, reaching the maximum between 30 and 50 TeV, and decreases for 500 and 1000 TeV. A direct quantitative comparison of (winter-summer) differences in our calculations with the NOvA-ND result [1] is being evaluated and will be reported in a forthcoming work.