Spectra and composition of ultrahigh-energy cosmic rays and the measurement of the proton-air cross section

The shape of the longitudinal development of the showers generated in the atmosphere by very high-energy cosmic ray particles encodes information about the mass composition of the flux, and about the properties of hadronic interactions that control the shower development. Studies of the shape of the depth of maximum distributions of showers with $E\gtrsim {10}^{17.3}\mathrm{eV}$ measured by the Pierre Auger Observatory, suggest, on the basis of a comparison with current models, that the composition of the cosmic ray flux undergoes a very important evolution, first becoming lighter and then rapidly heavier. These conclusions, if confirmed, would have profound and very surprising implications for our understanding of the high-energy astrophysical sources. Studies of the shape of the depth of maximum distribution in the same energy range have been used by the Auger and Telescope Array Collaboration to measure the interaction length of protons in air, a quantity that allows to estimate the $pp$ cross sections for values of $\sqrt{s}$ well above the LHC range. In this paper we argue that it is desirable to combine in a self-consistent way the studies of the cosmic ray composition with those aimed at the measurement of the $p$-air cross section. The latter necessarily provide information on both the fraction of protons in the flux and the properties of the showers generated by protons, that can be of great help in decoding the cosmic ray composition and in constraining the models of air shower development. At the same time a good determination of the cosmic ray composition, in particular of the helium component, is essential to infer correctly from the data the proton-air cross section.

Figure 1

All particle energy spectrum of very high-energy cosmic rays. The measurements are by Auger [19, 20], Telescope Array [21] and TALE [22]. The lines are fits to the data reported in the original publications.

Figure 2

Measurements of the average $⟨X_{\max }⟩$ (top panel) and dispersion $W=[⟨X_{\max }^2⟩-⟨X_{\max }{⟩}^2{]}^{1/2}$ (bottom panel) of the depth of maximum distributions measured by the Pierre Auger Observatory in different energy intervals [25]. The predictions for protons and iron nuclei particles are calculated using the hadronic models QGSJet II-04 [26], EPOS-LHC [27], and Sibyll 2.3c [28].

Figure 3

Allowed region in the plane $\{⟨X_{\max }⟩,W\}$ calculated using the hadronic model QGSJetII-04 [26] at the energy $E={10}^{17.5}\mathrm{eV}$ and considering the contributions of 5 nuclei ($p$, ${}^4\mathrm{He}$, ${}^{14}\mathrm{N}$, ${}^{28}\mathrm{Si}$, and ${}^{56}\mathrm{Fe}$). The lines show points that can be generated by the combinations of two nuclei.

Figure 4

Allowed regions in the plane $\{⟨X_{\max }⟩,W\}$ calculated using three hadronic models (QGSJetII-04, Epos-LHC, and Sibyll 2.3c) for two energies $E={10}^{17.5}$ and ${10}^{19.5}\mathrm{eV}$.

Figure 5

The Auger measurements of the average and dispersion of the depth of maximum distributions (shown in Fig. 2) are represented as points in the plane of the rescaled variables $x$ and $y$ [see Eqs. (5) and (6)]. Three panels show the results for the three hadronic models QGSJetII-04, Sibyll 2.3c and EPOS-LHC. The contour of the allowed region in the $\{x,y\}$ plane (for each model) has a weak energy dependence, and in each panel the contour is shown for three values of the energy (${\log }_{10}E=17.5$, 18.5 and 19.5), the difference between the contours remains always small. The broken line connects data points in adjacent energy intervals, with the highest energy point the one with the lowest $x$ and $y$ values. The last panel shows the trajectory in the space $\{x,y\}$ for the composition model discussed in the text (and shown in Fig. 8) calculated using the EPOS-LHC model (representative values of $\log [E(\mathrm{ev})]$ are labeled). In all panels the shaded area shows the region of parameters space allowed for the combination of five nuclei considered. The darker shaded areas indicate the parameter regions that are physically possible for a fixed proton fraction [with values $f_p=0.75$, $f_p=0.5$, $f_p=0.25$, and $f_p=0$ as marked, a pure proton composition ($f_p=1$) corresponds to the corner at the upper right].

Figure 6

Interpretation of the Auger measurements $⟨X_{\max }⟩$ and $W$ at energy $E={10}^{18.25}\mathrm{eV}$ in terms of a composition formed by protons and a second component of mass $A$, using the EPOS-LHC an Sibyll 2.3c models. The thick lines show the proton fraction needed to reproduce the measured value of $⟨X_{\max }⟩$ as a function of $A$ (with the shaded area a one sigma uncertainty interval). The ellipses show the (one standard deviation) the allowed region in the plane $\{A,f_p\}$ calculated taking into account the measurement of the width $W$ of the depth of maximum distribution.

Figure 7

Width of the depth of maximum distribution predicted at the energy $E={10}^{18.25}\mathrm{eV}$ for a composition formed by protons and nuclei of mass number $A$. For each value of $A$ the proton fraction $f_p$ is determined by the requirement to reproduce the value of $⟨X_{\max }⟩$ obtained by Auger [25] (the shaded area is a $1\sigma$ uncertainty band). The measured value of $W$ (with a $1\sigma$ error) is shown as the horizontal band. The top (bottom) panel uses the EPOS-LHC (Sibyll 2.3c) model.

Figure 8

Model of the CR energy spectrum and composition constructed to reproduce the Auger data. The spectrum below the ankle (shown as a thick solid line) is formed by two component, one of protons and the other of nitrogen and silicon (with equal abundances) that have both power law form with superexponential cutoffs. The spectrum above the ankle (thick dashed line) is formed by the contributions of five nuclei (protons, helium, nitrogen, silicon, and iron) that have a hard power law spectra (with the same slope), with rigidity dependent cutoffs (see main text for more details). The data points are from Auger [19, 20].

Figure 9

Fraction of protons in the CR flux as a function of energy. The thin lines are estimates of the proton fraction obtained from the measurements of the average depth of maximum obtained by Auger [25] comparing with the predictions of the EPOS-LHC model and assuming that the composition is formed by the combination of protons and one nuclear component (helium, nitrogen, silicon, and iron) [see Eq. (10)]. The thick solid line is the proton fraction for the model discussed in the text (and shown in Fig. 8). The dashed line is the proton fraction in the very high-energy component discussed by the Auger Collaboration [7, 8]. The Auger model does not include the subankle component, but is in good agreement with the superankle component discussed in this paper.

Figure 10

Total and elastic $pp$ cross sections plotted as a function of the c.m. energy $\sqrt{s}$. The points are measurements of the TOTEM detector at the LHC [37, 38, 39]. The solid lines are fits to the total and elastic cross sections that are quadratic in $\log s$ [39]. The shaded areas are estimates of the uncertainties.

Figure 11

Proton interaction length in air plotted as a function of the projectile laboratory energy. The central line and the shaded area are calculated using the best fits to the total and elastic $pp$ cross sections and the uncertainties shown in Fig. 10, and using the algorithms of Glauber and Matthiae [18] to estimate the proton-nucleus cross sections. The points are the estimates of the proton-air interaction length obtained from measurements of the longitudinal developments of cosmic ray showers by Fly’s Eye [13], Auger [14], and Telescope Array [15, 16]. The lowest energy point is calculated from the measurements at $\sqrt{s}=13\mathrm{TeV}$ by TOTEM at LHC [39].

Figure 12

Distributions of $X_{\max }$ for showers generated by particles with energy $E={10}^{18.25}\mathrm{eV}$. The distributions are calculated with Monte Carlo methods using the Sibyll 2.1 model and the $p$-air interaction length shown in Fig. 11, for four different nuclei: protons, ${}^4\mathrm{He}$, ${}^{16}\mathrm{O}$, and ${}^{56}\mathrm{Fe}$. The high $X_{\max }$ part of the distributions has been fitted with a simple exponential form: $dN/d{X}_{\max }\propto e^{-X_{\max }/\mathrm{\Lambda }}$.

Figure 13

Parameters of the distributions of $X_{\max }$ for the showers four different nuclei (protons, helium, oxygen, and iron) at energy $E={10}^{18.25}\mathrm{eV}$ shown in Fig. 12. The top panel shows the average $⟨X_{\max }⟩$, and the line is a linear fit for the relation $⟨X_{\max }⟩=X_0+D\log A$. The middle panels shows the r.m.s. of the distributions $W=\sqrt{⟨X_{\max }^2⟩–⟨X_{\max }{⟩}^2}$, and the two lines are analytical approximations of the $A$ dependence ($W_A=W_p{A}^{-0.25}$ and $W_A=W_p[1-a\log A+b(\log A{)}^2{]}^{1/2}$). The bottom panel shows the parameter $\mathrm{\Lambda }$ that fits the high-energy part of the distributions. The line is a polynomial fit to the $\log A$ dependence.

Figure 14

Distribution of $X_{\max }$ (bigger points) and $Y=X_{\max }-X_0$ (smaller points) calculated for proton showers with energy $E_0={10}^{18.25}\mathrm{eV}$ using the Sibyll code and a shower Monte Carlo model. The thin (black) line is a fit to the $Y$ distribution described in the main text. The thick (red) line is obtained convoluting the previous result with an exponential with slope equal to the interaction length ${\lambda }_p(E_0)$.

Figure 15

Slope $\mathrm{\Lambda }(X_{\max })$ of the $X_{\max }$ distribution obtained with a Monte Carlo calculation for protons of energy $E_0={10}^{18.25}\mathrm{eV}$ and shown in Fig. 14. The dashed line shows the $p$-air interaction length used in the Monte Carlo calculation. The shaded area shows the range of $X_{\max }$ and the best fit value for $\mathrm{\Lambda }$ in the study of Auger in [14].

Figure 16

The points show the depth of maximum distribution observed by Auger for showers with a reconstructed average energy around ${10}^{18.25}\mathrm{eV}$ [14] and [40, 41]. The thick solid line is the prediction of the Sibyll 2.1a (that is the model in [36] with the interaction lengths discussed in the main text) for a pure proton composition. The thick dashed line is the same distribution with the showers deeper by $30\mathrm{g}/{\mathrm{cm}}^2$ (the difference in average depth of maximum between showers simulated with the Sibyll 2.1 and Sibyll 2.3c models). In this case the proton fraction is of order $f_p\simeq 0.55$.

Figure 17

Depth of maximum distributions of showers detected by Auger in different energy intervals [40, 41]. The lines are fits discussed in the main text.

Figure 18

Fits to the depth of maximum distributions of the CR showers detected by Auger [40, 41] in six different energy intervals plotted together for comparison (all curves are normalized to unity). Comparisons of the fits with the data are shown in Fig. 17, the labeling of the curves ($\mathrm{a},\dots ,\mathrm{f}$) is the same as in the panels in that figure.

Figure 19

Best fits parameters for the depth of maximum distributions of the Auger data (see Fig. 17). The top panel shows $X_{\mathrm{peak}}$, the value of $X$ where the distribution has its maximum value. The middle panel shows ${\sigma }_{X,\mathrm{left}}$, the width of the Gaussian that fits the distribution for $X<X_{\mathrm{peak}}$. The bottom panel shows $\mathrm{\Lambda }$, the slope of the tail of the distribution for large $X$ values.