A quantitative analysis of model uncertainties for calculations of the maximum depth of protoninitiated extensive air showers EAS has been performed. Staying within the standard physics picture and using the conventional approach to the treatment of highenergy interactions we found that present uncertainties on the energy dependence of the inelastic cross section the rate of diffraction and the inelasticity of hadronic collisions allow one to increase the predicted average EAS maximum depth by about inlineformulammlmath displayinlinemmlmrowmmlmn10mmlmnmmlmtext mmlmtextmmlmtext mmlmtextmmlmi mathvariantnormalgmmlmimmlmommlmommlmsupmmlmrowmmlmicmmmlmimmlmrowmmlmn2mmlmnmmlmsupmmlmrowmmlmathinlineformula. Invoking more exotic assumptions regarding a potentially significant modification of the parton hadronization procedure by hypothetical collective effects we were able to change drastically the predicted energy dependence of the inelasticity of protonair interactions and to increase thereby the predicted EAS maximum depth by up to inlineformulammlmath displayinlinemmlmommlmommlmn30mmlmnmmlmtext mmlmtextmmlmtext mmlmtextmmlmi mathvariantnormalgmmlmimmlmommlmommlmsupmmlmrowmmlmicmmmlmimmlmrowmmlmrowmmlmn2mmlmnmmlmrowmmlmsupmmlmathinlineformula. However those latter modifications are disfavored by the data of the LHCf experiment regarding forward neutron production in protonproton collisions at the Large Hadron Collider and by measurements of the muon production depth by the Pierre Auger Observatory.
['Astrophysics and astroparticle physics', 'High Energy Physics - Phenomenology', 'Astrophysics - High Energy Astrophysical Phenomena']
Model uncertainties for the predicted maximum depth of extensive air showers
{'On the model uncertainties for the predicted maximum depth of extensive air showers': 'Sergey Ostapchenko and Günter Sigl  \nUniversität Hamburg, II Institut für Theoretische Physik, 22761 Hamburg, Germany  \nSeptember 10, 2024', 'Abstract': 'A quantitative analysis of model uncertainties for calculations of the maximum depth of proton-initiated extensive air showers (EAS) has been performed. Staying within the standard physics picture and using the conventional approach to the treatment of high energy interactions, we found that present uncertainties on the energy dependence of the inelastic cross section, the rate of diffraction, and the inelasticity of hadronic collisions allow one to increase the predicted average EAS maximum depth by about 10 g/cm 2 . Invoking more exotic assumptions regarding a potentially significant modification of the parton hadronization procedure by hypothetical collective effects, we were able to change drastically the predicted energy dependence of the inelasticity of proton-air interactions and to increase thereby the predicted EAS maximum depth by up to /similarequal 30 g/cm 2 . However, those latter modifications are disfavored by the data of the LHCf experiment, regarding forward neutron production in proton-proton collisions at the Large Hadron Collider, and by measurements of the muon production depth by the Pierre Auger Observatory.', '1 Introduction': 'One of the exciting directions of astroparticle physics research is connected to experimental studies of ultra-high energy cosmic rays (UHECRs). Because of the extremely low incoming flux of UHECRs, one is forced to employ indirect detection methods: inferring the properties of the primary particles from measured characteristics of extensive air showers (EAS) - nuclearelectromagnetic cascades initiated by interac-  \ntions of UHECRs in the atmosphere [1]. Consequently, the success of such studies depends on the quality of the description of EAS development by the corresponding numerical tools, like the CORSIKA program [2]. A special role in such EAS simulation procedures is played by Monte Carlo (MC) generators of high energy hadronic interactions, used to describe the backbone of an air shower - the cascade of interactions in the atmosphere of both primary cosmic ray (CR) particles and of secondary hadrons produced [3].  \nAll of the above fully applies to investigations of UHECR composition, based on studies of the longitudinal extensive air shower development, notably, on measurements of EAS maximum depth, X max , corresponding to the maximum of the charged particle profile of an air shower. Modern fluorescence detectors allow one to measure X max for individual events rather accurately (see [4] for a review). On the other hand, predictions of EAS simulation procedures for the average EAS maximum depth, 〈 X max 〉 , being driven by the treatment of primary particle interactions with air nuclei, are rather seriously constrained by experimental data from the Large Hadron Collider (LHC) (see, e.g. [5]). Therefore, somewhat surprising is a certain tension between the predictions of air shower simulations and the experimental data of the Pierre Auger Observatory, regarding distributions of the shower maximum depth [6, 7]. More specifically, Auger data demonstrate that both the elongation rate ER( E 0 ) = d 〈 X max ( E 0 ) 〉 /d lg E 0 and the width of X max distributions, σ ( X max ) , decrease with the primary energy E 0 faster than predicted by EAS simulations corresponding to an energy-independent UHECR composition, thereby indicating a change towards heav-  \nimaries at ultra-high energies [8, 9]. However, the observed energy dependence of σ ( X max ) implies a faster change of the composition, compared to the one deduced from ER( E 0 ) , e.g., if one employs the QGSJET-II-04 model [10, 11] for interpreting the measurements. As demonstrated in [7], a consistent interpretation of the data requires a significantly slower extensive air shower development, compared to current EAS simulation results. Consequently, of significant importance is to quantify the range of uncertainty for X max predictions, related to potential modifications of the treatment of high energy hadronic interactions, allowed by LHC data.  \nSuch a quantitative investigation is the subject of the current work. We consider all plausible changes of the interaction treatment, in the framework of the QGSJET-III model [12, 13], which may potentially improve the agreement of the model predictions with Auger data, while staying within the standard physics picture. Like in our previous study of model uncertainties for the predicted EAS muon content [14], our investigation is guided by three basic principles: 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.  \nThe outline of the paper is as follows. In Section 2, we identify basic characteristics of hadronic interactions, which are of direct relevance to X max predictions. In Section 3, we consider various modifications of the treatment of high energy interactions, which change such characteristics, and study the corresponding consequences both for the predicted EAS maximum depth and for comparisons with relevant LHC data. We summarize our investigation in Section 4.', '2 Characteristics of high energy interactions, relevant for X max predictions': "Unlike the EAS muon content depending on the whole history of the nuclear cascade in the  \natmosphere, the maximum depth of a protoninitiated extensive air shower is largely governed by an interaction of the primary particle. As demonstrated, e.g., in [14], modifying the treatment of secondary pion-air collisions has a rather weak impact on the calculated 〈 X max 〉 . Therefore, predictions of EAS simulations for the longitudinal extensive air shower development are rather seriously constrained by available LHC data on proton-proton and, to some extent, proton-nucleus interactions.  \nOf primary importance for 〈 X max 〉 predictions is the proton-air inelastic cross section σ inel p -air since it controls the mean free path of protons in air, λ p = m air /σ inel p -air , m air being the average mass of air nuclei, and thereby the starting point for a nuclear cascade. The energy rise of σ inel p -air is the main factor which causes the decrease of the elongation rate ER( E 0 ) with increasing primary energy E 0 . The inelastic proton-proton cross section σ inel pp has been measured at LHC with a high precision; the difference between the results of the TOTEM and ATLAS experiments, based on Roman Pot techniques, is /similarequal 2 . 6 % at the center-of-mass (c.m.) energy √ s = 13 TeV [15, 16]. Using the corresponding values for calculating σ inel p -air , within the Glauber-Gribov formalism [17, 18], one obtains even a smaller difference for the proton-air cross section since it is largely dominated by the nuclear size.  \nAlso of importance for calculations of EAS maximum depth is the treatment of inelastic diffraction. Diffractive interactions of primary protons are characterized by a small inelasticity K inel p -air , i.e., by a small relative energy loss of 'leading' (most energetic) secondary nucleons. This is especially so for a diffractive excitation of a target nucleus, in which case the incoming proton looses a tiny fraction of its initial energy, K inel p -air /similarequal M 2 X /s /lessmuch 1 , M X being the diffractive state mass, i.e., one essentially deals with a quasi-elastic collision. Therefore, the main effect of diffraction amounts to a 'renormalization' of the inelastic proton-air cross section, merely subtracting its diffractive part, σ diffr p -air : σ inel p -air → σ inel p -air -σ diffr p -air , and of proton mean free path:  \nλ p → λ p (1 + σ diffr p -air /σ inel p -air ) . (1)  \nIn addition, the inelastic diffraction is closely related to the inelastic screening effect [18]; a  \nhigher diffraction rate corresponds to a smaller σ inel p -air , for a given σ inel pp , and vice versa.  \nFinally, calculations of 〈 X max 〉 depend on the predicted energy dependence of the average inelasticity of proton-air interactions, which is governed by the treatment of nondiffractive collisions. A smaller K inel p -air corresponds to a slower energy dissipation for leading nucleons, hence, to a larger fraction of the energy of the primary particle, being retained in the hadronic cascade, and to a slower EAS development. An enhancement of secondary particle production, with increasing energy, due to the energy rise of the multiple scattering rate, inevitably leads to an increase of the inelasticity. However, the speed of the energy rise of K inel p -air is rather weakly constrained by LHC data and is highly modeldependent [19], being currently the main source of model uncertainties for 〈 X max 〉 predictions.  \nComing now to fluctuations of EAS maximum depth, characterized by the standard deviation σ ( X max ) , those are rather insensitive to the average inelasticity, being dominated by variations of the proton mean free path, hence, depending mostly on σ inel p -air . In addition, they are influenced somewhat by the rate of inelastic diffraction, where the effective 'renormalization' of the inelastic proton-air cross section and the inelastic screening effect work in the same direction, e.g., both effects contribute to enlarging σ ( X max ) in case of a higher diffraction. Overall, the uncertainties of predictions for σ ( X max ) are rather small, as noticed already in [20], notably, thanks to measurements of total, elastic, and diffractive pp cross sections by the TOTEM and ATLAS experiments at LHC [15, 16, 21, 22, 23, 24, 25, 26, 27, 28]. For example, the current /similarequal 3 % difference between the results of TOTEM and ATLAS for σ inel pp at √ s = 13 TeV translates itself into /similarequal 1 . 5 g/cm 2 variation of σ ( X max ) [cf. Eq. (1)]. On the other hand, the impact on the fluctuations of X max of present uncertainties for σ diffr pp is at the level of few g/cm 2 [29].  \nRegarding average characteristics of extensive air showers initiated by primary nuclei, e.g. 〈 X max 〉 , those are well described by the superposition model, i.e., coincide with a good accuracy with the corresponding contributions of A proton-induced EAS of A times smaller energy, for a primary nucleus of mass number A (e.g. [30, 31]). This follows from the relation between  \nthe inelastic nucleus-air cross section σ inel A -air and the mean number of interacting projectile nucleons 〈 ν A 〉 , per inelastic collision [32]:  \n〈 ν A 〉 = Aσ inel p -air σ inel A -air , (2)  \nand from the possibility to approximate forward production spectra of an A -air collision, for given 〈 ν A 〉 , by the ones of proton-air interactions:  \ndN X A -air ( E 0 , E ) dE ∣ ∣ ∣ ∣ ∣ 〈 ν A 〉 →〈 ν A 〉 dN X p -air ( E 0 /A, E ) dE , (3)  \nfor any secondary particle X . As a consequence of Eq. (2), mean number of projectile nucleons interacting in a given depth interval, for an air shower initiated by nucleus A , is the same as for a superposition of A proton-induced EAS of A times smaller energy [31].  \nYet the superposition model is invalid for fluctuations of EAS characteristics [30, 31, 33], e.g., for σ ( X max ) , which are dominated by variations of the impact parameter of nucleus-air collisions, causing large fluctuations of the number of interacting projectile nucleons. While those are well-defined in the Glauber-Gribov approach, an additional contribution to σ ( X max ) comes from a fragmentation of the spectator part of the projectile nucleus [30, 31]. For example, for an extensive air shower initiated by a primary iron nucleus, one obtains a factor two difference between the predicted values of σ ( X max ) , when considering two extreme assumptions: a full breakup of the nuclear spectator part into separate nucleons or keeping all non-interacting nucleons together, as a single secondary nucleus [31]. Since experimental data on nuclear fragmentation are available at fixed target energies only, this could have constituted a serious source of uncertainty for predictions of σ ( X max ) , for nucleus-induced EAS. However, the relative yields of various nuclear fragments are energy-independent, above few GeV per nucleon [34], which allows one to reliably calibrate the fragmentation procedures, based on fixed target data, and to safely extrapolate them to UHECR energies.  \nAs follows from the above discussion, the main model uncertainties for 〈 X max 〉 and σ ( X max ) stem from the treatment of proton-air inter-  \nns, while an extension to the case of nuclear primaries is rather well defined. Therefore we concentrate in the following on a study of proton-induced extensive air showers.", '3.1 Smaller cross section and higher diffraction rate': 'As discussed in Section 2, both a smaller protonair cross section and a higher diffraction rate lead to a larger 〈 X max 〉 predicted. Here we are going to study a combined effect of both, increasing the rate of low mass diffraction (LMD) in the QGSJET-III model. The LMD treatment in QGSJET-III is based on the Good-Walker (GW) formalism [35]: considering a hadron to be represented by a superposition of a number of GW Fock states characterized by different sizes and different (integrated) parton densities. Since such states undergo different absorption during a scattering process, one generally has transitions of the initial hadrons into various low mass excited states (see, e.g. [36], for the corresponding discussion). Obviously, to enhance the LMD rate, one thus needs to enlarge the difference between transverse sizes of the GW states. Using a twice smaller value for the ratio of the squared radii of the smallest and largest GW states of the proton (parameter d p in QGSJET-III [12]), we obtain cross sections for single diffractivelike (SD-like) events 1 , for different intervals of diffractive state mass M X , listed in Table 1. It is easy to see that the LMD rate, for M X < 3 . 4 GeV, is enhanced by /similarequal 30 %, compared to the default QGSJET-III model, while being still compatible with the observations of the TOTEM experiment [22]. Since the considered modification gives rise to a strong enhancement of the inelastic screening effect (see, e.g. [37], for the corresponding discussion), it leads also to a sizable  \nreduction of the calculated total, elastic, and inelastic pp cross sections plotted in Fig. 1, which now become compatible with the data of the ATLAS experiment [16, 24, 25].  \nFigure 1 : C.m. energy dependence of the total, inelastic, and elastic pp cross sections, calculated with the default QGSJET-III model (red solid lines) and with the option characterized by an enhanced diffraction (blue dashed lines), compared to experimental data (points) from [15, 16, 38].  \n<!-- image -->  \nFigure 2 : Laboratory (lab.) energy dependence of σ inel p -air , calculated with the default QGSJETIII model (red solid line) and with the option characterized by an enhanced diffraction (blue dashed line).  \n<!-- image -->  \nHowever, as anticipated in Section 2, the impact of such changes on the calculated inelastic  \nTable 1 : Cross sections (in mb) of SD-like pp collisions at √ s = 7 TeV, for different ranges of mass M X of diffractive states produced, calculated with the default QGSJET-III model and with the option characterized by an enhanced diffraction, compared to TOTEM data.  \nFigure 3 : Lab. energy dependence of the probability of diffractive-like interactions (left) and of the inelasticity (right), for p 14 N collisions, calculated with the default QGSJET-III model (red solid lines) and with the option characterized by an enhanced diffraction (blue dashed lines).  \n<!-- image -->  \n0  \n<!-- image -->  \n0  \nproton-air cross section plotted in Fig. 2 is much more moderate: reaching only /similarequal 1 % level at E 0 = 10 19 eV. On the other hand, one obtains a significant enhancement of diffraction in protonnucleus collisions. For example, for proton interactions with the most abundant air element, nitrogen, the increase of the rate of diffractive-like collisions characterized by a small inelasticity, K inel p N < 0 . 1 , reaches /similarequal 15 % level, as one can see in Fig. 3 (left). It is also noteworthy that such an enhancement of diffraction has a small impact on the average inelasticity of proton-nitrogen interactions [cf. Fig. 3 (right)], the latter being dominated by the treatment of nondifractive collisions.  \nAs is obvious from the above-presented results, the changes of the predicted 2 〈 X max 〉 and σ ( X max ) plotted in Fig. 4 are driven by the increased rate of diffractive-like proton-air collisions. The magnitude of these changes, up to /similarequal 8 g/cm 2 for 〈 X max 〉 and /lessorsimilar 4 g/cm 2 for  \nσ ( X max ) , agrees well with the corresponding results of [29].', '3.2 Slower energy rise of the inelasticity': "Thus, the only possibility to predict a substantially larger EAS maximum depth is to decrease the inelasticity of nondiffractive proton-air interactions. A violation of the Feynman scaling in hadronic collisions is a well established experimental fact, e.g., regarding a steep rise of central pseudorapidity η density of charged hadrons produced in pp interactions, dN ch pp /dη ∣ ∣ η =0 , with energy, which stems from an increase of multiple scattering rate. However, an approximate Feynman scaling of secondary hadron spectra at large values of Feynman x may still be allowed by experimental data.  \nGenerally, the rate of multiple scattering in hadronic collisions rises with energy, primarily, due to a fast increase of the rate of semihard processes leading to production of hadron  \nFigure 4 : Dependence on primary energy of the average maximum depth (left) and of the corresponding standard deviation (right), for proton-initiated EAS, calculated with the default QGSJET-III model (red solid lines) and with the option characterized by an enhanced diffraction (blue dashed lines).  \n<!-- image -->  \n(mini)jets, which is related, in turn, to a steep low x rise of parton momentum distribution functions of hadrons. Yet the rate of this increase may be tamed by nonlinear interaction effects, notably, regarding a copious production of minijets of relatively small transverse momenta p t . On the other hand, the inelasticity of high energy collisions depends strongly on the choice of momentum distributions of constituent partons for the interacting hadrons (nuclei), involved in numerous inelastic rescattering processes [19, 40, 41]. Choosing a softer distribution, ∝ x -α , for the fraction x of the initial light cone (LC) momentum of the parent hadron, taken by such a parton, i.e., using a larger value for α , one obtains a weaker impact of multiple scattering on forward hadron production, arriving to an approximate Feynman scaling at large x in the α → 1 limit [19].  \nTo investigate the impact of a reduced multiple scattering rate, we increase the strength of higher twist (HT) corrections to hard scattering processes in the QGSJET-III model, choosing a twice larger value, K HT = 5 , for the corresponding normalization parameter [12]. As one can see in Fig. 5, this way we decrease significantly, by up to /similarequal 30 %, the minijet production at small p t . It is noteworthy, however, that the considered change is an extreme one since it causes a tension with HERA data on the low x behavior, at small Q 2 , of the proton structure func-  \nFigure 5 : Transverse momentum spectra for (mini)jet production in pp collisions at √ s = 10 2 , 10 3 , and 10 4 GeV, as indicated in the plot, as calculated with the default QGSJETIII model (red solid lines) and with the option characterized by twice stronger HT corrections (blue dashed lines).  \n<!-- image -->  \ntion F 2 ( x, Q 2 ) plotted in Fig. 6. Additionally, we vary the parameter α sea which governs the LC momentum distribution of constituent sea (anti)quarks in the model ( ∝ x -α sea ) [13]: using α sea = 0 . 8 and α sea = 0 . 9 , in addition to the default value α sea = 0 . 65 . For all the considered options, we adjust other parameters of  \nFigure 6 : x -dependence of the proton structure function F 2 ( x, Q 2 ) , for different Q 2 , as indicated in the plots, as calculated with the default QGSJET-III model (red solid lines) and with the option characterized by twice stronger HT corrections (blue dashed lines), compared to HERA data [42] (points).  \n<!-- image -->  \nthe hadronization procedure of the model in order to keep a reasonable agreement with hadron production data, both from fixed target experiments and from LHC, as illustrated in Figs. 7 and 8. Noteworthy is the modification of the η -dependence of the charged hadron yield, for increasing α sea (cf. the plot on the right-hand side of Fig. 8): a steeper fall-down of the distribution at large η . However, the experimental accuracy does not allow one to discriminate between the different values of α sea .  \nThe impact of the considered modifications of the interaction treatment on the energy dependence of the inelasticity K inel of proton-nitrogen collisions is shown in Fig. 9, while the corresponding changes of the predicted shower maximum depth are plotted in Fig. 10. As one can see in Figs. 9 and 10, increasing the strength of HT effects, without modifying momentum distributions of constituent partons (the case α sea = 0 . 65 , shown by the black dotted lines in the Figures), does not produce any significant differences, compared to the default QGSJETIII predictions. This is because such HT corrections mostly affect relatively central (small impact parameter b ) collisions involving high parton densities and being characterized by high multiple scattering rates, hence, by a very high inelasticity. The overall reduction of multiple scattering due to stronger HT effects does not  \nhave an appreciable impact on the average K inel which is rather dominated by contributions of more peripheral (larger b ) proton-nucleus collisions.  \nThe situation changes noticeably when choosing softer LC momentum distributions of constituent partons (using α sea = 0 . 8 and α sea = 0 . 9 ), to which strings of color field are connected [see Fig. 11 (left)]. The larger α sea value is used, the shorter such strings are and the smaller fraction of LC momentum of incident proton goes into multiple hadron production resulting from the fragmentation of such strings. As one can see in Figs. 9 and 10, for α sea = 0 . 9 , one has up to /similarequal 6 % reduction of K inel and up to /similarequal 12 g/cm 2 larger 〈 X max 〉 at the highest energies.  \nWhile the considered modifications remain compatible with LHC data on secondary hadron production, as demonstrated in Fig. 8, one may wonder why their impact on the energy dependence of K inel and 〈 X max 〉 is so moderate. This is quite nontrivial, being related to the so-called initial state radiation (ISR) of partons in hard scattering processes. As discussed, e.g. in [47], in any partial semihard rescattering, the hardest (highest p t ) parton scattering process is typically preceded by multiple emission of 'softer' partons characterized by smaller transverse momenta but having larger fractions of LC momentum taken from the parent hadron. In fact, it  \nFigure 7 : x F -dependence of p t -integrated invariant cross section for proton production (left) and x F -distribution of neutrons (right) in c.m. frame, for pp collisions at 158 GeV/c, compared to NA49 data [43]. Black dotted, blue dashed, and green dash-dotted lines correspond to calculations with the modified QGSJET-III model characterized by twice stronger HT corrections, for the parameter α sea = 0 . 65 , 0.8, and 0.9, respectively.  \n<!-- image -->  \nis this ISR, i.e., the perturbative parton cascade preceding the hardest scattering process, which produces large collinear and infrared logarithms compensating the smallness of the strong coupling constant involved and gives rise to a steep energy rise of (mini)jet production [48]. The crucial thing here is the p t and x ordering: each previous parton in the 'ladder' formed by successive parton emissions is typically characterized by a much smaller transverse momentum p t than the next one, while having a much larger LC momentum fraction x . Therefore, the total LC momentum fraction taken by the first s -channel partons produced in such perturbative parton cascades constitutes a natural lower bound on the inelasticity, whatever soft distribution for nonperturbative constituent partons is chosen [47].  \nThus, we arrived here to the point that potential variations of K inel and the corresponding changes of X max are restricted by theoretical arguments, rather than by experimental data. A natural question is how robust are these restrictions. In principle, one can not exclude the possibility that the hadron production pattern corresponding to the above-discussed picture is modified at sufficiently high energies by collective effects. In particular, one popular option is to consider a rearrangement of 'color flows', which gives rise to different string configura-  \ntions, typically reducing both the total string length in the rapidity space and the fractions of LC momenta taken by string end partons from the interacting hadrons (nuclei) 3 (e.g. [50]), thereby overcoming the above-discussed lower bound on K inel .  \nHere we prefer to restrain from relying on a particular treatment of collective effects in secondary hadron production, rather using an effective approach which allows for extreme modifications of the predicted inelasticity. Namely, for each semihard rescattering, we generate the hardest scattering process according to the standard collinear factorization formalism of the perturbative quantum chromodynamics, while suppressing ISR, such that strings of color field are stretched between constituent partons of the interacting hadrons (nuclei) and final state partons produced in the hardest process, as sketched in Fig. 11 (right). In such a scheme, varying the α sea parameter corresponds effectively both to a modification of momentum distributions of constituent partons and to a varying strength of collective effects in the hadronization procedure. In particular, in the limit α sea → 1 , one ends up with very short strings of color field, concentrated in the central  \nFigure 8 : Calculated pseudorapidity η distributions of charged hadrons in c.m. frame. Left: for pp collisions at different √ s (from top to bottom: 13, 7, 2.36, and 0.9 TeV), for hadron transverse momentum p t > 0 . 5 GeV/c, compared to ATLAS data [44, 45]. Right: for pp collisions at √ s = 8 TeV, compared to the data of CMS and TOTEM [46]. The notations for the lines are the same as in Fig. 7.  \n<!-- image -->  \nFigure 9 : Lab. energy dependence of the inelasticity of p 14 N collisions, calculated with the default QGSJET-III model (red solid line) and with the option characterized by twice stronger HT corrections, for the parameter α sea = 0 . 65 , 0.8, and 0.9 - black dotted, blue dashed, and green dash-dotted lines, respectively.  \n<!-- image -->  \nrapidity region in c.m. frame, such that multiple semihard rescattering processes have a minor impact on K inel . However, it is important to remark that such a limit is nonphysical for two reasons. First, potential collective effects may be efficient in relatively central collisions character-  \nFigure 10 : Dependence on primary energy of the average maximum depth of proton-initiated EAS. The notations for the lines are the same as in Fig. 9.  \n<!-- image -->  \nized by high parton densities, while being rather weak in more peripheral collisions dominating the average inelasticity. Secondly, the chance that collective effects, however strong they are, eliminate all the partons from ISR should be vanishingly small.  \nAs previously, we perform a modeling of hadron-proton and hadron-nucleus interactions, using different values of α sea and adjusting other  \n-  \nFigure 11 : Schematic view of a single semihard scattering process. In the standard treatment (left), strings are formed between constituent partons (quarks and antiquarks) and/or partons produced by perturbative cascades, following the color and anticolor flows (thick green and blue lines). Neglecting ISR (right), strings are formed between constituent partons and partons emerging from the hardest scattering process.  \n<!-- image -->  \n<!-- image -->  \n-  \nparameters of the hadronization procedure in order to keep an agreement with accelerator measurements, a comparison with selected data sets being plotted in Figs. 12 and 13. As we can see in Fig. 14, even using the default value α sea = 0 . 65 , the inelasticity is reduced by up to /similarequal 10 % at the highest energies, compared to the QGSJET-III predictions, which demonstrates the importance of perturbative parton cascades both for hadron production in general and for the energy loss of leading nucleons in particular. On the other hand, for α sea = 0 . 9 , K inel is practically energy-independent above 1 PeV, where secondary hadron production is dominated by semihard processes. The corresponding energy dependence of the predicted 〈 X max 〉 , for α sea = 0 . 65 , 0.8, and 0.9, is plotted in Fig. 15. The reduction of the inelasticity leads to a noticeably larger elongation rate for proton-induced EAS; for α sea = 0 . 9 , the shower maximum is /similarequal 30 g/cm 2 deeper at E 0 = 10 20 eV, compared to the QGSJET-III predictions.  \nNow the crucial question is whether such modifications can be constrained by LHC data on forward hadron production. In Fig. 16, we compare production spectra of neutrons, calculated with the modified QGSJET-III model: suppressing ISR and using α sea = 0 . 65 , 0.8, and 0.9, to the corresponding measurements of the LHCf  \nexperiment. As one can see in Fig. 16, choosing softer momentum distributions of constituent partons, i.e., using larger α sea , one obtains larger neutron yields at high x F /similarequal 2 E/ √ s /greaterorsimilar 0 . 5 than observed by the experiment, thereby underestimating the inelasticity for neutron production. 4 Thus, the most extreme modifications of the model, leading to an approximate Feynman scaling in the fragmentation region, are somewhat disfavored by the LHCf data. It is noteworthy that more stringent constraints may arise from studying correlations between central and forward hadron production [19], e.g., from the ongoing combined measurements of hadron production in pp collisions by the ATLAS and LHCf experiments [52, 53].  \nIt is worth remarking that the modifications of the interaction treatment, leading to a substantially larger elongation rate, are rather strongly disfavored by observations of the Pierre Auger Observatory, regarding the muon production depth in air showers [54, 55], notably, by the measured values of 〈 X µ max 〉 - the average depth of maximum of the muon production profile [56]. As discussed in [5, 57], changes of a model treatment of proton-air collisions, which produce a larger 〈 X max 〉 , shift the average depth of maximum of the muon production profile deeper in the atmosphere by a comparable amount. Moreover, since the above-considered modifications, namely, the suppression of ISR and the use of softer momentum distributions of constituent partons, impact also pion-air interactions, the corresponding effect on 〈 X µ max 〉 is stronger, compared to the change of 〈 X max 〉 , as demonstrated in Fig. 17. 5 For the most extreme modifications, the calculated 〈 X µ max 〉 appears to be up  \n5 Plotted in Fig. 17 is the the average depth of maximum of muon production profiles produced by EAS simulations, for muon energies E µ > 1 GeV. This should be distinguished from the so-called apparent muon production depth derived by taking into consideration both the impact of muon propagation in the atmosphere and the effects of the corresponding experimental reconstruction procedures, regarding measurements with ground-based detectors (see, e.g. [54, 55] for the corresponding discussion).  \nFigure 12 : x F -dependence of p t -integrated invariant cross section for proton production (left) and x F -distribution of neutrons (right) in c.m. frame, for pp collisions at 158 GeV/c, compared to NA49 data [43]. Black dotted, blue dashed, and green dash-dotted lines correspond to calculations with the modified QGSJET-III model: suppressing ISR and using α sea = 0 . 65 , 0.8, and 0.9, respectively.  \n<!-- image -->  \nto /similarequal 40 g/cm 2 larger at the highest energies, than predicted by QGSJET-III, thereby creating a strong tension with the Auger data. 6", '4 Summary': 'We performed a quantitative analysis of model uncertainties for predicted maximum depth of proton-initiated extensive air showers, in the framework of the QGSJET-III hadronic interaction model, restricting ourselves to the standard physics picture. Using the conventional approach to the treatment of high energy interactions, we investigated a possibility to obtain larger values of 〈 X max 〉 , considering variations of the inelastic proton-proton cross section, the rate of inelastic diffraction, the strength of nonlinear interaction effects, and momentum distributions of constituent partons involved in multiple scattering processes, allowed by LHC data. The studied modifications of the interaction treatment allowed us to increase the predicted 〈 X max 〉 by only ∼ 10 g/cm 2 .  \nSuch a small variation of the predicted EAS maximum depth, when modifying σ inel pp and σ diffr pp within the range allowed by accelerator data, comes at no surprise, given extensive and precise measurements of proton-proton interaction cross sections at LHC. Yet one could have expected larger changes of the calculated 〈 X max 〉 to result from the other options studied, because of their potentially strong impact on the inelasticity of proton-proton and proton-nucleus collisions. However, potential variations of the inelasticity appeared to be limited by the initial state radiation of partons in semihard scattering processes: the total fraction of the incident hadron momentum, taken by all such perturbatively generated partons constitutes a natural lower bound on the inelasticity.  \nWe further investigated a more exotic scenario, considering a potentially significant modification of the parton hadronization procedure by hypothetical collective effects. That way, we were able to change drastically the predicted energy dependence of the inelasticity of protonair collisions and to increase thereby the predicted EAS maximum depth by up to /similarequal 30 g/cm 2 . However, those most extreme modifications appeared to be disfavored both by the data of the LHCf experiment, regarding forward neutron production in pp collisions at LHC, and by measurements of the muon production depth by the Pierre Auger Observatory.  \nFigure 13 : Calculated pseudorapidity η distributions of charged hadrons in c.m. frame. Left: for pp collisions at different √ s (from top to bottom: 13, 7, 2.36, and 0.9 TeV), for hadron transverse momentum p t > 0 . 5 GeV/c, compared to ATLAS data [44, 45]. Right: for pp collisions at √ s = 8 TeV, compared to the data of CMS and TOTEM [46]. The notations for the lines are the same as in Fig. 12.  \n<!-- image -->  \nFigure 14 : Lab. energy dependence of the inelasticity of p 14 N collisions, calculated with the default QGSJET-III (red solid line) and with the modified model: suppressing ISR and using α sea = 0 . 65 , 0.8, and 0.9 - black dotted, blue dashed, and green dash-dotted lines, respectively.  \n<!-- image -->', 'Acknowledgments': 'The work of S.O. was supported by Deutsche Forschungsgemeinschaft (project number 465275045). G.S. acknowledges support by the Bundesministerium für Bildung und Forschung, under grants 05A20GU2 and 05A23GU3.  \nFigure 15 : Dependence on primary energy of the average maximum depth of proton-initiated EAS. The notations for the lines are the same as in Fig. 14.  \n<!-- image -->', 'References': '- [1] M. Nagano and A. A. Watson, Observations and implications of the ultrahighenergy cosmic rays , Rev. Mod. Phys. 72 , 689 (2000).\n- [2] D. Heck, J. Knapp, J. N. Capdevielle, G. Schatz, and T. Thouw, CORSIKA: A Monte Carlo code to simulate extensive air  \nFigure 16 : Calculated neutron energy spectra in c.m. frame, for pp collisions at √ s = 13 TeV, compared to the data of the LHCf experiment [51] (points). Shown as red solid histograms are the results of the default QGSJET-III, while back dotted, blue dashed, and green dash-dotted histograms correspond to calculations with the modified model: suppressing ISR and using α sea = 0 . 65 , 0.8, and 0.9, respectively.  \n<!-- image -->  \nshowers , Forschungszentrum Karlsruhe Internal Report FZKA-6019 (1998).  \n- [3] R. Engel, D. Heck, and T. Pierog, Extensive air showers and hadronic interactions at high energy , Ann. Rev. Nucl. Part. Sci. 61 , 467 (2011).\n- [4] K.-H. Kampert and M. Unger, Measurements of the Cosmic Ray Composition with Air Shower Experiments , Astropart. Phys. 35 , 660 (2012).\n- [5] S. Ostapchenko, High energy interactions of cosmic rays , Adv. Space Res. 64 , 2445 (2019).\n- [6] P. Abreu et al. (Pierre Auger Collaboration), Interpretation of the Depths of Maximum of Extensive Air Showers Measured  \n- by the Pierre Auger Observatory , JCAP 02 , 026 (2013).\n- [7] A. Abdul Halim et al. (Pierre Auger Collaboration), Testing Hadronic-Model Predictions of Depth of Maximum of Air-Shower Profiles and Ground-Particle Signals using Hybrid Data of the Pierre Auger Observatory , Phys. Rev. D 109 , 102001 (2024).\n- [8] A. Aab et al. (Pierre Auger Collaboration), Depth of maximum of air-shower profiles at the Pierre Auger Observatory. I. Measurements at energies above 10 17 . 8 eV , Phys. Rev. D 90 , 122005 (2014).\n- [9] A. Aab et al. (Pierre Auger Collaboration), Depth of maximum of air-shower profiles at the Pierre Auger Observatory. II. Composi-  \nFigure 17 : Dependence on primary energy of the average depth of maximum of the muon production profile, for proton-initiated EAS, as calculated with the default QGSJET-III (red solid line) and with the modified model: suppressing ISR and using α sea = 0 . 65 , 0.8, and 0.9 - black dotted, blue dashed, and green dashdotted lines, respectively.  \n<!-- image -->  \n- tion implications , Phys. Rev. D 90 , 122006 (2014).\n- [10] S. Ostapchenko, Monte Carlo treatment of hadronic interactions in enhanced Pomeron scheme: QGSJET-II model , Phys. Rev. D 83 , 014018 (2011).\n- [11] S. Ostapchenko, QGSJET-II: physics, recent improvements, and results for air showers , EPJ Web Conf. 52 , 02001 (2013).\n- [12] S. Ostapchenko, QGSJET-III model of high energy hadronic interactions: The formalism , Phys. Rev. D 109 , 034002 (2024).\n- [13] S. Ostapchenko, QGSJET-III model of high energy hadronic interactions: II. Particle production and extensive air shower characteristics , Phys. Rev. D 109 , 094019 (2024).\n- [14] S. Ostapchenko and G. Sigl, On the model uncertainties for the predicted muon content of extensive air showers , Astropart. Phys. 163 , 103004 (2024).\n- [15] G. Antchev et al. (TOTEM Collaboration), First measurement of elastic, inelastic and total cross-section at √ s = 13 TeV by\n- TOTEM and overview of cross-section data at LHC energies , Eur. Phys. J. C 79 , 103 (2019).\n- [16] G. Aad et al. (ATLAS Collaboration), Measurement of the total cross section and ρ -parameter from elastic scattering in pp collisions at √ s = 13 TeV with the ATLAS detector , Eur. Phys. J. C 83 , 441 (2023).\n- [17] R. J. Glauber, High-energy collision theory , in: Lectures in theoretical physics, Ed. by W. E. Brittin and L. G. Dunham, Interscience Publishers, New York, 1959, vol. 1, p. 315.\n- [18] V. N. Gribov, Glauber corrections and the interaction between high-energy hadrons and nuclei , Sov. Phys. JETP 29 , 483 (1969).\n- [19] S. Ostapchenko, M. Bleicher, T. Pierog, and K. Werner, Constraining high energy interaction mechanisms by studying forward hadron production at the LHC , Phys. Rev. D 94 , 114026 (2016).\n- [20] R. Aloisio, V. Berezinsky, P. Blasi, and S. Ostapchenko, Signatures of the transition from Galactic to extragalactic cosmic rays , Phys. Rev. D 77 , 025007 (2008).\n- [21] G. Antchev et al. (TOTEM Collaboration), Luminosity-independent measurements of total, elastic and inelastic cross-sections at √ s = 7 TeV , Europhys. Lett. 101 , 21004 (2013).\n- [22] G. Antchev et al. (TOTEM Collaboration), Measurement of proton-proton inelastic scattering cross-section at √ s = 7 TeV , Europhys. Lett. 101 , 21003 (2013).\n- [23] G. Antchev et al. (TOTEM Collaboration), Luminosity-Independent Measurement of the Proton-Proton Total Cross Section at √ s = 8 TeV , Phys. Rev. Lett. 111 , 012001 (2013).\n- [24] G. Aad et al. (ATLAS Collaboration), Measurement of the total cross section from elastic scattering in pp collisions at √ s = 7 TeV with the ATLAS detector , Nucl. Phys. B 889 , 486 (2014).  \n- [25] M. Aaboud et al. (ATLAS Collaboration), Measurement of the total cross section from elastic scattering in pp collisions at √ s = 8 TeV with the ATLAS detector , Phys. Lett. B 761 , 158 (2016).\n- [26] G. Antchev et al. (TOTEM Collaboration), Double diffractive cross-section measurement in the forward region at the LHC , Phys. Rev. Lett. 111 , 262001 (2013).\n- [27] G. Aad et al. (ATLAS Collaboration), Measurement of differential cross sections for single diffractive dissociation in √ s = 8 TeV pp collisions using the ATLAS ALFA spectrometer , J. High Energy Phys. 02 , 042 (2020).\n- [28] F. Oljemark, Single Diffraction in protonproton scattering with TOTEM at the Large Hadron Collider , PhD thesis, University of Helsinki (2020).\n- [29] S. Ostapchenko, LHC data on inelastic diffraction and uncertainties in the predictions for longitudinal extensive air shower development , Phys. Rev. D 89 , 074009 (2014).\n- [30] J. Engel, T. K. Gaisser, T. Stanev, and P. Lipari, Nucleus-nucleus collisions and interpretation of cosmic ray cascades , Phys. Rev. D 46 , 5013 (1992).\n- [31] N. N. Kalmykov and S. S. Ostapchenko, The nucleus-nucleus interaction, nuclear fragmentation, and fluctuations of extensive air showers , Phys. Atom. Nucl. 56 , 346 (1993).\n- [32] A. Białas, M. Bleszynski, and W. Czyz, Multiplicity Distributions in NucleusNucleus Collisions at High-Energies , Nucl. Phys. B 111 , 461 (1976).\n- [33] N. N. Kalmykov and S. S. Ostapchenko, Comparison of Nucleus-Nucleus Interaction Characteristics in the Model of QuarkGluon Strings and in the Superposition Model , Sov. J. Nucl. Phys. 50 , 315 (1989).\n- [34] S. Fredriksson, G. Eilam, G. Berlad, and L. Bergström, High-energy Collisions With Atomic Nuclei. Part 1 , Phys. Rept. 144 , 187 (1987).\n- [35] M. L. Good and W. D. Walker, Diffraction disssociation of beam particles , Phys. Rev. 120 , 1857 (1960).\n- [36] V. A. Khoze, A. D. Martin, and M. G. Ryskin, Dynamics of diffractive dissociation , Eur. Phys. J. C 81 , 175 (2021).\n- [37] S. Ostapchenko and M. Bleicher, Taming the energy rise of the total proton-proton cross-section , Universe 5 , 106 (2019).\n- [38] R. L. Workman et al. (Particle Data Group), Review of Particle Physics , Prog. Theor. Exp. Phys. 2022 , 083C01 (2022).\n- [39] T. Bergmann, R. Engel, D. Heck, N. N. Kalmykov, S. Ostapchenko, T. Pierog, T. Thouw, and K. Werner, One-dimensional Hybrid Approach to Extensive Air Shower Simulation , Astropart. Phys. 26 , 420 (2007).\n- [40] S. S. Ostapchenko, Contemporary models of high-energy interactions: Present status and perspectives , J. Phys. G 29 , 831 (2003).\n- [41] R. D. Parsons, C. Bleve, S. S. Ostapchenko, and J. Knapp, Systematic uncertainties in air shower measurements from highenergy hadronic interaction models , Astropart. Phys. 34 , 832 (2011).\n- [42] F. D. Aaron et al. (H1 and ZEUS Collaborations), Combined measurement and QCD analysis of the inclusive e ± p scattering cross sections at HERA , J. High Energy Phys. 01 , 109 (2010).\n- [43] T. Anticic et al. (NA49 Collaboration), Inclusive production of protons, anti-protons and neutrons in p+p collisions at 158GeV/c beam momentum , Eur. Phys. J. C 65 , 9 (2010).\n- [44] G. Aad et al. (ATLAS Collaboration), Charged-particle multiplicities in pp interactions measured with the ATLAS detector at the LHC , New J. Phys. 13 , 053033 (2011).\n- [45] G. Aad et al. (ATLAS Collaboration), Charged-particle distributions in √ s = 13 TeV pp interactions measured with the ATLAS detector at the LHC , Phys. Lett. B 758 , 67 (2016).  \n- [46] S. Chatrchyan et al. (CMS and TOTEM Collaborations), Measurement of pseudorapidity distributions of charged particles in proton-proton collisions at √ s = 8 TeV by the CMS and TOTEM experiments , Eur. Phys. J. C 74 , 3053 (2014).\n- [47] S. Ostapchenko, Cosmic ray interactions in the atmosphere: QGSJET-III and other models , SciPost Phys. Proc. 13 , 004 (2023).\n- [48] L. V. Gribov, E. M. Levin, and M. G. Ryskin, Semihard Processes in QCD , Phys. Rept. 100 , 1 (1983).\n- [49] K. Werner, Core-corona procedure and microcanonical hadronization to understand strangeness enhancement in proton-proton and heavy ion collisions in the EPOS4 framework , Phys. Rev. C 109 , 014910 (2024).\n- [50] J. R. Christiansen and P. Z. Skands, String Formation Beyond Leading Colour , J. High Energy Phys. 08 , 003 (2015).\n- [51] O. Adriani et al. (LHCf Collaboration), Measurement of inclusive forward neutron production cross section in proton-proton collisions at √ s = 13 TeV with the LHCf Arm2 detector , J. High Energy Phys. 11 , 073 (2018).\n- [52] ATLAS and LHCf Collaborations, Measurement of contributions of diffractive processes to forward photon spectra in pp collisions at √ s = 13 TeV , Tech. Rep. ATLASCONF-2017-075, CERN, Geneva, 2017.\n- [53] M. Kondo et al. (LHCf Collaboration), Performance evaluation of LHCf-ATLAS ZDC joint measurement using proton beam , EPJ Web Conf. 283 , 05012 (2023).\n- [54] S. Andringa, L. Cazon, R. Conceicao, and M. Pimenta, The Muonic longitudinal shower profiles at production , Astropart. Phys. 35 , 821 (2012).\n- [55] L. Cazon, R. Conceicao, M. Pimenta, and E. Santos, A model for the transport of muons in extensive air showers , Astropart. Phys. 36 , 211 (2012).  \n- [56] A. Aab et al. (Pierre Auger Collaboration), Muons in Air Showers at the Pierre Auger Observatory: Measurement of Atmospheric Production Depth , Phys. Rev. D 90 , 012012 (2014).\n- [57] S. Ostapchenko and M. Bleicher, Constraining pion interactions at very high energies by cosmic ray data , Phys. Rev. D 93 , 051501(R) (2016).'}