QGSJET-III model of high energy hadronic interactions: II. Particle production and extensive air shower characteristics

The hadronization procedure of the QGSJET-III Monte Carlo (MC) generator of high energy hadronic interactions is discussed. Selected results of the model, regarding production spectra of secondary particles, are presented in comparison to experimental data and to the corresponding predictions of the QGSJET-II-04 MC generator. The model is applied to calculations of basic characteristics of extensive air showers initiated by cosmic ray interactions in the atmosphere and the results are compared to predictions of other MC generators of cosmic ray interactions.


Introduction
All current Monte Carlo (MC) generators of high energy hadronic interactions rely on the qualitative picture of quantum chromodynamics (QCD) in the sense that such interactions are assumed to be mediated by parton cascades developing between the colliding hadrons (nuclei).Particle production occurs when the coherence of (some of) these parton cascades is broken by the scattering process and final s-channel partons convert into hadrons.Since such a hadronization is essentially a nonperturbative process, the corresponding treatment relies on phenomenological approaches.The most popular ones are string fragmentation procedures implying a formation of strings of color field between final partons.With such partons flying apart, the tension of a string rises until it breaks up, with the color field being neutralized via a creation of quarkantiquark and diquark-antidiquark pairs from the vacuum, ending up with a formation of col-orless parton clusters identified with hadrons [1].
Needless to say, modeling such a process involves a substantial number of adjustable parameters and requires a thorough calibration with various experimental data.Nevertheless, this does not necessarily spoil the predictive power of hadronic interaction models since the mechanism is assumed to be universal.E.g., moving to higher energies or replacing the interacting hadrons, say, protons, by other hadrons or by nuclei, one arrives to a different parton configuration.However, the hadronization procedure remains unaffected by such changes, hence, the conversion of partons into hadrons can be described using one and the same approach, regardless the type and kinematics of the interaction.This allows one to test both the validity of an interaction model itself and the quality of the hadronization mechanism employed and, in case both perform satisfactorily enough, to use the complete model for predictive calculations, e.g., regarding the properties of socalled extensive air showers (EAS) initiated by interactions of high energy cosmic rays in the atmosphere (see [2] for a review).
In the particular case of the QGSJET-III MC generator [3], the treatment of the hadronization follows closely the one of the original QGSJET model [4,5], with a number of technical improvements, as discussed in Section 2. In Section 3, a comparison with selected accelerator data will be presented, followed in Section 4 by a discussion of basic properties of proton-nitrogen interactions, relevant for EAS predictions.In turn, Section 5 is devoted to calculations of extensive air shower characteristics.The conclusions will be given in Section 6.

Hadronization procedure of the QGSJET-III model
As described in [3], the procedure for simulating a particular hadronic (nuclear) collision event in the QGSJET-III model consists of two parts.First, a "macro-configuration" for an inelastic interaction is defined, i.e., the complete "net" of "elementary" parton cascades giving rise to secondary hadron production and being described by cut Pomerons is reconstructed, based on the corresponding partial cross sections.The t-and s-channel evolution of perturbative parton cascades corresponding to parton virtualities |q2 | > Q 2 0 , Q 2 0 being the "infrared" cutoff of the model, is treated explicitly. 1 At the second step, one considers a formation of two strings of color field per cut Pomeron, stretched between constituent partons of the interacting hadrons (nuclei) or/and between the final schannel partons resulting from the treatment of perturbative parton cascades.The breakup and hadronization of those strings are modeled by means of a string fragmentation procedure.
While the 4-momenta of final s-channel partons resulting from perturbative parton cascades are already defined at the first step, this is not the case for constituent partons.For a given projectile or target hadron h, to which 2n strings are attached, the sharing of the, respectively, light cone plus (LC + ) or light cone minus (LC − ) momentum fractions x ± between constituent partons (string "ends") is performed according to the distribution Here the exponent α sea is used for sea (anti)quarks, considering only light u and d (anti)quarks as string ends, while the small x limit of a valence (anti)quark distribution follows the Regge behavior (∝ x −α R 2n−1 ), with α R = α ρ (0) = 0.5.Similarly, the LC momentum distribution of the remaining valence parton [(anti)quark for h being a meson or (anti)diquark for a baryon] follows the corresponding Regge behavior (∝ x −αr 2n ) [7,8,9]: the corresponding Regge intercepts being compiled in Table 1.Regarding the LC momentum distribution of constituent sea quarks, in principle, it should also be characterized by the Regge behavior (∝ x −α R i ) in the low x limit [7], which was adopted in the QGSJET [4,5] and QGSJET-II [10,11] models.However, here we treat the respective exponent as an adjustable parameter, 0.5 ≤ α sea < 1.The reason is that we consider Pomeron-Pomeron interactions, introducing a minimal Pomeron rapidity "size" ξ [3], for the Pomeron asymptotics to be applicable.Therefore, the deviation of α sea from α R is assumed to take effectively into account a formation of small rapidity gaps y gap < ξ, not treated explicitly by the model.
In turn, transverse momenta of string ends are sampled according to the distribution: the corresponding parameters being listed in Table 2.
The fragmentation of a string into secondary hadrons is performed considering iteratively an emission of a hadron from either string end (see, e.g.[12]), with the hadron LC momentum fraction x (LC + or LC − for a projectile or target side end, respectively) and the transverse momentum p t in the center-of-mass (c.m.) frame of  the string being sampled according to the distribution: using Eqs.(2)(3)(4)(5) for the relations between the parameters α q and the intercepts of the corresponding Regge trajectories.Here the large x limit is defined by the probability to slow down (move away in rapidity space) the antiparton [(anti)quark or (anti)diquark] q′ of the vacuumcreated pair, while the small x limit is obtained by slowing down both the original parton q and the parton q ′ attached to it to form the final hadron (q + q ′ → h) [9] (cf.Fig. 1).Vari- x and small x limits of the proton LC momentum distribution correspond to a large rapidity y separation between the u and ū of the vacuum-created uū pair.Left: the large x limit (1 − x ∝ e −y ) is obtained by slowing down the ū antiquark.Right: the small x limit (x ∝ e −y ) corresponds to slowing down both the u quark and the original ud diquark.
ous parton-antiparton pairs q′ q ′ created from the vacuum are sampled according to the corresponding probabilities a q ′ , relative to ūu and dd pairs (a u = a d = 1/2).The iterative string fragmentation is continued until the mass of the string reminder falls below a specified threshold 2 M thr , followed by a modeling of a two particle decay of the remaining string.The string fragmentation parameters are compiled in Table 3.
The above-discussed treatment takes into account the production of final "stable" hadrons: pions, kaons, nucleons, lambda and sigma baryons, including the corresponding antiparticles, where relevant.Contributions of decays of short-lived resonances are assumed here to be accounted for in the stable hadron yields via the duality principle.The exceptions are the formation of ∆ ++ and ∆ − resonances (and the corresponding antiparticles) by a fragmentation of constituent uu and dd (anti)diquarks as well as the production of ρ and η mesons, which is also treated explicitly.Regarding the latter, the probability to produce an η meson, relative to π 0 , is chosen as w η = 1/9 from the quark counting rules.Similarly, we use w ρ = 1/3 for the corresponding probability of ρ mesons, relative to pions, assuming that ∼ 50% of final pions originate from ρ decays.For brevity, we restrain here from discussing other technicalities of the hadronization procedure.
In addition to the basic procedure described above, a treatment of the pion exchange contribution to forward hadron production, i.e., the Reggeon-Reggeon-Pomeron (RRP) process, with R = π, has been considered in the QGSJET-III model3 [14].As discussed in more detail in [14], due to the low intercept of the pion Regge trajectory, α π (0) ≃ 0, the process does in the plots), for proton production in pp (left) and pC (right) collisions at 158 GeV/c, as calculated using the QGSJET-III (solid lines) and QGSJET-II-04 (dotted lines) models, compared to NA49 data [16,17] (points).not lead to a formation of a large rapidity gap between the "leading", i.e., most energetic, secondary hadron and the remaining hadronic system, with the consequence that absorptive corrections to the corresponding "bare" triple-Reggeon contribution are essentially the same as for the single cut Pomeron process.This motivated us to sample the RRP process, for R = π, as a fixed fraction, with probability w π h , of the single cut Pomeron contribution, for a given hadron h.The hadronization procedure remains practically unaltered in such a case, except that the leading hadron is sampled according to the respective x-and p t -distribution for the π-exchange process [14], while the pair of strings is attached here to the quark and antiquark of the virtual pion, instead of constituent partons of the initial hadron h.In other words, after the production of the leading hadron, the rest of the final state corresponds to an inelastic interaction 4 of the virtual pion with the target (in case h is a projectile hadron), with the respectively reduced c.m. energy for the collision.Added in a similar way is the f -Reggeon exchange process (RRP contribution, with R = f ), with α f (0) = 0.7 [15] and with the corresponding probabilities w f h .

Results for secondary particle production
Regarding applications of MC event generators to EAS simulations, of primary importance is an accurate enough description of forward hadron production.Here comes a significance of experimental measurements performed at fixed target energies, where forward spectra of secondary particles have been studied in great detail.We start in Fig. 2 with invariant cross sections in c.m. frame for proton production in pp and pC collisions at 158 GeV/c, comparing the corresponding results of the QGSJET-III model to the ones of QGSJET-II-04 and to measurements by the NA49 experiment.The same comparison regarding Feynman x (x F ) distributions of charged pions, neutrons, and antiprotons, also for pp and pC interactions at 158 GeV/c, is presented in Figs. 3, 4, and 5; the partial contribution to neutron production due to the pion exchange mechanism in the QGSJET-III model is shown separately in Fig. 4. In turn, x F distributions of secondary charged kaons, for pp collisions at 158 GeV/c, calculated using the QGSJET-III and QGSJET-II-04 models, are compared to NA49 data in Fig. 6.Regarding the production of p, n, π ± , and K ± , we observe relatively small differences between QGSJET-III and QGSJET-II-04, both models being in a satisfactory agreement with the mea- GeV/c, as calculated using the QGSJET-III (solid lines) and QGSJET-II-04 (dotted lines) models, compared to NA49 data [16,17] (points).Partial contributions of the pion exchange process in the QGSJET-III model are shown by dash-dotted lines.surements.This is, however, not the case for secondary antiprotons: the predicted x F distributions are substantially harder than observed by the experiment, especially, in the case of the QGSJET-III model.One possible way to improve the agreement with the data is to allow for a separate treatment of the hadronization of the "remnant" of the incoming proton, which, however, requires to introduce additional adjustable parameters to a model [21].Generally, the problem is related to an extrapolation of the high energy treatment towards low energies [22].Correspondingly, its ultimate solution requires a full scale treatment of the transition into the low energy regime, taking into consideration Reggeon exchanges and the Reggeon-Reggeon-Reggeon (RRR) contributions.Further, in Fig. 7, we compare x F spectra of positively and negatively charged secondary hadrons, for π + p and K + p collisions at 250 GeV/c, calculated using the QGSJET-III and QGSJET-II-04 models, to measurements by the NA22 experiment [23].Additionally, in Fig. 8, we compare the results of the two models for momentum distributions of various secondary hadrons produced in π − C collisions at 158 and 350 GeV/c to recent measurements of the NA61 experiment [24].For both models, we observe an acceptable agreement with the data both for the charged particle spectra in Fig. 7 and for the pion distributions in Fig. 8.However, charged kaon production in π − C collisions appears to be underestimated by QGSJET-III by ∼ 30%, compared to the NA61 data.This is somewhat surprising since no such a "kaon deficit" was observed for pp collisions (cf.Fig. 6).Even more striking is the considerable underestimation by the QGSJET-III and QGSJET-II-04 models of proton and antiproton production in π − C collisions, compared to the NA61 measurements.While these discrepancies may, in principle, be attributed to potential deficiencies of the hadronization procedures of the models, it is noteworthy that no such a striking underestimation of secondary proton and antiproton yields is observed if we compare the predictions of QGSJET-III for p and p production in π − p collisions at 360 GeV/c to the data of the LEBC-EHS experiment [25], as one can see in Fig. 9.
As discussed in [13], of particular importance for model predictions regarding EAS muon content is a satisfactory description of forward production of ρ-mesons in pion-air collisions, which is dominated by the pion exchange process.In Fig. 10, we compare x F distributions of ρ 0 mesons produced in π − C collisions at 158 and 350 GeV/c, calculated using the QGSJET-III model, to the respective data of the NA61 experiment, while Fig. 11 contains a comparison of the model predictions for the spectra of ρ 0 and ρ + mesons produced in π + p interactions at 250 GeV/c to measurements of the NA22 experiment.Both in Fig. 10 and in Fig. 11, we show also the partial contribution of the pion exchange mechanism to ρ meson production, which clearly dominates the forward parts of the spectra.Overall, the agreement is quite good, notably, when comparing to the recent measurements of NA61.
Coming now to measurements at the Large Hadron Collider (LHC), those mostly provided information on hadron production in pp, pA, and AA collisions, at central (y ∼ 0) rapidities in c.m. frame.A comparison of the pseudorapidity density of charged hadrons, dN ch pp /dη, in pp collisions at √ s = 0.9, 2.36, 7, and 13 TeV, calculated using the QGSJET-III and QGSJET-II-04 models, to the respective ATLAS data 5 [28,29] is presented in Fig. 12.Since these experimental data have been used to tune the QGSJET-III model, the overall satisfactory agreement comes at no surprise.Much less trivial is the agreement of the predictions of both models to the results of a combined measurement 6 by the CMS and TOTEM experiments of dN ch pp /dη over a wide η range, in pp collisions at √ s = 8 TeV [30], see Fig. 13.As discussed in [31], this set of experimental data allows one to test very basic model assumptions regarding the momentum distributions of constituent partons.
We also show in Fig. 14 the calculated pseudorapidity density of charged hadrons, dN ch pPb /dη, for proton-lead collisions at √ s = 5 and 8 TeV, compared to the data of the CMS experiment [32].Here the results of the QGSJET-III model agree with the measurements substantially better, compared to the ones of QGSJET-II-04.
As an example comparison with identified secondary hadron spectra, we plot in Fig. 15 transverse momentum distributions of charged pions, charged kaons, and protons plus antiprotons, produced in pp collisions at √ s = 7 TeV, cal- 5 The experimental event selection requires at least one charged secondary hadron of pt > 0.5 GeV at |η| < 2.5. 6The experimental event selection requires at least one charged secondary hadron at 5.3 < |η| < 6.5.culated with the QGSJET-III and QGSJET-II-04 models, being confronted to the corresponding data of the ALICE experiment.The somewhat deficient description of these distributions by the two models may be related to the simplified string fragmentation procedure employed, based on the duality principle: not treating explicitly the production and decays of short-lived hadronic resonances.
However, regarding the relevance to EAS calculations, of primary significance is the energy-dependence of forward hadron production.In that respect, of particular importance are diffraction studies at LHC.In Table 4, we compare the predictions of the QGSJET-III and QGSJET-II-04 models, regarding partial cross sections for single diffraction-like (SD-like) events at √ s = 7 TeV, for different intervals of the diffractive state mass M X , to the respective data of the TOTEM experiment [34].While the cumulative cross section for the full range 3.1 < M X < 3100 GeV, studied by TOTEM, is well in agreement between both models and the data, the predicted rates of SD-like events for intermediate diffractive masses, 7.7 < M X < 1150 GeV, are somewhat lower than observed by TOTEM.On the other hand, the production of large mass states, 1150 < M X < 3100 GeV, is overestimated in both models by more than a factor of two, compared to the TOTEM mea- M X range, GeV  4: Predictions of the QGSJET-II-04 and QGSJET-III models for cross sections of SD-like events (in mb) at √ s = 7 TeV, for different ranges of mass M X of diffractive states produced, compared to TOTEM data [34].Values in brackets correspond to the respective contributions of nondiffractive events.surements.Taking into account that the event rate in this mass interval is dominated by a formation of random rapidity gaps in nondiffractive hadron production (see the corresponding values in brackets in Table 4), primarily, by RRP processes (see, e.g.[35] for the corresponding discussion), this may indicate a certain overestimation of the RRP contributions or/and some deficiencies of the hadronization procedures of both models.
Regarding the so-called inelasticity, i.e., the relative energy loss of most energetic secondary hadrons, important constraints come from measurements of forward neutron production by the LHCf experiment [36,37].When interpreting those experimental results, it is important to take into consideration the contribution of the t-channel pion exchange [38,39,40].Here the data of the NA49 experiment on neutron production in pp and pC collisions at 158 GeV/c [16,17] (cf.Fig. 4) can be combined with those of LHCf at LHC energies in order to test the implementation of the pion exchange mechanism in the QGSJET-III model [14], notably, regarding the predicted energy-dependence over a wide c.m. energy range, 17 < √ s < 13000 GeV.Comparing the calculated neutron energy spectra to the observations in Figs.16-17, we see a satisfactory agreement of the model with the measurement at √ s = 13 TeV, while there is a certain underestimation of the neutron yield at √ s = 7  frame, produced in pp collisions at different √ s (from top to bottom: 13, 7, 2.36, and 0.9 TeV), as calculated using the QGSJET-III (solid lines) and QGSJET-II-04 (dotted lines) models, compared to ATLAS data [28,29] (points).
TeV, compared to LHCf data. 7The latter is ac- 7 As one can see in Figs. 4 and 16, the agreement with the data on neutron production, both from the NA49 experiment at 158 GeV/c and from LHCf at √ s = 13 TeV, can be improved considering a smaller contribution of the pion exchange process, i.e., choosing a smaller value for the w π p parameter.This would, however, aggravate the tension with the LHCf data at √ s = 7 TeV. Figure 13: Pseudorapidity distribution of charged hadrons in c.m. frame, produced in pp collisions at √ s = 8 TeV, as calculated using the QGSJET-III (solid line) and QGSJET-II-04 (dotted line) models, compared to a simultaneous measurement by CMS and TOTEM [30] (points).tually surprising, given the small difference between the two c.m. energies for LHCf studies.
Another important benchmark from the LHCf experiment concerns the measured forward spectra of neutral pions [41], which is also of importance for astrophysical studies based on γrays (see, e.g.[42,43] for the corresponding discussion).The results of the QGSJET-III and QGSJET-II-04 models are confronted to the LHCf data in Fig. 18.While both models demonstrate a satisfactory agreement with the measurements, somewhat disturbing is a certain underestimation of π 0 production at p z ∼ 1 TeV by the QGSJET-III model, which may be related to its softer momentum distribution of constituent partons [larger α sea , cf.Eq. ( 1)].

Basic characteristics of proton-nitrogen collisions
Before proceeding with the model application to extensive air shower modeling, it is worth considering its predictions for some basic characteristics of hadron-nucleus collisions, relevant for EAS predictions.The first quantity to consider is the inelastic proton-air cross section σ inel p−air which defines the mean free path (m.f.p.) λ p = m air /σ inel p−air (m air being the average mass of air nuclei) of a primary CR proton in the atmosphere and impacts thereby many EAS characteristics, notably, the EAS maximum depth X max corresponding to the position of the maximum of charged particle profile of an air shower.The energy-dependence of σ inel p−air , calculated using the QGSJET-III, QGSJET-II-04, EPOS-LHC [44], and SIBYLL-2.3[21] models, is plotted in Fig. 19.Not surprisingly, there is a good agreement between the predictions of QGSJET-III and QGSJET-II-04, as well as with the ones of the other two models considered, since more or less the same set of accelerator data on the total and elastic protonproton cross sections had been used for the relevant parameter tuning in all the cases.The small differences between the results for σ inel p−air , all obtained using the Glauber-Gribov formalism [45,46], may stem from a choice of the nuclear ground state wave function or/and from a model treatment of the inelastic screening effects.With the most steep energy rise of σ inel p−air being predicted by EPOS-LHC, its cross section is ≃ 5% higher at E 0 = 10 10 GeV than the one of QGSJET-III.Yet the potential impact of that difference on X max at the highest energies is < 3 g/cm 2 .
Next, we consider the rate of the inelastic   diffraction, which may also impact X max predictions noticeably because of its typically small inelasticity K inel , i.e., the relative energy loss of leading (most energetic) secondary nucleons.This is particularly so for target diffraction characterized at very high energies by a tiny inelasticity K inel ≃ M 2 X /s, M X being the diffractive state mass.The effect of such "quasi-elastic" collisions is equivalent to a reduction of the inelas-tic proton-air cross section, giving rise to an increase of m.f.p of primary CR protons: In Fig. 20, we compare the energy-dependence of the rate of diffractive-like interactions, for proton-nitrogen collisions, i.e., the fraction of events characterized by K inel < 0.1, predicted by different models.It is worth reminding that, apart from diffractive interactions, a sizable contribution to this rate comes from a formation of random rapidity gaps in nondiffractive collisions, notably, by RRP processes.While the results of QGSJET-III and QGSJET-II-04 agree with each other within 10%, the overall spread of the model predictions reaches here 30%, stemming both from different approaches to the treatment of diffraction and from differences regarding the hadronization procedure, the latter potentially having a strong impact on the probability to create random rapidity gaps.Yet, since such a variation of the rate of diffractive-like interactions can produce only ≃ 3% change of λ p , the corresponding impact on X max is limited by ≃ 2 g/cm 2 [cf.Eq. ( 8)].
Further we consider the inelasticity K inel of general proton-nitrogen collisions, which impacts strongly the energy loss of leading secondary nucleons in the cause of the nuclear cascade development, influencing significantly the predicted X max (see, e.g.[47]).The energy dependence of K inel , predicted by different models, is plotted in Fig. 21.Despite large differences between QGSJET-II-04 and QGSJET-III, both regarding the new developments in the latter, notably, the treatment of higher twist corrections to hard parton processes [3,48] and the implementation of the pion exchange process [14], and concerning the values of important model parameters (e.g., the parameter α sea governing the momentum distribution of constituent partons), their predictions for K inel are remarkably identical, both for the energy dependence and for the absolute values.A qualitatively similar behavior of K inel is predicted by EPOS-LHC.In contrast, in case of the SYBYLL-2.3model, the inelasticity depends rather weakly on the collision energy, which is related to very basic model assumptions regarding the structure of constituent parton Fock states and impacts strongly the model results for X max [31].
Regarding the EAS muon content, the simple qualitative model of Heitler [49] implies a correlation of this observable with the number of "stable" secondary hadrons [(anti)nucleons, kaons, and charged pions], i.e., those which have significant chances to interact in the atmosphere, before decaying.More precisely, since the cascade nature of extensive air showers enhances the importance of forward hadron production, the relevant quantity is the second moment of the distribution dn stable /dx E of the energy fraction x E = E/E 0 of such hadrons [50] (see also [51,52] for recent studies), i.e., the total fraction of the parent hadron energy taken by all stable hadrons, ) In Fig. 22, we compare the energy-dependence of this quantity for pion-nitrogen collisions, Figure 22: Energy-dependence of x E n stable , for π + 14 N collisions, calculated using different models.The notations for the lines are the same as in Fig. 19.calculated using the QGSJET-III, QGSJET-II-04, EPOS-LHC, and SIBYLL-2.3models.While there are only minor differences between QGSJET-III and QGSJET-II-04 for the predicted x E n stable , the results of the other two models are quite different, which is related to a copious forward production of, respectively, (anti)baryons and ρ-mesons in the EPOS-LHC and SIBYLL-2.3models. 8verall, the results of the QGSJET-III and QGSJET-II-04 models, plotted in Figs.19-22, appear to be rather similar to each other, suggesting that EAS predictions of the two models should not differ significantly.

Predictions for EAS characteristics
In Fig. 23, predictions of the QGSJET-III, QGSJET-II-04, EPOS-LHC, and SIBYLL-2.3models for the energy-dependence of the maximum depth X max and of its fluctuation, σ Xmax , of proton-induced extensive air showers are compared.As anticipated in Section 4, the results of QGSJET-III and QGSJET-II-04 are rather similar to each other, their X max values differing by less than 10 g/cm 2 .Such an agreement is very impressive, given the substantial differences with the predictions of the other two models.In turn, there is a very good agreement between all the models, regarding the calculated values of σ Xmax .This is not surprising since the quantity is mainly defined by the inelastic proton-air cross section (see, e.g.[54]), while the impact of uncertainties regarding the rate of inelastic diffraction does not exceed few g/cm 2 [55].Further, the calculated energy-dependence of the muon number N µ (at sea level) of protoninitiated EAS, for muon energies E µ > 1 GeV, is compared between the models in Fig. 24.Here again, the results of QGSJET-III and QGSJET-II-04 agree with each other within 5%.On the other hand, higher N µ predicted by EPOS-LHC and SIBYLL-2.3comes at no surprise, given their larger values for x E n stable (cf.Fig. 22 and the corresponding discussion in Section 4), being due to an enhanced forward production of, respectively, (anti)baryons and ρ-mesons in those models.
Finally, plotted in Fig. 25 is the energydependence of the maximal muon production depth X µ max , i.e., the depth in the atmosphere corresponding to a maximal number of muons being produced by hadron decays, calculated using the different models.As discussed in [56], this quantity is very sensitive to various aspects of pion-air collisions and, thus, can be used to test and constrain the respective model treatment, based on CR data.Even here, despite the fact that the energy dependence of hadron production in pion-nucleus collisions is weakly constrained by available accelerator data, the agreement between the results of QGSJET-III and QGSJET-II-04 is surprisingly good.On the other hand, the substantially larger values of X µ max , predicted by the EPOS-LHC and SIBYLL-2.3models, appear to be in a strong contradiction to the corresponding measurements by the Pierre Auger Observatory [57].

Conclusions
In this work, we discussed in some detail the hadronization procedure of the QGSJET-III MC generator and presented selected results of that generator, regarding secondary hadron production, in comparison to experimental data and to the corresponding results of the previous model version, QGSJET-II-04.Overall, the quality of agreement with the data for QGSJET-III remained at approximately the same level, as in the case of QGSJET-II-04, certain improvements being mostly related to the treatment of the pion exchange process in hadron-proton and hadron-nucleus collisions.
Comparing the predictions of different MC generators, regarding basic characteristics of proton-induced extensive air showers, we observed an outstanding agreement between the results of QGSJET-III and QGSJET-II-04, in spite of an implementation of new theoretical approaches and a number of technical improvements in the new model.This may indicate that the treatment of relevant aspects of hadronic interaction physics is already sufficiently constrained by available accelerator data.On the other hand, in view of rather different results of the other MC generators, regarding, e.g., the air shower maximum depth or EAS muon content, one may rise a question whether the similarity of extensive air shower predictions of QGSJET-III and QGSJET-II-04 is rather a consequence of the underlying theoretical approaches shared by the two models or/and common deficiencies regarding the treatment of certain aspects of the interaction physics.Therefore, a general analysis of potential model uncertainties for EAS predictions is a warranted task.Such an analysis will be presented elsewhere [52,58].

Figure 1 :
Figure1: Schematic view of the fragmentation of a fast ud diquark into a proton: both large x and small x limits of the proton LC momentum distribution correspond to a large rapidity y separation between the u and ū of the vacuum-created uū pair.Left: the large x limit (1 − x ∝ e −y ) is obtained by slowing down the ū antiquark.Right: the small x limit (x ∝ e −y ) corresponds to slowing down both the u quark and the original ud diquark.

Figure 3 :Figure 4 :
Figure 3: x F -distributions of charged pions (in c.m. frame) in pp (left) and pC (right) collisions at 158 GeV/c, as calculated using the QGSJET-III model, compared to NA49 data [18, 19] (points): π + -solid lines and filled squares, π − -dashed lines and open squares.The corresponding results of the QGSJET-II-04 model are shown by dotted and dash-dotted lines for π + and π − , respectively.

Figure 6 :
Figure 6: x F -distributions of charged kaons (in c.m. frame) in pp collisions at 158 GeV/c, as calculated using the QGSJET-III model, compared to NA49 data [20] (points): K + -solid line and filled squares, K − -dashed line and open squares.The corresponding results of the QGSJET-II-04 model are shown by dotted and dash-dotted lines for K + and K − , respectively.

Figure 7 :
Figure 7: x F -dependence of production cross sections in c.m. frame for charged hadrons (h ± ) in π + p (left) and K + p (right) collisions at 250 GeV/c, as calculated using the QGSJET-III model, compared to NA22 data [23] (points): h + -solid lines and filled squares, h − -dashed lines and open squares.The corresponding results of the QGSJET-II-04 model are shown by dotted and dash-dotted lines for h + and h − , respectively.

Figure 8 :
Figure 8: Momentum distributions in laboratory frame of charged pions (top panels), charged kaons (middle panels), protons and antiprotons (bottom panels), produced in π − C collisions at 158 GeV/c (left) and 350 GeV/c (right), as calculated using the QGSJET-III model, compared to NA61 data [24] (points).Solid lines and filled squares correspond to positively charged hadrons while dashed lines and open squares refer to negatively charged hadrons.The corresponding results of the QGSJET-II-04 model are shown by dotted and dash-dotted lines for positively and negatively charged hadrons, respectively.

Figure 9 :
Figure 9: x F -dependence of p t -integrated invariant cross sections in c.m. frame, x E dσ/dx F , for proton and antiproton production in π − p collisions at 360 GeV/c, as calculated using the QGSJET-III model, compared to LEBC-EHS data [25] (points): p -solid line and filled squares, p -dashed line and open squares.The corresponding results of the QGSJET-II-04 model are shown by dotted and dash-dotted lines for p and p, respectively.

Figure 10 :
Figure 10: x F -distributions of ρ 0 -mesons (in c.m. frame) in π − C collisions at 158 GeV/c (solid line) and 350 GeV/c (dashed line), as calculated using the QGSJET-III model, compared to NA61 data [26] (points).Partial contribution of the pion exchange process for π − C collisions at 158 GeV/c is shown by a dash-dotted line.

Figure 11 :Figure 12 :
Figure11:x F -dependence of production cross sections in c.m. frame for ρ 0 (left) and ρ + (right) mesons in π + p collisions at 250 GeV/c, as calculated using the QGSJET-III model, compared to NA22 data[27] (points).Partial contributions of the pion exchange process are shown by dashdotted lines.

Figure 14 := 5
Figure 14: Pseudorapidity distributions of charged hadrons in c.m. frame, produced in pPb collisions at different √ s, as calculated using the QGSJET-III model, compared to CMS data [32] (points): √ s = 5 TeV -solid line and filled squares, √ s = 8 TeV -dashed line and open squares.The corresponding results of the QGSJET-II-04 model are shown by dotted and dash-dotted lines for √ s = 5 TeV and √ s = 8 TeV, respectively.

Figure 15 :
Figure 15: Transverse momentum distributions of identified hadrons at central rapidity (|y| < 0.5) in c.m. frame, for pp collisions at √ s = 7 TeV, calculated using the QGSJET-III model, in comparison to ALICE data [33] (points): π + +π − -solid line and filled circles, K + +K −dashed line and open squares, p+ p -dash-dotted line and filled squares.The corresponding results of the QGSJET-II-04 model are shown by dotted lines.

Figure 16 :Figure 17 : 4 <Figure 18 :
Figure 16: Neutron energy spectra (in c.m. frame) in pp collisions at √ s = 13 TeV, calculated using the QGSJET-III (solid lines) and QGSJET-II-04 (dotted lines) models, for different pseudorapidity intervals (as indicated in the plots), in comparison to LHCf data [37] (points).Partial contributions of the pion exchange process in the QGSJET-III model are shown by dash-dotted lines.

Figure 20 :
Figure20: Energy-dependence of the probability of diffractive-like interactions, for p 14 N collisions, calculated using different models.The notations for the lines are the same as in Fig.19.

Figure 21 :
Figure21: Energy-dependence of the inelasticity of p 14 N collisions, calculated using different models.The notations for the lines are the same as in Fig.19.

Figure 25 :
Figure25: Dependence on primary energy of the maximal muon production depth X µ max (E µ > 1 GeV) of proton-initiated EAS, calculated using different models.The meaning of the lines is the same as in Fig.23.

Table 2 :
Parameters of the energy-momentum sharing and hadronization procedures of the QGSJET-III model (excluding the parameters of the string fragmentation procedure).

Table 3 :
Parameters of the string fragmentation procedure of the QGSJET-III model.
Figure 2: x F -dependence of invariant cross sections in c.m. frame, x E dσ/dx F /d 2 p t (x E = 2E/