Probing the energy spectrum of hadrons in proton air interactions at ultrahigh energies through the fluctuations of the muon content of extensive air showers

We demonstrate that the shower-to-shower fluctuations of the muon content of extensive air showers correlate with the fluctuations of a variable of the first interaction of Ultra High Energy Cosmic Rays, which is computed from the fraction of energy carried by the hadrons that sustain the hadronic cascade. The influence of subsequent stages of the shower development is found to play a sub-dominant role. As a consequence, the shower-to-shower distribution of the muon content is a direct probe of the hadron energy spectrum of interactions beyond 100 TeV center of mass energies.


Introduction
Ultra-High-Energy Cosmic Rays (UHECRs) collide with nuclei of the Earth's atmosphere at center-of-mass energies that surpass the LHC energy scale, creating huge cascades of particles, the Extensive Air Showers (EAS) [1]. EAS constitute an unique experimental opportunity to explore hadronic interactions well beyond the energy scale attained by the largest human-made accelerator, the Large Hadron Collider (LHC).
On the other hand, the precise composition of UHECR is still unknown, as the nature of the primaries has to be inferred from EAS measurements themselves by using phenomenological models of the hadronic interactions tuned to a limited range in energy and kinematic phase-space. Breaking this degeneracy is one of the challenges currently faced by UHECR physics. This is to be achieved by increasing the number of independent observables taken from EAS.
Muons are direct messengers from the hadronic activity within EAS. The overall muon number can signal inconsistencies in the shower description. In fact, measurements at the Pierre Auger Observatory have shown that simulations, using the LHC-tuned hadronic interaction models, underestimate the muon content in air showers [2,3]. It has also been shown that the shower-to-shower variation in the muon content is use-Email address: friehn@lip.pt (Felix Riehn) ful information in discerning the composition of UHECR primary [4]. As a consequence, the effort to better measure EAS muons is increasing and several experiments are deploying upgrades to their detectors to enhance the sensitivity to them (see for instance [5]).
It is thus necessary to deepen research into the phenomenology of the muon component and in particular to understand the shower-to-shower distributions.
In this letter, we demonstrate that the fluctuations of the muon content of EAS directly correlate with the fluctuations of a variable that probes the hadron energy spectrum of the first interaction. The subsequent shower interactions give the overall scale of the muon content, but play a sub-dominant role in the shower-to-shower fluctuations. responding to the π 0 production rate. As the energy within the hadronic cascade decreases, it eventually becomes more likely that unstable hadrons decay rather than interact, this being the stage of the hadron cascade where most muons are formed [1,6]. In Fig. 1 the shower-to-shower distributions of the number of muons, N µ , for ensembles of simulated air showers initiated by different primaries are shown. The primary energy is 10 19 eV and the zenith angle is 67 • . We define N µ as the number of muons with E > 1 GeV arriving at ground. If one considers lower zenith angles the truncation of the muon production distribution and the muon propagation [7] change the resulting distributions, although their main characteristics remain the same (see Appendix A).
At 10 19 eV, air showers typically contain 10 7 muons with E > 1 GeV. In Fig. 1, one sees that the distribution differs substantially between different primaries. The relative fluctuations σ[N µ ]/ N µ for proton primaries are ∼ 0.15 (post-LHC model average 1 . In general, this value is much larger than 1/ N µ = 10 −4 , which is what one would obtain if the muons in the EAS were produced independently. Instead, the relative fluctuations indicate that the number of independent participants is 1/0.15 2 ∼ 44, a value of the same order as the number of particles typically produced in a single hadronic interaction. This suggests that the fluctuations could be dominated by phenomena in the first interaction. This idea has been tested with a modified air shower simulation wherein the first interaction has been fixed and showers were simulated by reinjecting all secondaries from this particular first interaction. The resulting effect on the distribution of the number of muons is shown in Fig. 1 by the dashed line. It can be seen that the width of the distribution is reduced while the average remains the same. The relative fluctuations in this modified simulation are much lower, around 0.05. A simple yet effective picture to understand this in terms of the shower dynamics is provided by Heitler-Matthews (HM) type models [8]. In such models, air showers are described as codeveloping EM and hadronic cascades, where all hadronic interactions yield a constant multiplicity. The produced particles are separated into two groups, those that decay into EM particles (π 0 ), feeding the EM cascade, and those that remain in the hadronic cascade (π ± ), accounting for the hadronic multiplicity. In each interaction, energy is equally shared by all produced particles. The variable describing the development of the shower is the generation number, i, which counts the number of interactions a particle has undergone. Particles arising from the first interaction are generation i = 1. Due to the equipartition of energy and to the transfer to the EM cascade, the energy in the hadronic cascade rapidly decreases. When it falls below a critical energy, E c , the cascade is terminated and the hadronic particles decay into muons. The average number of muons in the shower can be written as a function of the total (m tot ) and hadronic (m) multiplicity: Here g is the critical generation number, E is the energy of the primary particle, C = E −β c is a normalization constant, and β = ln m/ ln m tot . We have set β to 0.93, taken from post-LHC hadronic interaction models (Appendix B).
In order to understand the shower-to-shower fluctuations of N µ , one must consider that the multiplicity varies from interaction to interaction. The average hadronic multiplicity in generation i can be calculated as is the number of interacting hadrons in generation i (i − 1). The number of muons in the shower is then N µ = g i=1 m i . Note that, while the critical generation, g, is well defined for the HM model, allowing fluctuations of the hadronic multiplicity changes the overall energy budget of each sub-shower, leading to the fluctuation of the generation where the energy threshold E c is reached. Therefore, g must now be understood as an average parameter for the whole shower.
Assuming the hadronic multiplicities of all interactions arise from a common probability distribution with mean, m, and dispersion, σ(m), the fluctuations of m i are given by We thus find that the fluctuations of the average multiplicity in generation i are suppressed by the number of interactions resulting from the previous generation N i−1 . As the number of particles/interactions grows exponentially with the number of generations, the contributions to the fluctuations from later shower stages become increasingly smaller. This behaviour is typical for cascade processes with fixed multiplicity, as they occur in photomultiplier tubes [9,10]. In addition, for realistic multiplicity distributions, σ(m) decreases with energy [11]. As the interaction energy decreases from one generation to the next, fluctuations from later stages are further suppressed. We conclude that the overall fluctuations in the number of muons are dominated by the fluctuations in the first interaction, and that the contributions from further generations are exponentially suppressed.
So far, fluctuations of the shower were explained by fluctuations of the particle multiplicity. Nevertheless, an additional source of fluctuations comes from the fact that in each interaction, energy is shared among the emerging particle in a uneven and stochastic way. Taking into account only those fluctuations arising from the first interaction, one can write the total number of muons as the sum of the average number of muons that are produced in the m 1 subshowers that come out of the first interaction, i.e.
where the subscript 1 in N µ,1 denotes that only the fluctuations in the first interaction have been accounted for, and E j denotes the energy carried by the jth particle. We have used the fact that nucleon and pion initiated showers produce similar number of muons, as shown in Fig. 1. We define x j = E j /E as the fraction of the primary energy carried by particle j. Each sub-shower is thus weighted by x β j in the final number of muons. Defining we find that the number of muons in a shower is related to the average number of muons N µ,1 (E) = α 1 N µ (E) . For β = 1, the variable α 1 represents the fraction of energy that is passed on to the hadronic cascade. Its distribution is given by the hadronic energy spectrum [12]. In the opposite case (β = 0), α 1 becomes m 1 , so β shifts the weight between energy and multiplicity. In contrast to the HM model, both the multiplicity and energy fluctuations are included. We can also introduce the fluctuations induced by the second generation, writing N µ,2 (E) = m 1 j=1 jk where E jk = x j x jk E is the energy carried by second generation particles, where j and k run through the different combinations of 1st and 2nd generation particles respectively. This procedure can be generalized to account for any number of generations. Given α 1 we can define α i recursively for any generation as α i ≡ N µ,i /N µ,i−1 and translate the sum of sums into a product. Including fluctuations up to generation g the number of muons is given by As before, the number of particles in the hadronic cascade increases with generation so the fluctuations in α i decrease. By grouping fluctuations of all the generations but the first in a single parameter ω = N µ (E) g i=2 α i , the total number of muons can finally be written as

EAS Monte Carlo test
In order to verify our hypothesis, we have run 10 3 full 3Dsimulations of proton showers with Corsika (v7.6400) [13] at a primary energy of 10 19 eV and zenith angle of 67 • . The energy threshold for the thinning algorithm was set to thin = 10 −6 and weights were limited to w max = thin · (E 0 /1 GeV). We have also run 10 5 1D-simulations of air showers with Conex (v5.40) [14,15] for different primaries, energies 10 16 ,10 17 ,10 18 and 10 19 eV and zenith angles of 38 • and 67 • . High-energy interactions were simulated with the post-LHC models Epos-lhc [16], QGSjet II-04 [17] and Sibyll 2.3c [18,19], low-energies with Fluka v2011.2c [20]. The information about the first interactions was recorded [21]. Hadronic particles, i.e. particles that form the hadronic cascade and that determine α 1 , were defined as all hadrons other than π 0 and η. We refer to the number of muons at ground by N µ . In Fig. 2, the N µ and α 1 joint distribution, f α 1 ,N µ (α 1 , N µ ) derived with Conex is shown with events from Corsika superimposed. The correlation coefficient is approximately 0.8 2 . For comparison, the correlation with the fraction of hadronic energy and with the multiplicity of the first interaction is about 0.6 and 0.2 respectively. A detailed comparison of the correlations between variables of the first interaction and the number of muons is given in Appendix C. The correlation between N µ and the fraction of hadronic energy is discussed in Appendix D. Looking at the joint distribution shown in Fig. 2 one can see that the bulk of events follows a linear relation of the form: While approximately 1/3 of pions that are produced in a pp interaction are neutral, due to the preservation of the quantum numbers of the proton (leading particle effect) the fraction of events where most of the energy is carried by a neutral pion is much smaller. In these rare cases, α 1 → 0, and muon production is dominated by photo-pion production in the EM cascade. The additional term δ in Eq. (6) accounts for this (compare to Eq. (5)). Its value was found to be roughly 0.16 across different models (Appendix B). The distribution of the number of muons produced in a photon-initiated shower is shown in Fig. 1  vertical structure seen around α 1 = 1 in Fig. 2 corresponds to quasi-elastic/diffractive events, where a small amount of energy is transferred from the incoming primary to the target nucleus.
For convenience, we define n ≡ ln N µ , a ≡ ln(α 1 + δ) and w ≡ ln ω, and thus have n = a + w. The probability distribution of f n (n) is then given by f n (n) = f a (a) f w (n − a |a) da, where f w (w |a) is the w-distributions at different a, and can be accessed in simulations (see Appendix F). By neglecting the a dependence of f w (w |a), we arrive to the approximation f n (n) f a (a) f w (n − a) da . In Fig. 3, f a (a) and f w (w) are shown, where f w (w) ≡ f w (w |a) f (a) da was called average approach. Fig. 4 compares the average convolution approach with the exact distribution of a proton primary for different models. It can be observed that the general features are well reproduced: a) The low-n tails, a reflection of the different low-a tails in the models; b) The width σ(n) = 0.17, which is dominated by the width σ(a) = 0.14 as compared with σ(w) 0.11. Hence, w plays a sub dominant role in the N µ fluctuations (notice that f w (w) includes all contributions to the shower-to-shower fluctuations which are not accounted for by f a (a)); And finally c), the different average n-values are a reflection of the different average values of the w-distribution in Fig. 3. The high-n tail nevertheless would need the full correlation convolution to achieve high precision, which is out of the scope of this paper. Fig.  5, compares also a Gaussian convolution approach, where the average and width were extracted from a Gaussian fit to f w (w).
For a UHECR observatory to access the a-distribution of proton-air interactions experimentally, the shower effects would need to be unfolded from the observed n-distribution. This requires the use of a f w (w|a) distribution, where the a − w correlations should be accounted for if one wants to attain the max- imum precision. A contribution to the systematics comes from the differences in f w (w|a) across models. The present paper has laid the basis for these future analyses by demonstrating that the n-distribution is shaped by a more fundamental distribution of multiparticle production in the first interaction and therefore can be experimentally probed. The physics of the a-distribution itself and its relation with accelerator measurements is another matter for future research.

The scenario of nucleus-air interactions
In general, the UHECR composition may be composed of different types of nuclei. The shower-to-shower fluctuations of the number of muons also reflect the fluctuations of the primary UHECR mass on a shower-to-shower basis. In Fig. 5 a comparison between MC simulations for different primaries with the mass number A (number of nucleons) and the results of the superposition model are shown for the hadronic interaction model Epos-lhc . The n-distributions for nuclei were built from the Afold self-convolution of the nucleon n-distribution with energy E/A. The gradual change in the shape of the distribution in the transition from proton showers to heavier primaries is described well by the superposition model. Fluctuations are slightly underestimated, which is expected since superposition neglects, for example, the fluctuations of the energy between the nucleon sub-showers [22]. For other hadronic interaction models the results are similar (see Appendix E). Fig. 6 displays the joint-distribution of the number of muons and the modified hadronic energy fraction α 1 after the first interaction for a mixture of different primaries. Heavier nuclei display a larger α 1 value, reflecting the fact that, only a portion of the nucleus interacts, leaving the rest of the energy in the remaining nucleons. The number of muons increases in the same proportion, and as a consequence the populations separate. Fig. 7 shows the w-distribution and a-distribution. The wdistribution has a maximum at the same position, although the width becomes narrower as the mass number increases, while the a-distribution reflects the self-folding of A nucleons according to the superposition model and thus is shifted with the mass number. These figures are based on simulations with Epos-lhc . For other models the general behaviour of a shift of the populations along the diagonal in the N µ -α 1 plain with increasing mass remains the same, but details change for heavier nuclei (A > 4) (see Appendix E). While these details of evolution with mass may be interesting in their own right and deserve study elsewhere, they are outside the scope of this paper. The essential point here is that the tails in n and a are suppressed for nuclei, so that even for a mixed composition the signature tail in protons remains visible ( Fig. 7 and Fig. 5). While there is evidence for a mixed mass composition of cosmic rays in the energy range between 10 18.5 eV and 10 19 eV [23], discussion still remains over the value of the average mass reported by different experiments [24,25]. Nevertheless, by using a mass sensitive shower observable (like the depth at which the shower reaches its maximum developement, X max ) one could experimentally select samples of events which are proton enriched, and then follow the same procedure as for the case of a pure proton UHECR beam. The unfolding of the n-distribution would result in a a-distribution that reflects the composition of the sample. If protons are present, the adistribution is dominated by protons at low-a values, making it possible to extract direct information of the first interaction. The same strategy was used to measure the proton-Air cross section by using X max [26,27].

Conclusions
In this letter, it has been demonstrated that the number of muons in EAS is connected with a variable of the first interaction, α 1 . This variable is computed from the fraction of en- ergy carried by the hadronic particles that sustain the hadronic cascade, being a direct probe of the hadron energy spectrum. Using this knowledge, it was shown that the shower-to-shower distribution of the muon content can be explained by the α 1distribution of the first interaction of the UHECR. The subsequent shower interactions give the overall scale of the muon content, but play a sub dominant role in the fluctuations. It is worth noting that the p-Air cross section represents the only other example of a direct link between EAS measurements to a property of the first interaction of the UHECR, and which has been measured [26,27]. The present work demonstrates for the first time that it is possible to directly access multiparticle production variables of the p-Air interaction, which can occur at center-of-mass energies above the LHC scale. [3] A. Aab The zenith angle dependence of the correlation between the number of muons at ground and at production with α 1 is shown in Tab. A.1. In both cases, at ground and at production, there is only a small difference in the correlation between vertical and inclined showers. The effect of the zenith angle on the distribution of the number of muons is shown in Fig. A.8.  In general we find that the results for proton showers are very similar between the interaction models. In Tab. B.2 the first two moments of α 1 and N µ are compared between air shower simulations for proton primaries with the different post-LHC interaction models. The slope of the average number of muons β and the contribution to muon production from the EM cascade, δ, are shown as well. In Fig. B.9 the correlation between N µ and α 1 in proton showers is shown for the interaction models Sibyll 2.3c and QGSjet II-04 (Epos-lhc is shown in Fig. 2).  Defining the correlation coefficient ρ X,Y = cov(X, Y)/(σ X σ Y ), Tab. C.3 shows the correlation of different variables of the first interactions with the number of muons at ground and at production. The statistical uncertainty of the correlation coefficient is given by (1 − ρ 2 )/(N − 2) and is in the range of 0.02 − 0.04 accross variables and interaction models. The sample size was 10 3 and the simulations were done with Corsika .
The variables presented are: the hadronic multiplicity in the first interaction m 1 , the total multiplicity m tot and the fraction of energy carried by all hadrons E had /E. A variable mixing information from multiplicity and the hadronic energy fraction is realized by α 1 . It is defined as where β is the slope of the energy dependence of the average number of muons (see Tab. B.2) and x i is the fraction of energy E i /E carried by hadron i. The fraction of energy carried by the leading hadron relative to the total energy in the hadronic cascade E had is denoted by * . The inelasticity, i.e. the fraction of the primary energy that is carried by all other hadrons except the leading one is κ inel . The depth of the first interaction is X 0 .
Appendix D. Shower-to-shower fluctuations in the muon content: comparison of α 1 and E had /E The modified energy fraction α 1 contains more information on particle production and therefore exhibits a stronger correlation with muon production than the pure energy fraction of hadrons E had /E = m i x i . On the other hand E had /E is strictly constrained to the interval between 0 and 1 by energy conservation, while α 1 is only constrained by the multiplicity. Given a maximal multiplicity m max and assuming energy is shared equally (x i = 1/m max ), the maximum of α 1 is

Appendix E. The scenario of nucleus-air interactions: hadronic interaction models
In Tab. E.4 the evolution of the parameters of muon production for the hadronic interaction models Epos-lhc , Sibyll 2.3c and QGSjet II-04 as a function of the primary mass is shown. As the superposition model predicts the average number of muons increases with the mass while the relative fluctuations decrease. This behaviour is confirmed by Fig. E.12 (Fig. 5 for Epos-lhc ) which shows the full distribution of N µ for air showers at 10 19 eV. For α 1 the mean shows a similar monotonous increase. The fluctuations of α 1 on the other hand first decrease with the number of nucleons and then increase again. In Fig. E.13 (and Fig. 7 for Epos-lhc ) the evolution of the distributions of a = ln α 1 and w = ln ω with mass is shown. Since α 1 depends on energy and multiplicity, it is sensitive to the frag-