On the model uncertainties for the predicted muon content of extensive air showers

Motivated by the excess of the muon content of cosmic ray induced extensive air showers (EAS), relative to EAS modeling, observed by the Pierre Auger Observatory, and by the tension between Auger data and air shower simulations on the maximal muon production depth $X^{\mu}_{\max}$, we investigate the possibility to modify the corresponding EAS simulation results, within the Standard Model of particle physics. We start by specifying the kinematic range for secondary hadron production, which is of relevance for such predictions. We further investigate the impact on the predicted EAS muon number and on $X^{\mu}_{\max}$ of various modifications of the treatment of hadronic interactions, in the framework of the QGSJET-III model, in particular the model calibration to accelerator data, the amount of the"glue"in the pion, and the energy dependence of the pion exchange process. None of the considered modifications of the model allowed us to enhance the EAS muon content by more than 10\%. On the other hand, for the maximal muon production depth, some of the studied modifications of particle production give rise up to $\sim 10$ g/cm$^2$ larger $X^{\mu}_{\max}$ values, which increases the difference with Auger observations.


Introduction
Experimental studies of very high energy cosmic rays (CRs) are traditionally performed using extensive air shower (EAS) techniques: measuring various characteristics of nuclear-electromagnetic cascades initiated by interactions of primary CR particles in the atmosphere [1].Such investigations necessarily involve an extensive use of numerical simulation tools designed to describe the EAS development [2].An accurate reconstruction of the properties of CR primaries depends thus strongly on the precision of air shower modeling.This is especially so for CR composition studies since the relevant EAS characteristics are governed by interactions with air nuclei both of primary CR particles and of secondary hadrons produced in the shower cascade process.Consequently, a successful determination of the primary CR composition is only possible if such hadronic interactions are correctly described by the corresponding Monte Carlo (MC) generators.
One of the traditional methods for CR composition studies relies on measurements of the muon component of air showers by ground-based detectors.However, for ultra-high energy cosmic rays (UHECRs), the use of that technique remained hampered for nearly two decades by the persistent contradiction between experimental observations and EAS simulations, regarding the muon number at observation levels [3,4,5,6], the discrepancy reaching 50% level.Using present air shower simulation tools, the observed muon number N µ would correspond to primary particles being as heavy as uranium nuclei, which, given the exceedingly rare natural abundance of such nuclei, appears very unlikely.A number of theoretical approaches potentially capable to resolve this so-called "muon puzzle" [7] has been proposed, typically advocating to physics beyond the Standard Model [8,9,10,11].However, since such suggestions are of a rather speculative character, it may be important to perform a quantitative study of the uncertainty range for the predicted EAS muon content, within the standard physics picture.
Such an investigation is the subject of the current work.Using the recently developed QGSJET-III model [12,13], we study the dependence of the predicted N µ on uncertainties regarding the model calibration to relevant accelerator data and on possible changes of the underlying theoretical mechanisms.In addition to the total muon number at ground level, we consider also the maximal muon production depth X µ max measured by the Pierre Auger Observatory [14,15], for which significant discrepancies with predictions of EAS simulations have also been observed.Compared to other studies of uncertainties of EAS predictions (e.g.[16,17]), our investigation differs in three key aspects: i) the changes of the corresponding modeling are performed at a microscopic level; ii) the considered modifications are restricted by the requirement not to contradict basic physics principles; iii) the consequences of such changes, regarding a potential (dis)agreement with relevant accelerator data, are analyzed.
The paper is organized as follows.In Section 2, we study the kinematic range for secondary hadron production, which is of relevance to N µ and X µ max predictions, using the QGSJET-II-04 [18,19], EPOS-LHC [20], and SIBYLL-2.3[21] MC generators.In Section 3, we consider potential variations both of the parameters and of the underlying theoretical mechanisms of the QGSJET-III model, while checking the corresponding impact on the predicted N µ and X µ max .Finally, we discuss our results and conclude in Section 4.
2 Kinematics of secondary hadron production, relevant to N µ and X µ max predictions The muon component of extensive air showers is formed by a long nuclear cascade in the atmosphere, primarily driven by pion-air interactions. 1It is thus evident that the multiplicity of secondary pions or, more generally, of all kinds of relatively stable hadrons (i.e., those whose interaction length in the atmosphere is substantially shorter than the decay length), has a direct impact on the expected muon signal at ground, which can be qualitatively understood using the simple Heitler's picture of the cascade process [22] (see also [23]).A more delicate question addressed, e.g., in [24] (see also [25] for a recent study) is what range of energy fractions taken by secondary pions from their parent hadrons, x E = E/E 0 , is of highest importance for the predicted N µ .On the one side, secondary hadron production in high energy collisions is spread over the whole available rapidity range, 0 < y < ln s, in the laboratory (lab.)frame, s being the center-of-mass (c.m.) energy squared for the collision.Hence, the bulk of secondaries is characterized by vanishingly small x E ∝ e y /s.On the other hand, interactions of most energetic secondary hadrons initiate powerful enough subcascades, thereby providing a higher contribution to the observed N µ .The competition between these two trends can be easily quantified considering, for the moment, only contributions to muon production from interactions of secondary pions with air and making a rather reasonable assumption, being again motivated by Heitler's picture, that the energy-dependence of the muon content of a nuclear subcascade initiated by a charged pion can be approximated by a power law: Under these assumptions, N µ of a proton-initiated air shower is given as with dn π ± p−air (E 0 , x E )/dx E being the x E distribution of secondary charged pions in a proton-air collision at lab. energy E 0 .Such a distribution, for hadron h collisions with air nuclei, can be approximated as within the differences between the predictions of current MC generators, as illustrated in Fig. 1 for the case of a proton-nitrogen (pN) collision at E 0 = 10 17 eV.pN /dx E , at E 0 = 10 17 eV, as calculated using the QGSJET-II-04 (solid line) and SIBYLL-2.3(dash-dotted line) models.Shown by the dotted line is the approximation of the SIBYLL-2.3results, using the ansatz of Eq. ( 3), with ∆ = 0.4 and β = 4.5.
Inserting Eq. (3) into Eq.( 2) and making use of Eq. (1), we have immediately From Eq. (4), we get the "average" x π E for p-air interactions, weighted by the contributions of the corresponding pion subcascades to N µ : For α µ ≃ 0.9, using the parameters of the fit in Fig. 1, ∆ = 0.4 and β = 4.5, we get x π E ≃ 0.08.In Fig. 2, we plot the cumulative x E distribution, for these parameter values.As one can see in the Figure, only 1/4 of the contribution comes from the interval x E < 0.01.The same reasoning can be repeated for pion-air collisions, if we neglect the diffractive contribution which, as will be demonstrated in the following, is of minor importance for N µ , yielding a similarly large value for x π E .This is illustrated in Fig. 3, where the x E distribution of charged pions for π + N collisions at E 0 = 1 PeV is approximated by the ansatz of Eq. (3), with ∆ = 0.3 and β = 3.5.Using these values of ∆ and β in Eq. ( 5) yields x π E ≃ 0.12.Thus, the air shower muon content is governed by forward pion production 2 in hadron-air collisions, at the many steps of the nuclear cascade in the atmosphere.Let us further support this  6), for ∆ = 0.4, β = 4.5, and α µ = 0.9.x E distribution of charged pions for π + N collision, x E dn π ± π + N /dx E , at E 0 = 1 PeV, as calculated using the QGSJET-II-04 (solid line) and SIBYLL-2.3(dash-dotted line) models.Shown by the dotted line is the approximation of the SIBYLL-2.3results, using the ansatz of Eq. (3), with ∆ = 0.3 and β = 3.5.
conclusion by a numerical study, using a number of popular MC generators of CR interactions.To this end, we modify pion energy distributions provided by those models, replacing every second pion characterized by energy fraction x E bigger than some cutoff x cut by a pair of pions of half the energy.Importantly, doing so, we replace a charged pion by a pair of charged pions, while a π 0 is being "split" into a pair of neutral pions, such that the energy partition between nuclear and electromagnetic (e/m) cascades remains unaltered.
The modifications of secondary pion energy distributions by such a procedure are exemplified in Fig. 4. For large x cut ≃ 1, this impacts the diffractive peak only, leaving a half of it; such a peak reappears then at x E ≃ 0.5.Choosing a smaller x cut , a larger and larger part of the very forward pion spectrum is reduced, being accompanied by an enhancement of pion production at smaller x E .In the limit x cut → 0, the forward pion production is reduced considerably, while the total multiplicity of secondary pions is increased by 50%.
In Fig. 5, we plot the x cut dependence of the relative change of the calculated N µ (E µ > 1 GeV) at sea level, ∆N µ /N µ , for p-induced EAS of E 0 = 10 19 eV, applying such a "splitting" procedure PeV, as calculated using the QGSJET-II-04 model (solid line).The modifications of this distribution by the "splitting" procedure discussed in the text are shown by the dotted, dashed, and dash-dotted lines for x cut = 0.95, 0.3, and 0.1, respectively.to all hadron-air interactions in the shower, 3 which are treated by the QGSJET-II-04, EPOS-LHC, or SIBYLL-2.3models.The first thing to notice is that for x cut values close to unity, there is practically no change of N µ , i.e., pion diffraction is rather irrelevant for EAS muon content.Secondly, the largest gradient of the ∆N µ (x cut ) dependence is observed for x cut ∼ 0.1, thereby supporting our previous estimation, Eq. (5).In the limit x cut → 0, the EAS muon content is enhanced by ∼ 15 − 20%, depending on the interaction model.It should be stressed, however, that the performed study serves illustrative purposes only, not being a viable procedure for modifying the predicted N µ .For example, when the discussed "splitting" is applied to half the secondary pions (the case x cut → 0), one arrives to a strong contradiction with available accelerator data.This is illustrated in Fig. 6, where we plot momentum distributions of charged pions, for π − C collisions at 350 GeV/c, obtained using the QGSJET-II-04 model and applying the "splitting" procedure to half the secondary pions (the case x cut = 0), in comparison to NA61 data [27].
The same kind of analysis has been performed for the maximal muon production depth X µ max , the obtained ∆X µ max (x cut ) dependence being plotted in Fig. 7 x cut ∆X µ max ( g/cm 2 ) p-induced EAS (E 0 =10 19 eV) Figure 7: x cut dependence for the modification of X µ max by the pion "splitting" procedure, for p-induced EAS of E 0 = 10 19 eV, for different interaction models.The meaning of the lines is the same as in Fig. 5. 10 19 eV.Somewhat surprisingly, for this dependence, the maximal gradient is observed in almost the same x cut range: of maximal importance for X µ max are secondary pions carrying few percent of the parent hadron energy.As discussed in [28], X µ max roughly corresponds to the depth in the atmosphere, where the interaction and decay lengths for most of the pions in the cascade become comparable, i.e., where the bulk of the pions approaches the pion "critical" energy.Hence, this quantity is sensitive both to the pion-air inelastic cross section and to forward production spectrum of pions.What was not properly realised in that work is that the position of that "equilibrium" point can be changed most efficiently by modifying how fast the energy partition between the produced secondary hadrons is established, thereby speeding up the decrease of the average pion energy, with increasing depth, which can again be understood using the simple Heitler's picture for a nuclear cascade.Here, similarly to the case of N µ , one has a competition between a copious hadron production for x E → 0 and much less numerous secondaries characterised by finite x E , which, however, have a much stronger impact on the speed of decrease of the average pion energy.It is this competition which leads to the obtained range of pion energy fractions, corresponding to the maximal relevance for X µ max predictions.As a side remark, pion diffraction appears to have practically no impact on the predicted X µ max , contrary to what had been claimed in the literature [29] (cf.Fig. 7 for x cut → 1).
For completeness, we plot in Fig. 8  x cut Figure 8: x cut dependence for the modification of X max by the pion "splitting" procedure, for p-induced EAS of E 0 = 10 19 eV, for different interaction models.The meaning of the lines is the same as in Fig. 5.
maximum depth X max , under such a "splitting" procedure.As one can see in the Figure, the considered modifications of secondary pion energy spectrum have a very weak impact on the predicted X max , i.e., the treatment of pion-air interactions is of small importance for X max predictions: even in the extreme case, x cut → 0, the EAS maximum depth changes by less than 5 g/cm 2 .
To summarize this part of our analysis, both N µ and X µ max are governed by forward pion production in hadron-air collisions: with the corresponding pion energy fractions ranging between few per cent and few tens per cent.Obviously, these conclusions can be generalized by taking into consideration the production of all "stable" hadrons, i.e., the obtained ranges of energy fractions x E relevant for N µ and X µ max predictions apply to the production spectra of all such hadrons.To further support these conclusions, we plot in Fig. 9 the energy-dependence of the multiplicity of "stable" hadrons [(anti)nucleons, kaons, and charged pions] n πN stable and of the second moment of the corresponding energy fraction distribution,4 dn πN stable /dx E , i.e., the average energy fraction taken by all stable secondary hadrons, for pion-nitrogen collisions, SIBYLL-2.3predict noticeably higher values of x E n πN stable -due to a much more abundant forward production of, respectively, (anti)nucleons and ρ mesons in those models. 5The more abundant forward production of stable hadrons in the latter two models clearly explains their somewhat higher values of the predicted N µ , plotted in Fig. 10.Thus, it is the total fraction of the incident hadron energy, taken by all stable secondary hadrons, rather than the multiplicity, which governs the muon content of extensive air showers.
3 Potential variations of model predictions for N µ and X µ max

Enhancing the gluon content of the pion
As is evident from the analysis in Section 2, to increase the predicted N µ , one has to enhance forward production of stable hadrons.Generally, there exists no viable theoretical mechanism to produce a "hardening" of secondary hadron spectra, with increasing energy: the energy-rise of multiple scattering should rather give rise to a (rather moderate) opposite effect, leading to the energy partition between larger and larger numbers of secondaries.Therefore, to obtain an enhancement of forward production, one may try to increase the multiple scattering rate in pion-air collisions, desirably, without a significant modification of the shape of forward spectra: in order to have higher secondary particle yields for all values of x E .Since the energy-rise of multiple scattering is primarily driven by (semi)hard production processes giving rise to an emission of hadron (mini)jets, the only efficient and, probably, the only possible way to reach such a goal, without modifying substantially the treatment of proton-air collisions, rather seriously constrained by available data from the Large Hadron Collider (LHC), is to enhance the gluon content of the pion.
In principle, the fraction of the pion momentum, carried by its valence quarks, x qv , is already strongly constrained (e.g.[31]), notably, by experimental studies of the Drell-Yan and direct photon production processes in pion-proton interactions.Therefore, one may only change the partition of the remaining momentum fraction, 1 − x qv , between gluons and sea (anti)quarks, since experimental data on pion structure functions are rather scarce and correspond to relatively large values of Bjorken x.Shifting that momentum balance towards gluons, at the expense of sea (anti)quarks, would enhance the multiple scattering rate: because gluons participate in hard scattering with a larger color factor, C A = 3, compared to C F = 4/3 of (anti)quarks.
Here we wish to study an extreme variation of the glue in the pion.Therefore, we choose to modify the overall momentum fraction possessed by both gluons and sea (anti)quarks, x g + x qs , at the expense of x qv .I.e., for the moment, we disregard the respective experimental constraints and consider a factor of two reduction of x qv , 6 while enhancing correspondingly both x g and x qs .This way we change the gluon PDF rather radically in the x-range relevant for N µ and X µ max predictions, as demonstrated in Fig. 11.Yet, as one can see in Fig. 12 (left), such a modification gives rise to a noticeably higher multiplicity of stable secondary hadrons n πN stable , for pion-nitrogen collisions, only at the highest energies considered.Moreover, the average fraction of the parent pion energy, taken by all such hadrons, x E n πN stable , remains practically unaffected, see Fig. 12 (right), which means that the considered modification has a minor impact on forward hadron spectra.Not surprisingly, applying the so-modified treatment of pion-air interactions to EAS modeling, we observed rather insignificant changes for the calculated EAS muon content and for the maximal muon production depth, compared to QGSJET-III predictions ( 1% for N µ and few g/cm 2 for X µ max ).

Modifying the model calibration to experimental data
An alternative way to increase the predicted N µ is to enhance the yields of stable secondary hadrons, at the expense of neutral pions.In relation to that, let us pay attention to the fact that the QGSJET-III model seriously underestimates the production of kaons and (anti)nucleons in π − C collisions, compared to measurements by the NA61 experiment [27], as one can see in Figs. 13 and 6 This is achieved by modifying the original GRS parton distribution functions (PDFs) of valence quarks in the pion [31], employed in the QGSJET-III model, at the cutoff scale Q 2 0 = 2 GeV 2 for hard processes, by a factor A δ (1 − x) δ , and choosing the parameters δ and A δ in such a way that the desirable reduction of xq v is obtained without violating the valence quark sum rule.Additionally, we increase the value of the parameter βπ governing the "hardness" of the gluon PDF in QGSJET-III (∝ (1 − x) βπ -cf.Eq. ( 8) in [12]), from the default value 2 to 5, in order to shift the enhancement of the glue towards smaller x.PDFs at higher virtuality scales q 2 are obtained via a DGLAP evolution from Q 2 0 to q 2 .Figure 12: Energy-dependence of the multiplicity of "stable" hadrons (left) and of the energy fraction taken by all such hadrons (right) for π + N interactions, for the default QGSJET-III model (solid line) and for the enhanced gluon content of the pion, as discussed in the text (dashed line).
14 (solid lines).The agreement with the data can be improved performing suitable modifications of the hadronization procedure of the model: considering, respectively, 40% increase for the probability to create strange quark-antiquark (ss) pairs from the vacuum or 60% enhancement of the corresponding probability for diquark-antidiquark (udū d) pairs, when treating the hadronization of strings of color field, stretched between final partons produced [13], the corresponding results being shown in Figs. 13 and 14 by dashed lines.On the other hand, the agreement with NA61 data on ρ 0 meson production [32] can be further improved considering 50% enhancement for the probability of the fragmentation of light (anti)quarks into ρ mesons, compared to the default value 1/3 used in the model, see Fig. 15.
It is noteworthy, however, that the considered enhancements of kaon and (anti)nucleon production are at tension with the corresponding data on pion-proton collisions from other experiments, as one can see, e.g., in Fig. 16.Moreover, since the hadronization procedure is universal for all kinds of hadronic interactions, such changes lead to a strong contradiction with the corresponding particle yields measured in pp interactions, both at fixed target energies and at LHC (see [13]  for a comparison of the default QGSJET-III results with such data).In turn, for the considered enhancement of ρ meson production, the fraction of directly produced pions falls down to ≃ 30%, in a clear contradiction to observations [34].
Regarding the impact on EAS muon content, the considered enhancements of kaon and (anti)nucleon production give rise up to ∼ 10% increase of the predicted N µ , compared to the original results of QGSJET-III, as one can see in Fig. 17 (left).On the other hand, the considered increase of ρ meson yield appears to have only a minor effect: ∆N µ /N µ < 1%, which is a direct consequence of the isospin symmetry. 7Indeed, since such an enhancement of ρ meson production does not modify the energy partition between nuclear and e/m subcascades initiated by the, respectively, charged and neutral pions resulting from decays of ρ mesons, we are essentially back to the picture discussed in Section 2.Then, since this enhancement of ρ meson yield mainly affects central (in c.m. frame) production characterized by small energy fractions x E of the mesons, at sufficiently high energies, it is of minor importance for N µ predictions, as demonstrated in Section 2.
The effect of the considered modifications of particle production on the calculated maximal muon production depth X µ max is shown in Fig. 17  Figure 14: Momentum distributions (in lab.frame) of protons (top panels) and antiprotons (bottom panels) produced in π − C collisions at 158 GeV/c (left) and 350 GeV/c (right), as calculated using the default QGSJET-III model (solid lines) or considering an enhancement of (anti)nucleon production, as discussed in the text (dashed lines), compared to NA61 data [27] (points).and (anti)nucleons, one shifts the energy balance between hadronic and e/m subcascades in favor of the former, thereby slowing down the energy "leak" from the hadronic component of extensive air showers and elongating the nuclear cascade, as discussed in some detail in [28].The corresponding values of X µ max exceed the ones of the original QGSJET-III by up to ∼ 10 g/cm 2 .On the other hand, the considered increase of ρ meson production has a much weaker impact on X µ max .

Changing the energy-dependence of the pion exchange process
The very last remaining option regarding a potentially strong enhancement of EAS muon content is related to pion exchange process in pion-air collisions, i.e., to the so-called Reggeon-Reggeon-Pomeron (RRP) contribution, with R = π, to secondary hadron production.As discussed in [19], assuming a dominance of the pion exchange in the RRP configuration, due to the small pion mass, compared to other Reggeons, one obtains ≃ 20% enhancement of N µ , relative to the case when this contribution is neglected.This is because such a process changes the energy balance between secondary charged and neutral pions in favor of the former.Indeed, the dominant decay channels  considering an enhancement of (anti)nucleon production, as discussed in the text (dashed lines), compared to LEBC-EHS data [33] (points).
ρ ± → π ± π 0 and ρ 0 → π + π − give rise to for the products of forward ρ meson decays, in contrast to x π ± : x π 0 = 2:1 for the usual ("central") ρ meson production discussed in Section 3.2.However, the energy-dependence of that contribution depends crucially on the treatment of absorptive effects defining the probability for the corresponding "rapidity gap survival", i.e., that the (small) rapidity gap between the leading hadron (here, ρ meson) and other secondary particles produced by an interaction of the virtual pion with the target nucleus, is not filled by particle production due to additional multiple scattering processes [19,30,35,36,37]. of the modification of the maximal muon production depth X µ max (right), both for E µ > 1 GeV, for proton-initiated EAS, with respect to the corresponding predictions of the default QGSJET-III model, for the considered modifications of the model: enhancement of (anti)nucleon, kaon, and ρ meson production -solid, dashed, and dash-dotted lines, respectively.
Here we are going to consider an extreme case: neglecting such absorptive corrections and assuming that the pion exchange process constitutes a fixed fraction w lim π of the inelastic cross section.I.e., for any inelastic pion-air collision, we consider with the probability w lim π a production of a charged or neutral ρ meson, with equal probabilities (e.g.[38]), sampling its light cone momentum fraction x and transverse momentum p t according to the corresponding probability f π−air (x, p 2 t ) for reggeized pion exchange in the Born approximation (e.g.[39]), while treating the rest of hadron production for the event as an interaction of the virtual pion with air nucleus.Here with the squared momentum transfer t being related to p t as is the pion Regge trajectory with the slope α ′ π ≃ 0.9 GeV −2 , σ tot π−air is the total pion-air cross section,8 F (t) is the pion emission form factor,9 m π and m ρ are the pion and ρ meson masses, respectively, and s is the c.m. energy for the interaction.
Choosing w lim π =0.11, based on the data of the NA61 experiment [32], see Fig. 18, the resulting spectrum of ρ 0 mesons and the partial contribution of the pion exchange process for pion-nitrogen collisions at 1 PeV are shown in Fig. 19, in comparison to the corresponding predictions of QGSJET-III.In turn, the obtained energy-dependence of the relative change of N µ , corresponding to this alternative treatment, is plotted in Fig. 20 (left).
Somewhat surprisingly, such a treatment results in a lower muon number, compared to the predictions of the default QGSJET-III model.To understand this nontrivial result, let us remind ourselves that a considerable part of the cross section σ tot π−air in Eq. ( 9) corresponds to elastic scattering (see Fig. 21), with the corresponding fraction σ el π−air /σ tot π−air being ∼ 30% at fixed target energies and slowly rising with E 0 , approaching the "black disk" limit (σ el π−air /σ tot π−air → 0.5) in the E 0 → ∞ limit.For that elastic scattering contribution, the final state consists of a ρ meson and a pion only: with the virtual pion being put on-shell by the scattering process, see Fig. 21 (right).It is such a scarce hadron production which gives rise to the obtained decrease of N µ .
In contrast, in the treatment of the QGSJET-III model, the rise of absorptive corrections "pushes" the pion exchange process towards larger and larger impact parameters, causing a (slow) decrease of its relative contribution [30,37].At the same time, at large impact parameters, the above-discussed contribution of elastic scattering of the virtual pion is strongly suppressed and can be neglected [37].
Regarding the maximal muon production depth, in addition to a scarce hadron production Schematic view of the two contributions to the pion exchange process: for inelastic scattering of the virtual pion (left), the leading ρ meson is accompanied by multiple hadron production; elastic scattering (right) gives rise to only one secondary pion, in addition to the ρ meson.The light shaded croissant and ellipsis correspond to contributions of elastic rescattering processes.in case of elastic scattering of the virtual pion, one has a shift of the energy balance between hadronic and e/m subcascades in favor of the former [cf.Eq. ( 8)], both effects contributing to an elongation of the nuclear cascade profile.The resulting X µ max values exceed the ones of the original QGSJET-III by up to ∼ 10 g/cm 2 , see Fig. 20 (right).

Conclusions
In this work, we investigated the possibility to increase the predicted EAS muon content, within the framework of the QGSJET-III model, restricting ourselves to the standard physics picture.We started with specifying the kinematic range for secondary hadron production, which is of relevance for model predictions regarding the muon number N µ and the maximal muon production depth X µ max , and demonstrated that these quantities are governed by forward pion production in hadronair collisions: with the corresponding pion energy fractions, taken from their parent hadrons, ranging between few percent and few tens of percent.Further, we investigated the impact on the predicted N µ and X µ max of various modifications of the treatment of hadronic interactions, in particular the model calibration to accelerator data, the amount of the "glue" in the pion, and the energy dependence of the pion exchange process.The strongest N µ enhancement reaching ∼ 10% level has been obtained when we increased the yields of secondary kaons and (anti)nucleons in the model, adjusting its predictions to NA61 data on pion-carbon interactions at fixed target energies.However, such modifications should be considered as extreme ones, since they lead to a serious tension with results of other experiments, regarding the yields of such hadrons both in pion-proton and, especially, in proton-proton collisions.
As for the gluon content of the pion, even considering extreme changes of the momentum partition between valence quarks and gluons, we obtained rather insignificant modifications of the predicted N µ and X µ max .This is because such changes impact pion-air collisions at the highest energies only, in the atmospheric nuclear cascade, while being of minor importance for the bulk of such interactions at lower energies, at later stages of the cascade development.
Regarding the pion exchange process in pion-air collisions, the (slow) decrease with energy of the corresponding probability can only be tamed if one assumes that absorptive corrections to the process are rather weak.However, in such a case, a significant part of the cross section for the virtual pion interaction with the target nucleus corresponds to elastic scattering.The corresponding very scarce hadron production results in a decrease of the predicted N µ .
Regarding X µ max , for some of the considered modifications of particle production, we obtained up to ∼ 10 g/cm 2 changes of the maximal muon production depth, typically yielding larger values of X µ max .This aggravates the tension with observations of the Pierre Auger Observatory.Overall, we see no possibility to increase the predicted EAS muon content by more than 10%, without entering in a serious contradiction with relevant accelerator data, while staying within the standard physics picture.

Figure 6 :
Figure 6: Momentum distributions (in lab.frame) of charged pions produced in π − C collisions at 350 GeV/c, obtained using the QGSJET-II-04 model and applying the "splitting" procedure to half the secondary pions, compared to NA61 data [27] (points): π + -solid line and filled squares, π − -dashed line and open squares.
, again for p-induced EAS of E 0 = the x cut dependence for the change of the calculated EAS

Figure 11 :
Figure 11: Gluon PDF of the pion at Q 2 0 = 2 GeV 2 in the default QGSJET-III model (solid line) and for the enhanced gluon content of the pion, as discussed in the text (dashed line).

Figure 13 :
Figure13: Momentum distributions (in lab.frame) of K + (top panels) and K − (bottom panels) produced in π − C collisions at 158 GeV/c (left) and 350 GeV/c (right), as calculated using the default QGSJET-III model (solid lines) or considering an enhancement of kaon production, as discussed in the text (dashed lines), compared to NA61 data[27] (points).

Figure 15 :
Figure15: Feynman x distributions (in c.m. frame) of ρ 0 mesons produced in π − C collisions at 158 GeV/c (left) and 350 GeV/c (right), as calculated using the default QGSJET-III model (solid lines) or considering an enhancement of ρ meson production, as discussed in the text (dashed lines), compared to NA61 data[32] (points).

Figure 16 :
Figure16: Feynman x distributions (in c.m. frame) of protons (left) and antiprotons (right) in π − p collisions at 360 GeV/c, as calculated using the default QGSJET-III model (solid lines) or considering an enhancement of (anti)nucleon production, as discussed in the text (dashed lines), compared to LEBC-EHS data[33] (points).

Figure 17 :
Figure17: Energy dependence of the relative change of the muon number N µ at sea level (left) and of the modification of the maximal muon production depth X µ max (right), both for E µ > 1 GeV, for proton-initiated EAS, with respect to the corresponding predictions of the default QGSJET-III model, for the considered modifications of the model: enhancement of (anti)nucleon, kaon, and ρ meson production -solid, dashed, and dash-dotted lines, respectively.

Figure 18 :
Figure18: Feynman x distribution (in c.m. frame) of ρ 0 mesons produced in π − C collisions at 158 GeV/c, calculated using the default QGSJET-III model (solid line) or applying the alternative treatment for the pion exchange process, as discussed in the text (dashed line).The corresponding partial contributions of the pion exchange process are shown by dotted and dash-dotted lines, respectively.

Figure 20 :
Figure20: Energy dependence of the relative change of the muon number N µ at sea level (left) and of the modification of the maximal muon production depth X µ max (right), both for E µ > 1 GeV, for proton-initiated EAS, with respect to the corresponding predictions of the default QGSJET-III model, for the alternative treatment of the pion exchange process, discussed in the text. π