Ultrahigh-energy cosmic rays: Anomalies, QCD, and LHC data

. Measurements of proton and nuclear collisions at the Large Hadron Collider at nucleon-nucleon c.m. energies up to √ s NN = 13 TeV have improved our understanding of hadronic interactions at the highest energies reached in collisions of cosmic rays with nuclei in the earth atmosphere, up to √ s NN ≈ 450 TeV. The Monte Carlo event generators ( epos , qgsjet , and sibyll ) commonly used to describe the air showers generated by ultrahigh-energy cosmic rays (UHECR, with E CR ≈ 10 17 –10 20 eV) feature now, after parameter retuning based on LHC Run-I data, more consistent predictions on the nature of the cosmic rays at the tail of the measured spectrum. However, anomalies persist in the data that cannot be accommodated by the models. Among others, the total number of muons (as well as their maximum production depth) remains signiﬁcantly underestimated (overestimated) by all models. Comparisons of epos , qgsjet , and sibyll predictions to the latest LHC data, and to collider MC generators such as pythia , indicate that improved description of hard multiple minijet production and nuclear e ﬀ ects may help reduce part of the data–model discrepancies, shed light on the UHECR composition approaching the observed E CR ≈ 10 20 eV cuto ﬀ , and uncover any potential new physics responsible for the observed anomalies.


Introduction
Ultrahigh-energy cosmic rays (UHECR), with energies E CR ≈ 10 17 -10 20 eV, are produced and accelerated in (poorly-known) extreme astrophysical environments.Their exact extragalactic sources and their nature, protons or heavier ions, remain still open questions today [1,2].When reaching earth, they collide with N,O nuclei in the upper atmosphere at c.m. energies, • E CR (eV) ≈ 14-450 TeV, up to 30 times larger than those ever reached in any human-made collider [3].UHECR produce gigantic cascades of secondary particles in the atmosphere, called extensive air showers (EAS) [4], measured in dedicated observatories, such as the Pierre Auger Observatory [5] and Telescope Array (TA) [6], that combine the information from (i) the lateral distributions of secondary particles in surface detectors, and (ii) the fluorescence photons produced by nitrogen molecules excited along the shower.Key EAS observables are the average depth of the shower maximum X max and the RMS width of its fluctuations σ X max , and the number and total energy of electrons (N e , E e ) and muons (N µ , E µ ) on the ground for various shower zenith angle (θ) inclinations.By comparing these shower properties to the predictions of Monte Carlo (MC) models, the energy and mass of the incoming particles can be inferred.
The comparison of the measured EAS properties to the MC predictions is commonly done by interfacing the cor- * e-mail: david.d'enterria@cern.chsika [7] or conex [8] air transport programs to event generators such as epos-lhc [9], qgsjet 01 [10], qgsjet-ii-04 [11], sibyll 2.1 [12], or sibyll 2.3c [13] for the hadronic interactions, plus egs4 [14] for the electromagnetic (e.m.) cascade evolution.The hadronic interaction models describe particle production in high-energy proton and nuclear collisions based on the Gribov's Reggeon Field Theory (RFT) framework [15], extended to take into account perturbative quantum chromodynamics (pQCD) scatterings in (multiple) hard parton collisions via "cut (hard) Pomerons" (identified diagrammatically as a ladder of gluons).With parameters tuned to reproduce the existing accelerator data [16], all those MC generators are able to describe the overall EAS properties, although anomalies persist in the UHECR results that cannot be accommodated by the models.On the one hand, the X max and σ X max dependence on E CR indicates a change of cosmicray composition from proton-dominated to a heavier mix above E CR ≈ 10 18.5 eV [17,18], but quantitative differences among the X max and σ X max predictions exist that are not fully understood, even though independently each MC generator globally reproduces the LHC data [19][20][21].On the other hand, in the same range of E CR energies, many measurements [22] indicate significantly larger muon yields on the ground in particular at large transverse distances from the shower axis [23,24], as well as lower µ production depths [25], than predicted by all models.
The difficulties in reconciling aspects of the UHECR data with the MC generator predictions suggest that the best models of hadronic interactions are missing some  physics ingredient.For the muon component, either they do not account for processes that produce harder muons, such as from, e.g., jets or heavy-quark (in particular charm [26]) decays, and/or they do not feed enough energy into the EAS hadronic component (such as, e.g., through an increased production of baryons [27]).More speculative explanations have been suggested based on changes in the physics of the strong interaction at energies beyond those tested at the LHC [28,29], or on the production of electroweak sphalerons leading to the final production of many energetic muons [30].This writeup compares the RFT Monte Carlo predictions to a basic set of QCD observables measured recently in proton and nuclear collisions at the LHC (Sec.2), and summarizes the results of a recent study [31] that, for the first time, interfaced conex with the standard p-p collider MC pythia 6 [32,33] to assess the impact of heavy-quark and hard jet production on the EAS muon production (Sec.3).

LHC data versus UHECR Monte Carlo generators
The depth of the maximum of the EAS depends mainly on the characteristics of multiparticle production in the first few generations of hadronic interactions in the atmosphere of the incoming UHECR.Among the key observables in hadronic collisions with biggest impact on the longitudinal development of EAS [34] are: (i) the pp inelastic cross section σ inel (and its Glauber extension to proton-air interactions, σ pAir ), (ii) the density of charged particles at midrapidity per unit rapidity dN ch /dη| η=0 , and (iii) the fraction of energy carried at forward rapidities.
In addition, the average transverse momentum of hadrons p T is a sensitive probe of the underlying modeling of pQCD processes.The level of agreement between the MC predictions for such observables and the first LHC data (mostly p-p collisions at 7 and 8 TeV) was carefully reviewed in [16], where it was found that the pre-LHC models overall bracketed the experimental measurements and only small modifications, such as those that led to the new qgsjet-ii-04, epos-lhc and sibyll 2.3c releases, were required.A few representative recent experimental results from LHC, in particular, pp at √ s = 13 TeV, not included in the current Run-1 MC parameter tuning, are compared to the MC predictions below.

Inelastic p-p cross section
A fundamental quantity of all MC models and key in cosmic-ray physics, as it chiefly determines the UHECR penetration in the atmosphere, is the total p-p cross section σ tot and its separation into elastic and inelastic (and, in particular, diffractive) components, as well as its extension to proton-air interactions.Although non-computable from the QCD Lagrangian, the values of σ tot,el,inel are constrained by fundamental quantum mechanics relations such as the Froissart bound, the optical theorem, and the dispersion relations.The pre-LHC models predicted σ tot ≈ 90-120 mb at √ s = 14 TeV depending on whether they preferred to reproduce the (moderately inconsistent) E710 or CDF measurements at Tevatron [35].A large number of σ tot,el,inel measurements have been carried out in the last years at the LHC.At √ s = 13 TeV, TOTEM measures σ tot = 110.6 ± 3.4 mb, with ∼72% inelastic, and ∼28% elastic components [36].The MC generators tuned to the experimental results at 7 and 8 TeV reproduce well all existing data (Fig. 1 left).The value of σ inel was mostly overestimated by the pre-LHC MC models, leading to an increased σ pAir cross section, and thereby a smaller X max penetration than predicted now by the retuned models.In addition, the validity of the geometric Glauber multiple scattering model, used to extrapolate from p-p to p-nucleus cross sections [37], has been confirmed by LHC heavy- TeV measured by ATLAS compared to the predictions of various MC generators [41].Right: Evolution of the charged-particle pseudorapidity density at midrapidity, dN ch /dη| η=0 , as a function of CR energy: Data points (from the compilations [16,39]) are compared to the latest MC predictions [31].ion measurements.Figure 1 (right) shows the inelastic p-Pb cross section as a function of collision energy with the Glauber prediction (solid black) reproducing the σ inel,pPb data points measured by CMS at √ s NN = 5.02 TeV [38].

Central particle multiplicity
Inclusive hadron production in high-energy proton and nuclear collisions peaks at central rapidities and receives contributions from "soft" and "hard" interactions between the partonic constituents of the colliding hadrons.Soft (respectively hard) processes involve mainly t-channel partons of virtualities Q 2 typically below (resp.above) a scale Q 2 sat of a few GeV 2 .At LHC energies, Q sat ≈ 2 GeV in pp collisions, and MC models predict that ∼70% of the final hadrons at midrapidity issue from the fragmentation of multiply scattered gluons sharing a low fraction of the col-liding hadron momenta, x = p parton /p had ≈ p T,had / √ s NN  energies and, at the highest CR energies of 10 20 eV, vary within dN ch /dη| η=0 ≈ 10-25.

Forward particle production
The fraction of the primary CR energy transferred to secondary particles after removing the most energetic "leading" hadron emitted at very forward rapidities, called inelasticity K = 1 − E lead /E CR , has an important impact on the EAS development.The model predictions for the CRenergy dependence of the inelasticity are shown in Fig. 3 (left): epos-lhc, qgsjet-ii-04, and qgsjet 01 feature a ∼25% increase between 10 14 eV and 10 20 eV, at variance with the almost flat behaviour observed for pythia 6.Since X max and its average fluctuation σ X max are mostly driven by the values of σ inel and inelasticity, proton-EAS simulated with RFT models feature smaller penetration in the atmosphere than those predicted by pythia 6.Of particular interest is the forward particle (in particular baryon [27]) production because of its key role feeding the muonic component [43].The LHCf experiment has recently measured neutrons at very forward rapidities (η > 10.76) in p-p collisions at 13 TeV [42] finding that none of the current MC models can reproduce their total nor differential yields (Fig. 3 right).

Transverse momentum spectra
In hadronic collisions, the partonic cross sections peak around the saturation scale, p T ≈ Q sat , where multiple minijets (of a few GeV) are produced that subsequently fragment into lower-p T hadrons.This regime dominated by gluon saturation is modeled in various effective ways in the MC generators.The pythia MC has an energy-dependent p T cutoff separating hard from soft scatterings that mimics the (slow) power-law evolution of Q sat [40], whereas sibyll uses a logarithmically-running Q sat (s) scale, and epos and qgsjet have a fixed p T cut-off and low-x saturation is implemented through corrections to the multi-Pomeron dynamics.These different behaviours are seen in the energy evolution of the predicted average p T in Fig. 4 (left): pythia 6 shows a faster p T increase than the rest of models, approaching p T ≈ 1 GeV at the highest energies.All RFT models predict a small rise of p T plateauing at ∼0.6 GeV, except epos-lhc that includes final-state collective parton flow and thereby p T increases with multiplicity (and thus √ s) as seen in Fig. 4 (right).Whereas epos-lhc and pythia 6 reproduce the N chdependence of p T (through different final-state mechanisms: collective flow and colour reconnection), qgsjetii-04 clearly misses the data.If one included more pQCD activity in the RFT models, at the expense of smaller finalstate effects in the case of epos-lhc, and/or a nuclear-size dependent Q sat value in the initial state (as expected on general grounds) [40], p T would increase and approach the pythia p-p predictions.

Muon anomalies: Data vs. MC simulations
Starting at E CR ≈ 10 17 eV and with increasing primary CR energy, various measured characteristics of muon production (accessible in inclined EAS) are difficult to reconcile with the MC models: (i) the observed number of muons is underestimated by ∼30% (and even larger values far away from the shower axis), and (ii) the measured muon penetration depth is overestimated (indicating inconsistent primary composition as derived from the total and µ shower components) [22].In [31], we studied the production of decay muons from hard QCD jets and heavy-quarks, absent in the RFT models, using pythia 6 to generate proton-EAS in a hydrogen atmosphere with the same density as air.In general, the pythia 6 results for the muon densities and energies at sea level are found in between those predicted by the RFT MCs, although switching-off heavy-  quark production leaves more room for charged pion and kaon production, leading to a 15% muon increase on the ground (Fig. 5 left).Also, whereas pythia 6 showers feature a bit less muons close to the core ( 200 m), their yield is 10-30% larger at 600-1000 m from the shower axis (Fig. 5 right).These results indicate that the main µ sources in pythia 6 are the decays of light-quark mesons (charged π, k) from minijet fragmentation, and that heavyquarks account for a negligible fraction of the inclusive muons.Although part of the muon excess observed at large radii can be solved adding harder minijet activity in the RFT models, for real UHECR collisions on air, and at variance with the results found in the proton-hydrogen setup, epos-lhc produces more µ ± than qgsjet-ii-04.Thus, additional nuclear (and/or other) effects are still needed to explain the final production of muons observed in the data [44].

Summary
The determination of the identity (mass) of the highestenergy cosmic rays reaching earth relies heavily on MC hadronic generators with parameters tuned to collider data.Representative measurements of basic QCD observables (inelastic cross sections, hadron multiplicity, forward particle production, and mean hadron transverse momenta) measured recently in proton and nucleus collisions at the LHC have been compared to the predictions of cosmic-ray MC models.The agreement is overall good, but improvements are needed in the modeling of very forward particle production and of semi-hard multiparton interactions.Solving the current muon anomalies observed in the UHECR data requires improved multiple minijet (but, seemingly, not heavy-quarks, at least for proton-induced showers) production plus a better overall control of nuclear effects before any potential new-physics contribution can be isolated.

Figure 2 .
Figure2.Left: Per-event charged-particle multiplicity probability in p-p at √ s = 13 TeV measured by ATLAS compared to the predictions of various MC generators[41].Right: Evolution of the charged-particle pseudorapidity density at midrapidity, dN ch /dη| η=0 , as a function of CR energy: Data points (from the compilations[16,39]) are compared to the latest MC predictions[31].

10 − 4 .Figure 3 .
Figure 3. Left: MC predictions for the p-p inelasticity as a function of CR energy [31].Right: Differential neutron cross section measured by the LHCf experiment at η > 10.76 in p-p at √ s = 13 TeV compared to MC predictions [42].

Figure 4 .
Figure 4. Left: Data-model evolution of p T at midrapidity as a function of primary CR energy [31].Right: Data-model comparison of the mean transverse momentum vs. charged-particle multiplicity in p-p collisions at √ s = 13 TeV [41].

Figure 5 .
Figure 5. Left: Number of muons (normalized by E 0.9 CR ) at shower maximum for inclined proton-induced showers as a function of CR energy predicted by different MC generators [31].Right: Number of muons at sea-level as a function of radial distance to the shower axis for proton-induced EAS of E CR = 10 19 eV predicted by different MC generators over the qgsjet-ii-04 value [31].