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

TL;DR

This paper reviews how recent LHC data has improved models of hadronic interactions relevant to ultrahigh-energy cosmic rays, but persistent anomalies suggest possible new physics or the need for further model refinement.

Contribution

It provides an updated comparison of Monte Carlo models with LHC data, highlighting remaining discrepancies and potential avenues for understanding cosmic ray composition and new physics.

Findings

Models now better match LHC data after retuning.

Significant muon number discrepancies remain.

Potential new physics or model improvements needed.

Abstract

Measurements of proton and nuclear collisions at the Large Hadron Collider at nucleon-nucleon c.m. energies up to $\sqrt{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 $\sqrt{s_{NN}}\approx 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}\approx 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 significantly 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 effects may help reduce part of the data--model discrepancies, shed light on the UHECR composition approaching the observed $E_{CR}\approx 10^{20}$ eV cutoff, and uncover any potential new physics responsible of the observed anomalies.

Full text

1 Introduction

Ultrahigh-energy cosmic rays (UHECR), with energies $E_{{}_{\rm CR}}\approx 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 Mollerach:2017idb ; Dawson:2017rsp . When reaching earth, they collide with N,​ O nuclei in the upper atmosphere at c.m. energies, $\sqrt{s}=\sqrt{\rm\,m_{{}_{\rm CR}}^{2}+m_{{}_{\rm N,O}}^{2}+2\cdot E_{{}_{\rm CR}}\cdot m_{{}_{\rm CR}}}\approx\sqrt{\rm\,2\cdot 10^{9}\cdot E_{{}_{\rm CR}}(\rm eV)}\approx 14$–450 TeV, up to 30 times larger than those ever reached in any human-made collider Bluemer:2009zf . UHECR produce gigantic cascades of secondary particles in the atmosphere, called extensive air showers (EAS) Kampert:2012mx , measured in dedicated observatories, such as the Pierre Auger Observatory ThePierreAuger:2015rma and Telescope Array (TA) AbuZayyad:2012ru , 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 $\left<X_{\rm max}\right>$ and the RMS width of its fluctuations $\sigma_{X_{\rm max}}$, and the number and total energy of electrons ($N_{\rm e},E_{\rm e}$) and muons ($N_{\rm\mu},E_{\rm\mu}$) on the ground for various shower zenith angle ($\theta$) 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 corsika Heck:1998vt or conex Bergmann:2006yz air transport programs to event generators such as epos-lhc Pierog:2013ria , qgsjet 01 Kalmykov:1993qe , qgsjet-ii-04 Ostapchenko:2010vb , sibyll 2.1 Ahn:2009wx , or sibyll 2.3c Riehn:2017mfm for the hadronic interactions, plus egs4 Nelson:1985ec 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 Gribov:1968fc , 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 dEnterria:2011twh , 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 $\left<X_{\rm max}\right>$ and $\sigma_{X_{\rm max}}$ dependence on $E_{{}_{\rm CR}}$ indicates a change of cosmic-ray composition from proton-dominated to a heavier mix above $E_{{}_{\rm CR}}\approx 10^{18.5}$ eV Aab:2014kda ; Abbasi:2014sfa , but quantitative differences among the $\left<X_{\rm max}\right>$ and $\sigma_{X_{\rm max}}$ predictions exist that are not fully understood, even though independently each MC generator globally reproduces the LHC data Ostapchenko:2014mna ; Pierog:2015epa ; Ostapchenko:2016wtv . On the other hand, in the same range of $E_{{}_{\rm CR}}$ energies, many measurements Dembinski18 indicate significantly larger muon yields on the ground in particular at large transverse distances from the shower axis AbuZayyad:1999xa ; Aab:2014pza , as well as lower $\mu$ production depths Aab:2014dua , 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 Engel:2015dxa ) decays, and/or they do not feed enough energy into the EAS hadronic component (such as, e.g., through an increased production of baryons Pierog:2006qv ). More speculative explanations have been suggested based on changes in the physics of the strong interaction at energies beyond those tested at the LHC Allen:2013hfa ; Aloisio:2017ooo , or on the production of electroweak sphalerons leading to the final production of many energetic muons Brooijmans:2016lfv . 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 dEnterria:2018kcz that, for the first time, interfaced conex with the standard p-p collider MC pythia 6 Sjostrand:2006za ; Skands:2010ak to assess the impact of heavy-quark and hard jet production on the EAS muon production (Sec. 3).

2 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 Ulrich:2010rg are: (i) the p-p inelastic cross section $\sigma_{\rm inel}$ (and its Glauber extension to proton-air interactions, $\sigma_{\rm pAir}$), (ii) the density of charged particles at midrapidity per unit rapidity $dN_{\rm ch}/d\eta|_{\eta=0}$, and (iii) the fraction of energy carried at forward rapidities. In addition, the average transverse momentum of hadrons $\rm\left<p_{{}_{\rm T}}\right>$ 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 dEnterria:2011twh , 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 $\sqrt{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 $\sigma_{\rm 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 $\sigma_{\rm 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 $\sigma_{\rm tot}\approx 90$–120 mb at $\sqrt{s}=14$ TeV depending on whether they preferred to reproduce the (moderately inconsistent) E710 or CDF measurements at Tevatron dEnterria:2011gsx . A large number of $\sigma_{\rm tot,el,inel}$ measurements have been carried out in the last years at the LHC. At $\sqrt{s}=13$ TeV, TOTEM measures $\sigma_{\rm tot}=110.6\pm 3.4$ mb, with $\sim$72% inelastic, and $\sim$28% elastic components Antchev:2017dia . The MC generators tuned to the experimental results at 7 and 8 TeV reproduce well all existing data (Fig. 1 left). The value of $\sigma_{\rm inel}$ was mostly overestimated by the pre-LHC MC models, leading to an increased $\sigma_{\rm pAir}$ cross section, and thereby a smaller $\left<X_{\rm max}\right>$ 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 Auger:2012wt , has been confirmed by LHC heavy-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 $\sigma_{\rm inel,pPb}$ data points measured by CMS at $\sqrt{s_{{}_{\rm NN}}}=5.02$ TeV Khachatryan:2015zaa .

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_{\rm sat}^{2}$ of a few GeV2. At LHC energies, $Q_{\rm sat}\approx 2$ GeV in p-p collisions, and MC models predict that $\sim$70% of the final hadrons at midrapidity issue from the fragmentation of multiply scattered gluons sharing a low fraction of the colliding hadron momenta, $x=\rm p_{\rm parton}/p_{\rm had}\approx p_{\rm T,had}/\sqrt{s_{{}_{\rm NN}}}\lesssim 10^{-4}$. At such low-$x$ and moderate-$Q^{2}$ values, the parton densities are subject to non-linear QCD (gluon saturation) effects, in particular in the dense nuclear case Albacete:2014fwa . Particle multiplicities are thus particularly sensitive to the modeling of the gluon density at low $x$ and of multiparton interactions (MPI). The models tuned to LHC Run-1 data overall bracket the latest p-p and p-A distributions measured (Fig. 2 left), although improvements are still needed in the tails of the multiplicity distributions either at low (diffraction dominated) and/or at high (large number of MPI) values. The evolution of the charged hadron pseudorapidity density at $\eta=0$ as a function of CR energy is shown in Fig. 2 (right). All models retuned to LHC Run-1 data reproduce well the value $dN_{\rm ch}/d\eta|_{\eta=0}\approx 5.5$ measured in p-p collisions at 13 TeV. The different MC predictions, however, start to rapidly diverge above LHC energies and, at the highest CR energies of 1020 eV, vary within $dN_{\rm ch}/d\eta|_{\eta=0}\approx 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_{\rm lead}/E_{{}_{\rm CR}}$, has an important impact on the EAS development. The model predictions for the CR-energy dependence of the inelasticity are shown in Fig. 3 (left): epos-lhc, qgsjet-ii-04, and qgsjet 01 feature a $\sim$25% increase between 1014 eV and 1020 eV, at variance with the almost flat behaviour observed for pythia 6. Since $\left<X_{\rm max}\right>$ and its average fluctuation $\sigma_{X_{\rm max}}$ are mostly driven by the values of $\sigma_{\rm 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 Pierog:2006qv ) production because of its key role feeding the muonic component Ostapchenko:2016ytp . The LHCf experiment has recently measured neutrons at very forward rapidities ($\eta>10.76$) in p-p collisions at 13 TeV Adriani:2018ess 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, $\rm\left<p_{{}_{\rm T}}\right>$ $\approx Q_{\rm sat}$, where multiple minijets (of a few GeV) are produced that subsequently fragment into lower-$p_{{}_{\rm 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_{{}_{\rm T}}$ cutoff separating hard from soft scatterings that mimics the (slow) power-law evolution of $Q_{\rm sat}$ Albacete:2014fwa , whereas sibyll uses a logarithmically-running $Q_{\rm sat}(s)$ scale, and epos and qgsjet have a fixed $p_{{}_{\rm T}}$ cutoff 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_{{}_{\rm T}}$ in Fig. 4 (left): pythia 6 shows a faster $\rm\left<p_{{}_{\rm T}}\right>$ increase than the rest of models, approaching $\rm\left<p_{{}_{\rm T}}\right>\approx 1$ GeV at the highest energies. All RFT models predict a small rise of $\rm\left<p_{{}_{\rm T}}\right>$ plateauing at $\sim$0.6 GeV, except epos-lhc that includes final-state collective parton flow and thereby $\rm\left<p_{{}_{\rm T}}\right>$ increases with multiplicity (and thus $\sqrt{s}$) as seen in Fig. 4 (right). Whereas epos-lhc and pythia 6 reproduce the $N_{\rm ch}$- dependence of $\rm\left<p_{{}_{\rm T}}\right>$ (through different final-state mechanisms: collective flow and colour reconnection), qgsjet-ii-04 clearly misses the data. If one included more pQCD activity in the RFT models, at the expense of smaller final-state effects in the case of epos-lhc, and/or a nuclear-size dependent $Q_{\rm sat}$ value in the initial state (as expected on general grounds) Albacete:2014fwa , $\rm\left<p_{{}_{\rm T}}\right>$ would increase and approach the pythia p-p predictions.

3 Muon anomalies: Data vs. MC simulations

Starting at $E_{{}_{\rm CR}}\approx 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 $\sim$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 $\mu$ shower components) Dembinski18 . In dEnterria:2018kcz , 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 ($\lesssim$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 $\mu$ sources in pythia 6 are the decays of light-quark mesons (charged $\pi,\,k$) from minijet fragmentation, and that heavy-quarks 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 $\mu^{\pm}$ 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 Pierog:2017nes .

4 Summary

The determination of the identity (mass) of the highest-energy 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.

Acknowledgments – Common work and multiple discussions with Tanguy Pierog are gratefully acknowledged.