Hadronic interactions and cosmic rays at ultra high energies

Observations of ultra high energy cosmic rays allow in principle to obtain information about the properties of hadronic interactions in an energy range not accessible with particle accelerators. This task is however complicated by the fact that the composition of the primary CR is not known and must be estimated from the same data set. Solving the ambiguity between composition and hadronic interaction modeling is a central problem for UHECR observations. Recently the Pierre Auger collaboration has presented an estimate of the p–air interaction length at a laboratory energy E


INTRODUCTION
The study of the longitudinal development of the showers generated by ultra high energy (UHE: E 0 10 18 eV) cosmic rays in the Earth's atmosphere allows in principle to measure the chemical composition of these particles, and at the same time to obtain information about the properties of hadronic interaction in an energy range not accessible to accelerators.These two goals are unfortunately in conflict with each other.Uncertainties in the modeling of the properties of hadronic interactions at high energy limit the precision that can be achieved in the measurement of the composition of the primary radiation; viceversa, uncertainties in the composition of the primary particles do not allow an easy determination of the properties of the showers generated by protons (or other nuclear species).The best strategy to solve this ambiguity is the object of an active field of investigations.
At the present moment there is still a significant uncertainty both on the composition of cosmic rays in the UHE range, and on the quality of the description of their showers by the existing montecarlo codes.The analysis of the current situation is difficult also because the observations made by the Pierre Auger [1,2], HiRes [3] and Telescope Array [4] collaborations are not in agreement.The HiRes and Telescope Array results appear to be consistent with a composition that is proton dominated in the entire energy range, while the data of the Pierre Auger Observatory seem to imply a composition that becomes gradually enriched in heavy nuclei with increasing energy.This disagreement is a serious obstacle in the progress of the field, since a consistent interpretation of all data is currently impossible.It is clear that a satisfactory resolution of these discrepancies is urgent and important.
Recently the Pierre Auger collaboration [5] has presented a measurement of the interaction length for protons with energy E 0 10 18.24 eV (that corresponds to a c.m. energy √ s 57 TeV for nucleonnucleon interactions, that corresponds to a p-air cross section 505 ± 23(stat) + 28 − 36(sys) mb (this must be understood as the "production" cross section, that refers to interactions where at least one additional pion is produced, excluding the elastic and target fragmentation processes).Converting the a e-mail: paolo.lipari@roma1.infn.itThis is an Open Access article distributed under the terms of the Creative Commons Attribution License 2.0, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.p-air cross section into a a pp one using Glauber theory [8], the Auger collaboration arrives to the estimate: pp inel = 90 ± 7(stat) + 9 − 11(sys) ± 1.5(Glauber) mb, where the last error corresponds to the uncertainty in the conversion.It should be noted that the measurement of the p-air cross section was performed without a precise estimate of the fraction of protons at the energy considered, and that the systematic error (estimated at 6-7%) claims to include the effect of uncertainties in the modeling of the properties of particle production in high energy hadronic interactions.

EPJ Web of Conferences
The fundamental idea behind this measurement has been introduced by the Fly'e Eye collaboration in 1984 with the publication [7] of an estimate of the p-air cross section at √ s = 30 TeV ( prod pAir = 530 ± 66 mb).The method relies on the observation of the longitudinal development of the cosmic ray showers using the fluorescence light technique.In principle the direct observation of the distribution of X 0 , the position of first interaction of the primary particle could directly measure the p-air interaction length, however X 0 is not observable.On the other hand, for each individual shower it is possible to obtain an estimate of the depth of maximum development X max .The Fly's Eye group suggested that the p-air interaction length can be inferred from the shape of the tail of the X max distribution.Two arguments suggests that this estimate does not require a precise knowledge of the CR composition, and is also only weakly dependent on the properties of hadronic interactions.
1. Protons are the most penetrating component of the primary particles, and therefore they will become the dominant source of the observed showers for a sufficiently high X max (unless of course they constitute only a very small fraction of the CR). 2. Large positive fluctuations in the X max distributions are (likely) controled by large fluctuations in the position of the first interaction point.Accordingly, the tail of the X max distribution should asymptotically approach an exponential form with a slope equal to the p-air interaction length, These points are illustrated in Fig. 1 that shows the X max distributions for showers generated by protons and nuclei (Helium, Oxygen and Iron) of the same energy (E 0 = 10 18.25 eV).The most penetrating showers are generated by protons that, for large high X max have an exponential distribution.For a more quantitative analysis of the relation between the shape of the X max distribution for proton induced showers and the interaction length p (E 0 ) we note that X max , the maximum of the longitudinal development of a shower, can be decomposed as the sum of 2 contributions: X max = X 0 + Y , where X 0 is the position of first interaction, and Y the maximum of the shower development measured from this origin.At very high energy (when particle decays are suppressed) it is a good approximation to consider the distributions of X 0 and Y as independent.The distribution dN/dX max , that we will denote as F (X max ), can then be expressed as the simple convolution of two independent functions.Using the 07004-p.2UHECR 2012 notation dN/dY ≡ G(Y ) and the well known exponential form of the X 0 distribution one obtains: where p is the p-air interaction length (and we have left implicit the dependence on the energy E 0 ).For large X max the distribution is well approximated by an exponential of slope p , if the term in square parenthesis in equation ( 1) converges sufficiently rapidly to a constant when the upper limit of the integration becomes large.This is the case if the distribution G(Y ) vanishes more rapidly, than e −Y/ p when Y → ∞.It can be shown (see below) that the high Y behaviour of G(Y ) is in fact also nearly exactly proportional to e −Y/ p , so that the function F (X max ) does indeed become asymptotically (for large X max ) a simple exponential of slope p , but the convergence is only slow.Differentiating and reorganizing the result, one can recast equation ( 1) in the form: This relates the observable quantity and the interaction length p .The two quantities are equal when (and if) the ratio G(X)/F (X) is neglibly small.More in general to infer the value of p from the observed shape of the distribution one needs to apply a correction that requires good theoretical control of the shape of the function G(Y ) that describes the development of the shower from the point of first interaction.
In [7] the Fly's Eye group fitted the tail of the distribution as a simple exponential and assumed that the relation between the slope and p is approximately linear ( = K p ), introducing a (model dependent) "K-factor".A simple linear relation between is in fact not valid, as can be understood inspecting equation (2) and noting that the shape of the distribution G(Y ) also depends on the hadron interaction lengths (for energies E ≤ E 0 ).Changing the p-air interaction length at energy E 0 , implies also a change also at lower energies, and a modification of G(Y ).One possible form for the (modified) energy dependence of the interaction length is: The single parameter , and accordingly p (E 0 ), can be obtained fitting the shape of the data.Note that for a complete description of the shower development one also needs to define the meson interaction lengths.Perhaps the simplest (but by no mean unique) possibility (as adopted in [5])) is to rescale the nucleon and meson interaction lengths with the same factor.The Auger collaboration [5] has fitted the shape of the X max distribution of showers with energy around E 0 10 18.24 eV, in the interval X max ∈ [768, 1004] g cm −2 as an exponential, obtaining a slope = 55.8 ± (stat ± 1.6 (sys) g cm −2 .Converting this quantity into a cross section using the expression: = / m (with m the average mass of a target nucleus in air) one obtains a cross section 429 mbarn.This quantity is not the physical cross section because one needs to apply a correction to take into account that in the region considered one has not yet reached the expected asymptotic behaviour p .This correction was estimated by the Auger collaboration using Montecarlo methods.Using the QGSjet01 [10], QGSJetII [11], SIBYLL [9] and EPOS [12] codes resulted in estimates for the p-air production cross section of 523.7, 502.9, 496.7 and 497.7 mbarn respectively (this corresponds to a correction factor of order 1.16-1.22).

MONTECARLO CALCULATIONS
A crucial problem is if the method of using the spread among the available model to estimate the systematic uncertainty associated with the modeling of hadronic interactions is adequate.This procedure is certainly reasonable but lacks a really robust justification, and could potentially lead to an underestimate of the systematic error considering that in fact the models used for comparison have many common assumptions.To explore this issue we have we have performed a montecarlo of shower development study using the Sibyll code, and two other models that have been constructed modifiying the inclusive spectra of final state particles.In model A the pion spectra are significantly harder than in Sibyll and consequently the showers are more penetrating.In model B the spectra are significantly softer, and consequently the showers develop faster with a smaller X max (detailed of the codes will be given elsewhere).In all 3 calculations we have used the same interactions lengths.For the hadronnucleon cross sections we have used the parametrizations suggested by the Particle Data Group [6], converting them to hadron-nucleus cross sections using Glauber theory [8].For E 0 = 10 18.24 eV this gives: p (E 0 ) = 45.6 g cm −2 .For each individual shower we have calculated the longitudinal profile, and estimated the X max position fitting the region around the maximum with a parabola.The resulting X max distributions are shown in Fig. 2, and compared with the Auger data at the same energy.For the Sibyll montecarlo one has X max = 748 g cm −2 and a dispersion X max = 62 g cm −2 ; for model A X max = 779 g cm −2 and X max = 63 g cm −2 ; for model B X max = 705 g cm −2 and X max = 67 g cm −2 .
A more detailed view of the Sibyll calculation is shown in 3, that shows the distributions for the quantities X max and Y = X max − X 0 .Inspecting the shape the X max distribution one can see that the tail does in fact take aymptotically an exponential form with a slope approximately equal to p (E 0 ).The Y distribution also appears to have an exponential tail with similar slope.In fact the montecarlo distribution shown in Fig. 3 is well fitted (for Y ≥ 930 g cm −2 ) with an exponential of slope p (E 0 ).The exponential behaviour at high Y of the distribution can be understood observing that the showers with the longest Y development are mostly events where the first interaction is a diffractive scattering that result in a final state that contains one nucleon that carries away most of the primary energy (and therefore remains with an energy approximately equal to E 0 ).The analytic expression used to fit the montecarlo results for the Y distribution was then convoluted (using Eq. ( 1), with the X 0 distribution to obtain the form of F (X max ).The result (also shown in Fig. 3) is a very good description of the X max distribution, showing the consistency of the procedure.The same procedure was performed for the MC calculations that uses distorted spectra for final state particles (model A and model B).The relation between the observable (X) and the physical interaction length p (E 0 ) is shown in Fig. 4 where (X max ) is plotted as a function 07004-p.4 UHECR 2012 Figure 3. 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 eV using the Sibyll code and a shower montecarlo model.The thin (black) line is a fit to the Y distribution (as 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 p (E 0 ), and provides an excellent fit to the MC results. of X max .In all cases the slope (X) decreases gradually approaching the asymptotic value p (E 0 ).In the range X max 800-1000 g/cm 2 one needs to apply a correction of a factor 1.2-1.4.This result (based on large and arbitrary spectral distortions) suggests that the estimate of the Auger collaboration based on the available models is perhaps too small, but not by a large factor.At this point one can observe that this line of arguments can be naturally developed in interesting ways, that can simultaneously give information on the UHECR composition and on the properties of hadronic interactions.An illustration of this idea is given in Fig. 5 that again shows the Auger data for the X max distribution at E 0 10 18.25 eV, together with the MC predictions for the 3 models of proton showers discussed above (Sibyll and the hard and soft models A and B).In this case the MC predictions normalized to the data in the region X max > 800 g cm −2 .The 3 models have the same interaction lengths, and therefore approximately the same slope at large X max (always f order 50 g/cm 2 as in the data), however, differ significantly at lower X max .If one assumes that the Sibyll model is correct one remains with little room for adding other nuclear species to the CR flux at the energy considered, because a pure proton composition does give (in first approximation) a good description of the entire distibution.On the other hand, model A (with harder spectra of secondary particles and more penetrating showers) requires the addition of a significant fraction of nuclei to explain the data.The proton fraction can in fact be obtained fitting the high X max tail of the distribution, while a fit of the low X max edge gives information about the heaviest nuclear species that (assuming the validity of the model) is abundant in the flux.Finally a model with very soft spectra (and very rapidly developing showers) such as model B is in fact already excluded from the data, because if one fits the high X max tail of the distribution, the predictions of the rate of showers with small X max is larger than the data, and the addition of nuclei can only worsen the situation.The conclusion is therefore that a range of hadronic models can already be excluded with present data.
A study of this type for showers in different energy bands has the potential to give simultaneously valuable information about (i) the energy dependence of the p-air interaction length; (ii) the evolution with energy of the chemical composition of UHECR; and (iii) some global properties (such as the inclusive spectra) of hadronic interactions in an energy range not accessible with accelerators.
I would like to thanks the organizers for a timely, stimulating and very interesting workshop.

Figure 1 .
Figure 1.Distribution of X max for showers induced by p and nuclei (He, O and Fe) of energy E 0 = 10 18.25 eV.The distributions (normalized to unit area) are calculated with a MC code that uses the Sibyll [9] interaction model.

Figure 2 .
Figure 2. The lines are the X max distributions for showers generated by p and of energy E 0 = 10 18.25 eV (all distributions are normalized to the same area).For protons three different models have been used.The (black) line labeled S uses the Sibyll code.The (red dashed) line labeled A uses hard pion spectra in the final state and has very penetrating showers.The (blue dashed) line labeled B uses soft pion spectra and has fastly developing showers.The 3 models have identical interaction lengths.

Figure 4 .
Figure 4. Ratios (X max )/ p for the 3 MC calculations shown in fig. 2 and fig. 5.The lines are calculated from analytic fits to the montecarlo results that have asymptotic slope p (E 0 ).

Figure 5 .
Figure 5.The points are the measurement of the X max distribution for showers with E 0 10 18.24 eV obtained by the Auger observatory [5].The lines are MC distributions for proton showers of the same energy.The (black) line labeled S uses the Sibyll code; the (red dashed) line labeled A uses hard pion spectra, the (blue dashed) line labeled B uses soft pion spectra.The MC distributions are normalized to the data for the interval X max ∈ [800, 1100] g/cm 2 .