Reconstruction of the muon production longitudinal proﬁles in extensive air showers

. Muons produced in extensive air showers have large decay lengths and small radiative energy losses, enabling a large fraction of them to reach surface and underground detector arrays while keeping relevant information about the hadronic interactions that occurred high in the atmosphere. We can relate a muon’s arrival time and position at the detector to its production depth in the atmosphere. The total delay of muons with respect to the shower plane is primarily due to their geometric path and energy, we call these contributions the geometric and kinematic delays, respectively. We are working on the improvement of the current kinematic delay parameterizations using Deep Neural Networks for muons arriving at surface and underground detector arrays. We aim to reconstruct the longitudinal proﬁle of muons for future arrays of buried scintillator detectors at energies from around the second knee to the ankle of the cosmic ray spectrum, where there is an overlap with the nominal energies at the LHC. Given the low acceptance of scintillator detectors to inclined air showers and the richness of the forward physics near the shower core, we aim at applying a radial cut of 200 m instead of the usual 1000 m used in previous works.


Introduction
Cosmic rays with energies above 10 15 eV can only be detected indirectly through the secondary particles produced from their interaction in the Earth's atmosphere. Their energy and, particularly, nuclear mass composition are inferred by comparing the air-shower particles' time, lateral, and longitudinal distributions with those of Monte Carlo simulations using hadronic interaction models. However, the most relevant hadronic interactions of the most energetic cosmic rays occur at √ s ∼ 100 TeV, one order of magnitude above the nominal energy at the LHC, at phasespace regions not well-covered by most detectors. Lacking a theory to calculate multi-particle production from first principles, hadronic interaction models are our most significant source of systematic uncertainty [1].
Muons, mostly produced from the decay of charged pions, are the best messengers of the hadronic cascade. After production, they propagate in nearly straight lines to the ground without suffering significant energy losses, keeping relevant information about their mother particle [2][3][4]. In [5], from a compilation of the eight largest cosmic-ray experiments, a deficit of muons in simulations with respect to data is reported for E > 10 17 eV with a significance of more than 8σ, right where the center-of-mass energy equals the nominal energy at the LHC. Also, in [6], the Pierre Auger Collaboration observed a tension between the EPOS-LHC model predictions and data for the maximum depth of the muon production longitudinal profile.
In this work, we extend the method of the reconstruction of the muon production depth (MPD) described in [2,3], used in [6], for an array of underground muon detectors, as is the case of AMIGA [7], currently being deployed at the Pierre Auger Observatory.

Muon Production Depth
In [2,3], a method relating the arrival time of muons with their production distance is proposed. It is assumed that muons are produced in the shower axis and propagate to the ground in straight lines, as shown in Figure 1. Ac-  by an angle α, arriving at the ground with coordinates (r, ζ) expressed on the shower-front plane coordinate system. Here r is the muon distance from the shower axis, and ζ is its azimuthal angle. The total distance traveled by the muon is defined as l = r 2 + (z − Δ) 2 , where z is the muon production distance, measured along the shower axis, Δ = r cos ζ tan θ is the z-coordinate of the muon at the detector, expressed on the shower-front plane coordinate system, and θ is the shower zenith angle. Given the longer path traveled by muons, they will arrive delayed with respect to the arrival of the shower front. Assuming that muons are massless particles that travel at c, using a purely geometric relation, we can calculate the muon production distance, z, as where t g is the geometric delay of the muon. However, there are additional sources of delay, being the most relevant one the so-called kinematic delay, which is intrinsically associated with the energy of the muon [2][3][4]. Therefore, Equation (1) is rewritten as t g t − t , where t is the total delay, i.e., the one we can measure, and t is the kinematic delay. The latter is usually calculated from a parametrization assuming the average energy of muons arriving at a given distance from the shower axis, see Appendix B in [3]. Finally, the Muon Production Depth (MPD) is calculated from the muon production distance z as where ρ is the atmospheric density. The method proposed in [2,3] is intended for a detector array, as the one of the Pierre Auger Observatory [8], for θ ∼ 60 • showers, applying a radial cut of r > 1000 m. The reason for these cuts is two-fold; namely, they ensure negligible electromagnetic contamination in the total signal, while the contribution of the geometric delay, from which the MPD is calculated, is the dominant source of delay. In [6], this method was applied to the Pierre Auger data for events with E > 10 19.3 eV and 55 • ≤ θ ≤ 65 • , using detectors at 1700 < r [m] < 4000. Given the 1500 m spacing of the Auger detectors, the array proved too sparse, limiting the method's applicability to fewer than 500 events.
Our goal is to extend the present method's applicability to an array of buried scintillator detectors, as in the one of AMIGA [7]. Another advantage of AMIGA is its 750 m and 433 m detector spacing, allowing us to encompass the 10 17 eV energy region, where overlap with the LHC energies exists. That enables us to reconstruct the MPD in an energy region where we expect reduced systematic uncertainties from hadronic interaction models. However, since typical scintillator detectors have ∼ 5 cm width, their acceptance is limited to θ < 55 • , [9]. Also, the shower footprint is smaller at low zenith angles and energies, requiring a significantly lower radial cut of r 200 m.

Towards the reconstruction of the MPD
In this subsection, we report on our assessments of the feasibility of the extension of the present method to θ < 60 • , applying a radial cut of 200 m, as required for an array of buried scintillator detectors as it the case of AMIGA [7]. The AMIGA detectors are being deployed at 2.3 m depth, resulting on an additional mass overburden of 540 g cm −2 for a local soil density is ρ = 2.4 g cm −3 [7]. For our calculations we assume that all muons behave as minimum ionizing particles, i.e., losing a constant fraction of their energy per unit of crossed superficial density of a (E) = 1.8 MeV/g cm −2 [10], from which we calculate the minimum kinetic energy for the muons to reach the AMIGA detectors as where h = 2.3 m is the depth at which the detectors are buried, and θ is the zenith angle of the shower. From Equation (3), we expect that the minimum kinetic energy for a muon to arrive to the underground detectors to increase as a function of the zenith angle, resulting it to be E μ kin,min = 1 GeV (2 GeV) at θ = 0 • (θ = 60 • ). Throughout our work, we will use the results of Equation (3) to select the muons from which we reconstruct the MPD at 2.3 m depth.

CORSIKA shower libraries
For our studies, we produced two CORSIKA [11] shower libraries, one with a standard thinning of 10 −6 , combined with a radial thinning of r < 50 m, and a second one in which we did not activate the thinning algorithm. We set the ground level at 1452 m, the average height of the Pierre Auger Observatory. The second library was aimed to study the impact of thinning on the muon distributions close to the shower axis. In both cases, we use a fixed energy of 10 17 eV, and several zenith angles, assuming proton and iron-initiated showers. We have simulated 200 showers per nuclear species, zenith angle, and high-energy hadronic interaction model.

Testing the original model assumptions
In [2,3], it was assumed that muons are produced in the shower axis. In an update in [4], the average distance of the production of muons from the shower axis was estimated to be about a few tens of meters for 10 19 eV showers with θ = 60 • , using pre-LHC model predictions. We check if these findings hold for the current models and at lower-energy and zenith angles. In this section, for illustration, we show our results for 20 unthinned Sibyll 2.3d iron-initiated showers at 10 17 eV and θ = 0 • . In the text, we compare the results from the figures with the outcome of proton-initiated showers and also at θ = 60 • , when pertinent.
In Figure 2, we plot the distance from the shower axis at which the muons are produced. We observe that, for both primary species, at θ = 0 • , muons are produced at r prod ∼ 40 m, a distance that reduces to ∼ 30 m at θ = 60 • , confirming the findings of [4]. We also see that at the largest distances we are dominated by a population of low-energy muons. Quantitatively, at 0 • , muons with 1 < E i [GeV] < 2 (2 < E i [GeV] < 5) are produced at an average distance of 45 m (30 m) from the shower axis. In a second step, we studied the muon spectra and arrival time distributions of muons at the ground level as a function of the distance to the shower axis, as shown in the left and right panels of Figure 3, respectively. The arrival time of muons was calculated as t = t − t 0 , where t is the elapsed time since the first interaction, and t 0 = h 1 − h g / (c cos θ) is the time needed for a hypothetical massless particle to travel, parallel to the shower axis, from the first interaction point to the ground. In the equation, h 1 is the height of the first interaction, as given by CORSIKA, h g is the height of the ground level, and c is the speed of light. From the top panel of Figure 3, we observe that, near  the core, we have a very rich muon spectrum, ranging from hundreds of GeV to less than 1GeV. However, as we move further away from the shower axis, we become dominated by a less varied and less energetic muon component. We also find that the muon spectrum hardens with the zenith angle, most likely due to the decay of low-energy muons in the increased air overburden, which grows by a factor of two at θ = 60 • . We further observe that for θ = 0 • , there is a much richer energy spectrum at r < 200 m than the one for 60 • . On the bottom panel of Figure 3, we see that the arrival time distribution of muons is narrower closer to the shower axis, broadening as we move further away from it. There is a change in behavior in the arrival time distributions as a function of the zenith angle for the near and far regions from the shower axis. Namely, at 0 • , we have a broader arrival time distribution of muons arriving at r < 200 m, than at 60 • .
Finally, in Figure 4, ure 4, we observe that the most (least) energetic muons dominate the early (late) part of the distribution. However, at the beginning of the arrival time distribution, the muon energy spectrum is very diverse, ranging from the lowest to the highest energies, making the proper estimation of the kinematic delay challenging.

Estimating the kinematic delay
In [2,3], the estimated kinematic delay of muons traveling a distance l is estimated calculated as: where m 0 is the muon rest mass, a = 2 MeV/g cm −2 , assuming that muons behave as minimum ionizing particles [10]. E i , and E f are the muon energy at production, and at the detector, respectively. We will use Equation (4) to explicitly calculate the kinematic delay of muons t . In Figure 5, we plot the total (black), geometric (blue), and kinematic (red) delays of muons as a function from the shower axis. The estimate the total delay as the t = t − t 0 , as explained above, and the geometric delay as t g = t − t . From the middle and bottom panels of Figure 5, we confirm that the kinematic delay is intrinsically related with the energy of muons, and it is larger for low-energy muons. We also observe that the kinematic delay is dominant close to the shower axis, prevailing over the geometric delay until r = 200 m for E < 1 GeV.

Calculating the MPD using DNNs
For the MPD reconstruction, we used a Fully Connected Neural Network with 6 hidden layers in a pyramid scheme, where the number of nodes decreases from 1024 to 32 from the first to the last hidden layer. We use ReLU as the activation function and the ADAM optimizer. Our input features to the network are the sec θ of the shower, the ground coordinates of the muon (r, ζ), and its arrival time t , as explained above. We also add a boolean flag that tells us if the muons had enough energy to reach the underground detectors, according to Equation (3). Our loss function was derived empirically as where z DNN and z MC are the predicted and Monte Carlo values of muon production distance, respectively. ρ and γ 3 are the Pearson correlation coefficient and the distribution skewness, respectively, and c is a normalization constant. Finally, we convert the estimated muon production distance z into the MPD by applying Equation (2). We used 10% of the showers for the training of the DNN.
Our main results are shown in Figure 6. In the top panel, we show the correlation matrix between the reconstructed and Monte Carlo values of the MPD for all muons reaching the ground independently if they had enough energy to reach 2.3 m depth, from which observe a Pearson correlation of 0.936 between the reconstructed and true MPD. On the bottom panel, we show one MPD profile reconstructed at 2.3 m depth for one 10 17 eV proton-initiated shower at θ = 0 • , assuming a radial cut of 200 m. We see a good agreement between the reconstructed and Monte Carlo profiles up to a production depth of 500 g cm −2 . However, the DNN seems to perform worse for deeper MPD, probably due to the estimation of the kinematic delay since we expect that, at the deeper regions, the profile is dominated by low-energy muons.

Conclusions
We extended the MPD reconstruction method presented in [2,3] to an array of buried scintillator detectors. The applicability region of our work is intended for θ < 60 • , and r > 200 m, a fully complementary to the previous method. We were able to successfully reconstruct the MPD both for surface and underground muons using a Fully Connected DNN. As we expected, the estimation of the kinematic delay is the largest source of systematic uncertainty in the model. Our preliminary results show that, for the current DNN architecture, we achieve a muon-by-muon reconstruction of the MPD profiles at 2.3 m depth of ΔX < 10 g cm −2 , and σ (ΔX) ∼ 80 g cm −2 for E = 10 17 eV proton and iron-initiated showers, for 0 • < θ < 60 • , when a radial cut of r > 200 m is applied. These results can be improved if a different type of Neural Network or loss function is used. We plan a journal publication with updated results.