Constraints on high energy interaction models from LHC and cosmic ray data

Predictions of popular cosmic ray interaction models for some basic characteristics of cosmic ray-induced extensive air showers are analyzed in view of experimental data on proton-proton collisions, obtained at the Large Hadron Collider. The differences between the results are traced down to different approaches for the treatment of hadronic interactions, implemented in those models. Potential measurements by LHC and cosmic ray experiments, which could be able to discriminate between the alternative approaches, are proposed.


Introduction
Modeling of high energy hadronic interactions is of considerable importance for experimental studies at the Large Hadron Collider (LHC) and, especially, in high energy cosmic ray (CR) field.In the latter case, one traditionally relies on the extensive air shower (EAS) techniques: the properties of primary CR particles are reconstructed from measured characteristics of nuclear-electro-magnetic cascades induced by CR interactions in the atmosphere.This naturally implies the importance of detailed Monte Carlo simulations of EAS development, particularly, of its backbone -the cascade of nuclear interactions of both the primary particles and of secondary hadrons produced.Thus, the success of the experimental studies depends on the validity of hadronic interaction models used in the analysis.
One usually chooses between two main experimental procedures [1].In the first case, one deals with the information obtained with scintillation detectors positioned at ground.The energy of the primary particle is reconstructed from the measured lateral density of charged particles (mostly electrons and positrons) while the particle type is inferred from the relative fraction of muons, compared to all charged particles at ground.Alternatively, one may study the longitudinal EAS development by measuring fluorescence light produced by excited air molecules at different heights in the atmosphere, for which purpose dedicated fluorescence telescopes are employed.The primary energy is then related to the total amount of fluorescence light emitted.In turn, the particle type may be determined from the measured position of the shower maximum X max -the depth in the atmosphere (in g/cm 2 ) where the number of ionizing particles reaches its maximal value.
Not surprisingly, the observables used to determine the primary particle type -the lateral muon density ρ μ and the EAS maximum position X max -appear to be very sensitive to details of high energy hadronic interactions.More precisely, X max depends strongly on the properties of the primary particle √ s-dependence of the total, inelastic, and elastic pp cross sections, as calculated using the QGSJET-II-04 [5], EPOS-LHC [3], SIBYLL-2.3[7], and QGSJET [8] models (solid, dashed, dash-dotted, and dotted lines respectively).Experimental data are from Refs.[10][11][12].
interaction with air nuclei: the inelastic cross section and the forward spectra of secondary hadrons produced.In turn, the EAS muon content is formed in a multi-step cascade process, driven mostly by interactions of secondary pions and kaons with air.Hence, ρ μ depends strongly on the properties of pion-air collisions, most importantly, on the multiplicity and spectral shape for charged hadrons produced.
In the following, we shall mainly deal with the most important and best-measured EAS observable -the shower maximum position.The corresponding EAS simulation results are most directly related to the modeling of very high energy proton-proton collisions, hence, are expected to be strongly constrained by experimental studies at LHC.We shall compare the results of X max calculations, obtained using the most recent versions of the EPOS [2,3], QGSJET-II [4,5], and SIBYLL [6,7] models, which have been updated using experimental data from Run 1 of LHC.Additionally, we shall use the QGSJET model [8] which, though being already outdated physics-wise, demonstrated a generally good agreement with the data from Run 1 of LHC on minimum-bias collisions [9].Hence, it may be used to study the range of potential variations of model predictions for X max , in view of current LHC data.Our primary goals are to analyze the differences between the model results, to trace their origin to the underlying approaches for the treatment of hadronic collisions, implemented in those models, and to propose potential measurements by LHC and cosmic ray experiments, which could be able to discriminate between the alternative approaches.

Model predictions for X max
The uncertainties of model predictions for X max have been greatly reduced thanks to the LHC studies of very high energy proton-proton collisions.Of particular importance are precise measurements of the total, inelastic, and elastic pp cross sections by the TOTEM and ATLAS experiments [10,11], which strongly constrained the corresponding model results in the ultra-high energy limit, as illustrated in Fig. 1.In turn, this contributes to the convergence of the model predictions for the inelastic proton-air and nucleus-air cross sections which determine the depth on the first interaction of the primary particle in the atmosphere.It is worth also stressing the importance of the measurements of the differential elastic pp cross section which provides information on the spatial parton distribution in proton.The latter impacts both the calculations of the proton-air cross section and the magnitude of nonlinear interaction effects in proton-proton and proton-nucleus collisions.
Yet there remain considerable differences between the model predictions for X max , 1 as illustrated in Fig. 2. One might relate those to present experimental uncertainties concerning the rate of the inelastic diffraction in high energy pp collisions.Indeed, diffraction has a significant influence on the shape of very forward spectra of secondary hadrons, which in turn makes a strong impact on the longitudinal EAS development.This has been investigated in Ref. [14] in the framework of the QGSJET-II-04 model; the obtained characteristic uncertainty for X max amounted to 10 g/cm 2 , being thus considerably smaller than the differences between the results plotted in Fig. 2. One might then ask if the analysis of Ref. [14] was not general enough or X max depends on some other characteristics of hadronic interactions, not well constrained by present LHC data.

Impact of constituent parton Fock states
An important insight into the problem is provided by recent combined measurements by the CMS and TOTEM experiments of the pseudorapidity η density dn ch pp /dη of produced charged hadrons in pp collisions at √ s = 8 TeV [15].In contrast to the EPOS-LHC and QGSJET-II-04 models, hadronic Monte Carlo generators used in the collider field turn to underestimate considerably secondary hadron production at large η.As illustrated in Fig. 3, a similar problem arises for the SIBYLL-2.3model.This appears to be related to model assumptions concerning the structure of constituent parton Fock states in hadrons, as discussed in more detail in Ref. [16].Let us consider a proton-proton interaction in the center of mass (c.m.) frame, which involves multiple binary interactions of partons from two protons' parton clouds formed prior to the collision.We are interested here in the momentum distribution of these multiparton states, which is related to the evolution of parton cascades in hadrons.The picture behind the corresponding treatment of the SIBYLL model and of event generators used 1 Here and in the following the calculations of extensive air shower development are performed using the CONEX code [13].
TeV, as calculated using the QGSJET-II-04, EPOS-LHC, and SIBYLL-2.3models (solid, dashed, and dash-dotted lines respectively) for the nondiffractive event selection of TOTEM: at least one charged hadron produced both at −6.5 < η < −5.3 and at 5.3 < η < 6.5. in the collider field is shown schematically on the left-hand side (lhs) of Fig. 4. At large Feynmanx, one usually starts from the same universal parton Fock state.Additional partons (sea quarks and gluons) giving rise to new branches of the parton cascade, which take part in the multiple scattering processes, result from the evolution of the parton density corresponding to this initial state and their momentum fractions are typically distributed as ∝ 1/x in the very high energy limit.Such a picture reflects itself in the hadron production pattern predicted by such models: multiple scattering mostly affects central particle production, while having a weak influence of forward hadron spectra.Indeed, the latter are formed by the hadronization of partons emerging from the initial part of the underlying parton cascade, which starts from the same initial conditions and covers a short rapidity interval, being thus weakly dependent on the further development of the cascade.In particular, the steep fall-down of dn ch pp /dη predicted by SIBYLL-2.3(dash-dotted line in Fig. 3) for large η values reflects the quick decrease of the number of constituent partons when parton momentum fraction increases.In the alternative approach, implemented in the EPOS and QGSJET(-II) models, a proton is represented by a superposition of a number of Fock states containing different numbers of large x constituent partons, as shown schematically on the right-hand side (rhs) of Fig. 4. Further cascading of these partons "dresses" them with low x parton clouds.As the overall parton multiplicity in the central rapidity region is roughly proportional to the number of initial constituent partons, stronger multiple scattering is typically associated with larger Fock states.Thus, there is a strong long-range correlation between central and forward particle production; higher multiplicity in the central rapidity region reflects stronger multiple parton scattering.In turn, this implies that bigger numbers of large x constituent partons are involved in the process, which has a strong impact on forward hadron spectra.
The two approaches give rise to quite different predictions for the energy dependence of the inelasticity K inel , i.e. the relative energy loss of leading nucleons, in proton-proton and proton-nucleus collisions.Both schemes predict an enhancement of multiple scattering with the increase of collision energy -as a consequence of the enlarged kinematic space and of an increase of parton densities in the low x limit.However, in the approach corresponding to the schematic picture in the lhs of Fig. 4, this results in a substantial enhancement of secondary particle production in the central rapidity region only, while having a weak impact on forward hadron spectra.Hence, one expects a weak energy dependence of K inel in that case, which we observe indeed for SIBYLL 2.3 in Fig. 5, where its prediction for the inelasticity K inel pp for pp collisions is shown by dash-dotted line.In the alternative approach corresponding to the picture in the rhs of Fig. 4, one obtains a substantial energy-rise of the inelasticity, which is clearly seen in Fig. 5 for QGSJET-II-04 and EPOS-LHC (solid and dashed lines, respectively).The reason for this rise is twofold.First, for any given Fock state, increasing multiple scattering implies that bigger numbers of large x constituent partons are involved in the interaction, thus leaving smaller fractions of the initial proton momentum for spectator partons which finally form the "leading" (most energetic) nucleons.Additionally, Fock states with bigger and bigger numbers of large x constituent partons come into play.The momentum sharing between these partons gives rise to a smaller fraction of the initial proton momentum, possessed by each parton, which thus enhances the energy loss of the leading nucleons.This analysis clearly explains the larger X max values predicted by SIBYLL 2.3 (dash-dotted lines in Fig. 2) -as a smaller inelasticity results in a slower EAS development (see, e.g.Ref. [17]).As is evident from the discussion above, the best way to discriminate between the two theoretical approaches depicted schematically in Fig. 4 is via studies of correlations of the signal strength in central and forward-looking detectors at LHC, which may be performed by the CMS and TOTEM experiments.This is illustrated in Fig. 6, where we plot for √ s = 8 TeV the η-density of charged hadrons dn ch pp /d|η| (for transverse momentum p t > 0) at |η| = 6 (averaged over the interval 5.5 < |η| < 6.5) as a function of the central η-density of charged hadrons of p t > 0.1 GeV (averaged over the interval |η| < 1).Both EPOS-LHC and QGSJET-II-04 predict a strong correlation of the signal strength in CMS and TOTEM.The respective results of the two models practically coincide with each other, apart from the region corresponding to the tails of the multiplicity distributions.In contrast, for SIBYLL 2.3 such a correlation is twice weaker, being thus a "smoking gun" signature for the desirable discrimination.
Another possible way for the model discrimination is via measurements of very forward particle production, e.g. by the LHCf experiment at LHC, when supplemented by triggering different hadronic activities in the central detectors of the ATLAS experiment, as discussed in Ref. [16].

Relevance of pion-air interactions
The analysis in Section 3 does not explain the differences in X max predictions of EPOS-LHC, QGSJET-II-04, and QGSJET -as all the three models employ essentially the same treatment of constituent parton Fock states.To reveal the interaction features which are responsible for these differences, we employ the "cocktail" model approach: using QGSJET-II-04 to describe some selected interactions of hadrons in air showers or some particular features of the primary interaction, while treating the rest with one of the other two models (see Ref. [18] for more details).As the first step, we apply QGSJET-II-04 to determine the position of the primary particle interaction in the atmosphere and to describe the production of secondary nucleons in this interaction; all other characteristics of the first proton-air collision and all the subsequent interactions of secondary hadrons in the cascade are treated using EPOS-LHC.This way we check the sensitivity of the calculated shower maximum depth to the model differences concerning the proton-air cross section and the predicted nucleon spectra, which thus comprise the effects of the inelastic diffraction.The calculated X max values shown by the blue dash-dotted line in Fig. 7 differ from the original EPOS-LHC results by not more than 5 g/cm 2 , which is well within the uncertainty range obtained in Ref. [14].
Next, we apply QGSJET-II-04 to describe all the characteristics of the primary interaction, while treating the rest of the hadron cascade using EPOS-LHC.The obtained X max shown by the blue dashed line in Fig. 7 is shifted further towards the QGSJET-II-04 results by up to 5 g/cm 2 .This additional shift is explained by somewhat harder spectra of secondary mesons, most importantly, of secondary pions in EPOS-LHC, compared to the QGSJET-II-04 model.We also repeat the same calculation describing secondary hadron interactions in the cascade with the QGSJET model, the results being plotted by the green dashed line in Fig. 7.In this case, the difference with the pure QGSJET-based calculation does not exceed 3 g/cm 2 , which is due to the fact that forward particle spectra in proton-air collisions are rather similar in QGSJET and QGSJET-II-04.
There remains a large difference between the two dashed lines in Fig. 7 and the results of QGSJET-II-04, which arises from the model treatments of pion-air and kaon-air interactions.In the particular case of QGSJET, this is mainly related to larger pion-air cross section and softer production spectra for secondary mesons, predicted by that model, compared to QGSJET-II-04.The sensitivity of X max to the difference in pion-air cross sections is illustrated by the green dash-dotted line in Fig. 7, obtained by using QGSJET-II-04 results both for the primary interaction and for the inelastic cross sections for all the secondary hadron-air collisions.In turn, applying QGSJET-II-04 to describe in addition pion and kaon spectra in pion-air collisions gives rise to the energy-dependence of X max , shown by the green dotted line in Fig. 7.
In turn, for EPOS-LHC the remaining difference with the QGSJET-II-04 results is due to a copious production of baryon-antibaryon pairs in pion-air and kaon-air collisions, also due to harder (anti-)baryon spectra in EPOS-LHC [19].This slows down the energy dissipation from the hadronic cascade and thus contributes to the elongation of the shower profile.Indeed, if we apply QGSJET-II-04 to describe both the primary interaction and the production of nucleons and antinucleons in all the secondary pion-air and kaon-air collisions, while treating the rest with EPOS-LHC, the obtained X max shown by the blue dotted line practically coincides with the QGSJET-II-04 results.

QUARKS-2016 5 Relation to muon production depth
As demonstrated in Section 4, a large part of the model uncertainty for the predicted X max is related to the treatment of pion-air collisions at very high energies.While there exist no experimental data for such interactions above fixed-target energies, the corresponding model treatment may be constrained indirectly by studying other characteristics of cosmic ray-induced air showers.Recently, the Pierre Auger collaboration reported measurements of the maximal muon production depth in air showers, X μ max -the depth in the atmosphere (in g/cm 2 ) where the rate of muon production via decays of pions and kaons reaches its maximal value [20].In particular, one observed a strong contradiction between the results of EAS simulations with EPOS-LHC and the experimental data: the maximum of the muon production profile was observed substantially higher in the atmosphere than predicted by that model for the heaviest possible primary cosmic ray nuclei.
There are both similarities and differences concerning the relation of X max and X μ max to the properties of hadron-air interactions.Obviously, both characteristics are sensitive to the position X 0 of the primary particle interaction in the atmosphere, which depends on the respective inelastic cross section: fluctuations of X 0 shift the whole cascade upwards and downwards in the atmosphere and thus do so for X max and X μ max for a particular shower.However, in contrast to X max , the maximal muon production depth is much less sensitive to hadron production in the primary interaction.The EAS muon content rather depends on the multi-step hadronic cascade in which the number of pions and kaons increases in an avalanche way until the probability of their decays becomes comparable to the one for interactions.This happens when the energies for most of pions approach the pion critical energy, E π ± crit 80 GeV [21].The position of the maximum of the muon production profile is actually close to this turning point.
As a consequence, X μ max is very sensitive to the forward spectral shape of secondary mesons in pion-air collisions: producing in each cascade step a meson of a slightly higher energy would mean that a larger number of cascade branching steps is required for reaching the critical energy, with the result that the maximum of the muon production profile will be observed deeper in the atmosphere.A similar effect may be produced by a smaller pion-air cross section as this would increase the pion mean free pass and thereby elongate the muon production profile.In addition, in the particular case of the EPOS model, its predictions for X μ max may be influenced by the copious production of (anti-)baryons in pion-air collisions.Indeed, unlike pions and kaons, (anti-)nucleons continue to take part in the hadronic cascade down to the GeV energy range, producing additional generations of secondary hadrons.Muons emerging from decays of secondary pions and kaons created in interactions of low energy nucleons and antinucleons contribute to the elongation of the muon production profile and give rise to larger values of X μ max .Performing calculations of X μ max for muon energies E μ ≥ 1 GeV for proton-induced EAS, we observe substantially larger differences between the results of the different interaction models, compared to the case of X max , as demonstrated in Fig. 8.To reveal the reasons for these differences, we use the same "cocktail" model approach as in Section 4. First, we apply QGSJET-II-04 to describe all the characteristics of the primary interaction, while treating the rest of the hadron cascade using either EPOS-LHC or QGSJET, the results shown respectively by the blue and green dashed lines in Fig. 8.As expected, the obtained X μ max values deviate only slightly from the original model calculations -as the differences between the model predictions for X μ max are mainly due to secondary (mostly pion-air) interactions in the cascade.For example, the smaller X μ max values predicted by QGSJET are mostly due to somewhat larger inelastic pion-air and kaon-air cross sections and softer meson spectra produced by that model, compared to QGSJET-II-04.The first effect is illustrated by the green dash-dotted line in Fig. 8, obtained by applying QGSJET-II-04 to describe both the primary interaction and the inelastic cross sections for all the secondary hadron-air collisions, while treating hadron production max (E μ ≥ 1 GeV) for p-induced vertical EAS, as calculated using the EPOS-LHC, QGSJET-II-04, and QGSJET models (top blue, middle red, and bottom green solid lines respectively), or applying mixed model descriptions, as explained in the text (dashed, dash-dotted, and dotted lines).
in the secondary interactions with QGSJET.On the other hand, using QGSJET-II-04 results also for the pion and kaon spectra in pion-air collisions results in the energy-dependence of X μ max , shown by the green dotted line in Fig. 8.
In turn, a large part of the difference between EPOS-LHC and QGSJET-II-04 is due to the copious production of baryon-antibaryon pairs in the former model.This is illustrated by the blue dash-dotted line in Fig. 8, obtained by applying QGSJET-II-04 to describe both the primary interaction and the production of nucleons and antinucleons in all the secondary pion-air collisions, while treating the rest with EPOS-LHC.The remaining difference between the two models is due to harder spectra of secondary mesons in EPOS-LHC for pion-air and kaon-air collisions.Indeed, using QGSJET-II-04 results both for the primary interaction and for hadron spectra in secondary pion-air and kaon-air collisions, we get the energy-dependence of X μ max , shown by the blue dotted line in Fig. 8, which is very close to the pure QGSJET-II-04 calculation.

Outlook
Experimental studies of proton-proton collisions at LHC furnished valuable information for the CR field, substantially reducing the range of uncertainties for numerical simulations of CR-induced air showers.Nevertheless, there exists yet a significant spread between the results of calculations for the basic EAS observable -the shower maximum position.While a part of this spread may be explained by the present experimental uncertainty concerning the diffraction rate in very high energy pp collisions, a much bigger difference between certain model predictions for X max is related to the model assumptions concerning the momentum structure of constituent parton Fock states in hadrons.The "smoking gun" signature for discriminating between the different theoretical approaches is the correlation of the signal strength between central and forward-looking particle detectors, which may be measured by the CMS and TOTEM experiments at LHC.
Meanwhile, the dominant source of uncertainties for EAS simulations is presently related to the model treatments of very high energy pion-air interactions.This can be constrained indirectly by QUARKS-2016 studying other air shower observables, notably, by measurements of the maximal muon production depth.In particular, the results on X μ max from the Pierre Auger experiment strongly disfavor the copious production of baryon-antibaryon pairs in pion-air collisions, predicted by the EPOS model.
It is noteworthy that even for QGSJET-II-04 there is a certain tension between the data of the Pierre Auger experiment on X max and X μ max [20]: the latter point towards a heavier CR composition, compared to the former.According to the analysis presented in this work, one may try to achieve a consistency between the two results by modifying the treatment of pion-air collisions.However, as the potential changes would impact X μ max much stronger than X max , this would imply to aim at higher inelastic cross section and/or softer hadron production spectra for such interactions, which would push one towards a predominantly light, presumably proton-dominated CR composition.
The author acknowledges the support by Deutsche Forschungsgemeinschaft (project OS 481/1).

Figure 4 .
Figure 4. Schematic view of the initial part of the parton cascade in proton.Left: the cascade starts from the same universal parton Fock state; new partons participating in multiple scattering processes emerge from the cascade development, being characterized by ∝ 1/x distributions for the momentum fraction.Right: proton is represented by a superposition of Fock states consisting of different numbers of large x constituent partons; the more abundant multiple scattering the larger Fock states get involved in the process.

DOI: 10 Figure 8 .
Figure 8.Primary energy dependence of X μmax (E μ ≥ 1 GeV) for p-induced vertical EAS, as calculated using the EPOS-LHC, QGSJET-II-04, and QGSJET models (top blue, middle red, and bottom green solid lines respectively), or applying mixed model descriptions, as explained in the text (dashed, dash-dotted, and dotted lines).
Primary energy dependence of X max for p-induced vertical EAS, as calculated using the EPOS-LHC, QGSJET-II-04, and QGSJET models (top blue, middle red, and bottom green solid lines respectively), or applying mixed model descriptions, as explained in the text (dashed, dash-dotted, and dotted lines).