Interpretation of the Depths of Maximum of Extensive Air Showers Measured by the Pierre Auger Observatory
Abstract
To interpret the mean depth of cosmic ray air shower maximum and its dispersion, we parametrize those two observables as functions of the first two moments of the distribution. We examine the goodness of this simple method through simulations of test mass distributions. The application of the parameterization to Pierre Auger Observatory data allows one to study the energy dependence of the mean and of its variance under the assumption of selected hadronic interaction models. We discuss possible implications of these dependences in term of interaction models and astrophysical cosmic ray sources.
Prepared for submission to JCAP
Interpretation of the Depths of Maximum of Extensive Air Showers Measured by the Pierre Auger Observatory
The Pierre Auger Collaboration

Centro Atómico Bariloche and Instituto Balseiro (CNEAUNCuyoCONICET), San Carlos de Bariloche, Argentina
Centro de Investigaciones en Láseres y Aplicaciones, CITEDEF and CONICET, Argentina
Departamento de Física, FCEyN, Universidad de Buenos Aires y CONICET, Argentina
IFLP, Universidad Nacional de La Plata and CONICET, La Plata, Argentina
Instituto de Astronomía y Física del Espacio (CONICETUBA), Buenos Aires, Argentina
Instituto de Física de Rosario (IFIR)  CONICET/U.N.R. and Facultad de Ciencias Bioquímicas y Farmacéuticas U.N.R., Rosario, Argentina
Instituto de Tecnologías en Detección y Astropartículas (CNEA, CONICET, UNSAM), Buenos Aires, Argentina
National Technological University, Faculty Mendoza (CONICET/CNEA), Mendoza, Argentina
Observatorio Pierre Auger, Malargüe, Argentina
Observatorio Pierre Auger and Comisión Nacional de Energía Atómica, Malargüe, Argentina
Universidad Tecnológica Nacional  Facultad Regional Buenos Aires, Buenos Aires, Argentina
University of Adelaide, Adelaide, S.A., Australia
Centro Brasileiro de Pesquisas Fisicas, Rio de Janeiro, RJ, Brazil
Universidade de São Paulo, Instituto de Física, São Carlos, SP, Brazil
Universidade de São Paulo, Instituto de Física, São Paulo, SP, Brazil
Universidade Estadual de Campinas, IFGW, Campinas, SP, Brazil
Universidade Estadual de Feira de Santana, Brazil
Universidade Federal da Bahia, Salvador, BA, Brazil
Universidade Federal do ABC, Santo André, SP, Brazil
Universidade Federal do Rio de Janeiro, Instituto de Física, Rio de Janeiro, RJ, Brazil
Universidade Federal Fluminense, EEIMVR, Volta Redonda, RJ, Brazil
Rudjer Bošković Institute, 10000 Zagreb, Croatia
Charles University, Faculty of Mathematics and Physics, Institute of Particle and Nuclear Physics, Prague, Czech Republic
Institute of Physics of the Academy of Sciences of the Czech Republic, Prague, Czech Republic
Palacky University, RCPTM, Olomouc, Czech Republic
Institut de Physique Nucléaire d’Orsay (IPNO), Université Paris 11, CNRSIN2P3, Orsay, France
Laboratoire de l’Accélérateur Linéaire (LAL), Université Paris 11, CNRSIN2P3, France
Laboratoire de Physique Nucléaire et de Hautes Energies (LPNHE), Universités Paris 6 et Paris 7, CNRSIN2P3, Paris, France
Laboratoire de Physique Subatomique et de Cosmologie (LPSC), Université Joseph Fourier Grenoble, CNRSIN2P3, Grenoble INP, France
Station de Radioastronomie de Nançay, Observatoire de Paris, CNRS/INSU, France
SUBATECH, École des Mines de Nantes, CNRSIN2P3, Université de Nantes, France
Bergische Universität Wuppertal, Wuppertal, Germany
Karlsruhe Institute of Technology  Campus North  Institut für Kernphysik, Karlsruhe, Germany
Karlsruhe Institute of Technology  Campus North  Institut für Prozessdatenverarbeitung und Elektronik, Karlsruhe, Germany
Karlsruhe Institute of Technology  Campus South  Institut für Experimentelle Kernphysik (IEKP), Karlsruhe, Germany
MaxPlanckInstitut für Radioastronomie, Bonn, Germany
RWTH Aachen University, III. Physikalisches Institut A, Aachen, Germany
Universität Hamburg, Hamburg, Germany
Universität Siegen, Siegen, Germany
Dipartimento di Fisica dell’Università and INFN, Genova, Italy
Università dell’Aquila and INFN, L’Aquila, Italy
Università di Milano and Sezione INFN, Milan, Italy
Università di Napoli ”Federico II” and Sezione INFN, Napoli, Italy
Università di Roma II ”Tor Vergata” and Sezione INFN, Roma, Italy
Università di Catania and Sezione INFN, Catania, Italy
Università di Torino and Sezione INFN, Torino, Italy
Dipartimento di Matematica e Fisica ”E. De Giorgi” dell’Università del Salento and Sezione INFN, Lecce, Italy
Istituto di Astrofisica Spaziale e Fisica Cosmica di Palermo (INAF), Palermo, Italy
Istituto di Fisica dello Spazio Interplanetario (INAF), Università di Torino and Sezione INFN, Torino, Italy
INFN, Laboratori Nazionali del Gran Sasso, Assergi (L’Aquila), Italy
Benemérita Universidad Autónoma de Puebla, Puebla, Mexico
Centro de Investigación y de Estudios Avanzados del IPN (CINVESTAV), México, Mexico
Universidad Michoacana de San Nicolas de Hidalgo, Morelia, Michoacan, Mexico
Universidad Nacional Autonoma de Mexico, Mexico, D.F., Mexico
IMAPP, Radboud University Nijmegen, Netherlands
Kernfysisch Versneller Instituut, University of Groningen, Groningen, Netherlands
Nikhef, Science Park, Amsterdam, Netherlands
ASTRON, Dwingeloo, Netherlands
Institute of Nuclear Physics PAN, Krakow, Poland
University of Łódź, Łódź, Poland
LIP and Instituto Superior Técnico, Technical University of Lisbon, Portugal
’Horia Hulubei’ National Institute for Physics and Nuclear Engineering, Bucharest Magurele, Romania
University of Bucharest, Physics Department, Romania
University Politehnica of Bucharest, Romania
J. Stefan Institute, Ljubljana, Slovenia
Laboratory for Astroparticle Physics, University of Nova Gorica, Slovenia
Institut de Física Corpuscular, CSICUniversitat de València, Valencia, Spain
Universidad Complutense de Madrid, Madrid, Spain
Universidad de Alcalá, Alcalá de Henares (Madrid), Spain
Universidad de Granada and C.A.F.P.E., Granada, Spain
Universidad de Santiago de Compostela, Spain
School of Physics and Astronomy, University of Leeds, United Kingdom
Argonne National Laboratory, Argonne, IL, USA
Case Western Reserve University, Cleveland, OH, USA
Colorado School of Mines, Golden, CO, USA
Colorado State University, Fort Collins, CO, USA
Colorado State University, Pueblo, CO, USA
Fermilab, Batavia, IL, USA
Los Alamos National Laboratory, Los Alamos, NM, USA
Louisiana State University, Baton Rouge, LA, USA
Michigan Technological University, Houghton, MI, USA
New York University, New York, NY, USA
Northeastern University, Boston, MA, USA
Ohio State University, Columbus, OH, USA
Pennsylvania State University, University Park, PA, USA
University of Chicago, Enrico Fermi Institute, Chicago, IL, USA
University of Hawaii, Honolulu, HI, USA
University of Nebraska, Lincoln, NE, USA
University of New Mexico, Albuquerque, NM, USA
University of Wisconsin, Madison, WI, USA
University of Wisconsin, Milwaukee, WI, USA
Institute for Nuclear Science and Technology (INST), Hanoi, Vietnam(‡) Deceased
(a) Now at Konan University
(b) Also at the Universidad Autonoma de Chiapas on leave of absence from Cinvestav
(c) Now at University of Maryland
(d) Now at NYU Abu Dhabi
Keywords: cosmic ray experiments, ultra high energy cosmic rays
Contents
1 Introduction
The most commonly used shower observables for the study of the composition of Ultra High Energy Cosmic Rays (UHECR) are the mean value of the depth of shower maximum, , and its dispersion, . Inferring the mass composition from these measurements is subject to some level of uncertainty. This is because their conversion to mass relies on the use of shower simulation codes which include the assumption of a hadronic interaction model. The various interaction models [1] have in common the ability to fit lower energy accelerator data. However, different physical assumptions are used to extrapolate these low energy interaction properties to higher energies. Consequently they provide different expectations for and . The first aim of this paper is to discuss how the mean value of the depth of shower maximum and its dispersion can be used to interpret mass composition even in the presence of uncertainties in the hadronic interaction modeling.
Furthermore, we discuss the different roles of the two observables, and , with respect to mass composition. In the interpretation of data they are often used as different, and independent, aspects of the same phenomenon. However it is not true to say that both parameters reflect the cosmic ray composition to the same extent. According to the superposition model [2] is linear in and therefore it actually measures mass composition for both pure and mixed compositions. But, we will show that the behaviour of is more complex to interpret as there is no onetoone correspondence between its value and a given mean log mass. Only in the case of pure composition is this correspondence unique.
In this paper we refine the analysis method originally proposed by Linsley [3, 4] and apply it to the Auger data. The Pierre Auger Collaboration has published results on the mean and dispersion of the distribution at energies above eV [5, 6]. In this work we apply the proposed method to convert those observables to the first moments of the log mass distribution, namely and .
The paper is organized as follows. In Section 2 we discuss the parameterization for and . In Sec. 3 we test the method with shower simulations assuming different mass distributions. Sec. 4 describes the application of the method to data. The discussion of the results and the conclusions follow in sections 5 and 6 respectively. The details of the parameterization and the best fit values for the hadronic interaction models are summarized in Appendix A.
2 A method to interpret and
The interpretation of and can be simplified by making use of an analysis method based on the generalized Heitler model of extensive air showers [7]. In this context is a linear function of the logarithm of the shower energy per nucleon:
(2.1) 
where is the mean depth of proton showers at energy and is the elongation rate [8, 9, 10], i.e., the change of per decade of energy. The High Energy hadronic interaction models used in this work are EPOS 1.99 [11], Sibyll 2.1 [12], QGSJet 01 [13] and QGSJet II [14]. Simulated data show that eq. (2.1) gives a fair description of EPOS and Sibyll results in the full range of interest for this work, 10 to 10 eV, but does not reproduce accurately QGSJet models. For this reason we generalize the original representation as:
(2.2) 
where the parameters and are expected to be zero if the model predictions are compatible with the superposition result (2.1).
For nuclei of the same mass one expects the shower maximum to be on average:
(2.3) 
and its dispersion to be only influenced by showertoshower fluctuations:
(2.4) 
Here denotes the mean depth at maximum of proton showers, as obtained from either eq. (2.1) or (2.2), and is the variance for mass , . The energy dependent parameter appearing in (2.3) is:
(2.5) 
The values of the parameters , , , depend on the specific hadronic interaction model. In this work they are obtained from CONEX [15] shower simulations as described in Appendix A.
In the case of a mixed composition at the top of the atmosphere, the mean and variance of depend on the distribution. There are two independent sources of fluctuations: the intrinsic showertoshower fluctuations and the dispersion arising from the mass distribution. The first term gives rise to , the average variance of weighted according to the distribution. The second contribution can be written as where is the variance of the distribution. We can finally write for the two profile observables:
(2.6) 
(2.7) 
The two equations depend on energy through the parameters but also via and the possible dependence of the two moments of the distribution.
To obtain an explicit expression for we need a parameterization for . We assume a quadratic law in :
(2.8) 
where is the variance for proton showers. The evolution of with energy is included in and the parameter :
(2.9) 
The parameters , , , , , depend on hadronic interactions: the values used in the paper are given in Appendix A.
Using measurements of and , equations (2.6) and (2.7) can be inverted to get the first two moments of the distribution. From eq. (2.6) one gets:
(2.10) 
Averaging eq. (2.8) on one obtains:
(2.11) 
Substituting in eq. (2.7) we get:
(2.12) 
But by definition . Solving in one finally obtains:
(2.13) 
Equations (2.10) and (2.13) are the key tools used throughout this work for interpreting Pierre Auger Observatory data in terms of mass composition and assessing the validity of available hadronic interaction models.
3 Testing the method with simulation
Equations (2.6) and (2.7) can be tested with simulations. They contain parameters depending on the hadronic interaction properties and on the mass distribution of nuclei. The mass distribution of nuclei refers to those nuclei hitting the Earth’s atmosphere: it does not matter what is the source of the mass dispersion, either a mixed composition at injection or the dispersion caused by propagation. So, in order to test the method we will simply use different test distributions of the masses at the top of the atmosphere.
For this purpose we have chosen three different mass distributions:

A distribution uniform in from to and independent of energy. The values of and are respectively 2.01 and 1.16.

A Gaussian distribution with increasing linearly with from at 10 eV to at 10 eV and independent of energy. The Gaussian is truncated to less than 2 sigmas to avoid unphysical mass values. In this case the dispersion is fixed and equal to 0.66 but varies with energy.

Two masses, H and Fe, with proton fraction decreasing linearly with from 1 at 10 eV to 0 at 10 eV. In this case, both and vary with energy.
Figure 1 shows the result of the test for the three mass distribution hypotheses. To generate the distributions we have used CONEX [15] showers with Sibyll 2.1 [12] as the hadronic interaction model. These distributions do not include detector effects. For each test mass hypothesis, the mean and RMS are retrieved from the resulting distribution obtained from the simulations. These are shown as full circles, and in left and right panels respectively. The dashed lines are calculated using equations (2.6) and (2.7) for the three different mass hypotheses by using only the first two moments and .
One can see that, despite the simple assumptions made, good agreement is achieved for all the three mass distributions. The dotdashed line refers to the contribution of the first term in eq. (2.7). The comparison between the two lines (dashed vs. dotdashed) highlights how different the interpretation of data can be if one does not take into account the mass dispersion term.
The inverse equations (2.10) and (2.13) have also been tested using Monte Carlo simulation. In this case and have been obtained as a function of directly from the input mass distributions. These values are shown as full circles in Figure 2. The and retrieved from the corresponding distributions are used in equations (2.10) and (2.13) to get and . These are shown in Fig. 2 as dashed lines. Also in this case, the comparison is quite successful.
The simulated data sample can also be used to estimate the systematic uncertainty in the calculation of the moments of the () distribution induced by the missing knowledge of the hadronic interaction mechanism. This study is pursued using simulated showers generated with a given model together with parameters of another model in equations (2.6) and (2.7) for the profile variables, and (2.10) and (2.13) for the log mass variables. An example of this procedure is shown in Fig. 2 where the dotted lines show the calculation with the parameters of QGSJet II and the full circles refer to data simulated with Sibyll 2.1. As a summary of these crossmodel checks, we find mean absolute deviations of 4 to 27 g cm for and 1 to 5.4 g cm for , where the maximum deviations are obtained crossing EPOS with QGSjetII. The same study done for the moments of the log mass distribution gives mean absolute deviations of 0.2 to 1.2 for and 0.02 to 0.5 for . In this case the maximum values refer to EPOS vs. QGSJet 01 for the first moment and QGSJet II vs. QGSJet 01 for the second.
4 Application to data
At ultrahigh energies, shower development can be directly measured using fluorescence and Cherenkov light profiles. Mean data as a function of energy are available from Fly’s Eye [16], HiRes [17, 18], Auger [5], Yakutsk [19] and Telescope Array [20]. data were complemented with fluctuation measurements as early as 1980s (see e.g. [21] and references therein) but only recently have precise optical detector measurements become available [18, 5, 19].
The Pierre Auger Collaboration has published results on the mean and dispersion of the distribution at energies above 10 eV [5]. Here we apply the method presented in this work to an updated dataset available in [6, 22]. These data are shown in Figure 3.
In the Auger analysis [5], the events are selected using fiducial volume cuts based on the shower geometry. This ensures that the viewable range for each shower is large enough to accommodate the full distribution. Also, the detector resolution is accounted for by subtracting in quadrature its contribution to the measured dispersion. This allows the direct conversion to the moments of the distribution using equations (2.10) and (2.13) without the need of more complex treatment, such as is required in the presence of acceptance biases [23, 24].
The moments of the log mass distribution, and , as obtained using equations (2.10) and (2.13), are shown (full circles) as a function of in Figures 4 and 5 respectively. Error bars show the statistical errors obtained from the propagation of data errors and the errors of the fitted parameters. Shaded bands are the systematic uncertainties obtained by summing in quadrature the different individual contributions.
The systematic uncertainties on and data points have different sources: calibration, atmospheric conditions, reconstruction and event selection [5]. Another source of systematics is related to the uncertainty of the FD energy scale [25], 22 %, which induces an uncertainty in and via the parameters of the models. All these uncertainties contribute approximately at the same level and independently of energy. The figures show the results for the moments of the log mass distribution for EPOS 1.99 [11], Sibyll 2.1 [12], QGSJet 01 [13] and QGSJet II [14].
Despite the uncertainties and the different mass offsets of the models, the overall features are similar in all the cases. So far as the energy dependence is concerned, the data imply an increasing above 10 eV from light to intermediate masses and a decreasing over the whole energy range.
Looking more specifically to the different hadronic models we notice a slight change in the log mass scale. The highest masses are obtained for EPOS 1.99. Sibyll 2.1 and QGSJet II show intermediate values, whereas the lowest masses are obtained for QGSJet 01. In particular at the mean log mass, , is 1.10, 0.70, 0.60 and 0.12 respectively for EPOS 1.99, Sibyll 2.1, QGSJet II and QGSJet 01 with statistical errors of about 0.08 and systematic uncertainty of about 0.6. The Pierre Auger Collaboration has recently published the measurement of the protonair cross section for the energy interval 10 to 10 eV [26]. That measurement is done using the showers with g cm, corresponding to 20% of the total distribution. Even in the most unfavourable case, (the and predicted by EPOS), one finds that several realizations obtained from the allowed and have enough protons in the most deeply penetrating showers to fulfill the selection criteria adopted in the Auger analysis.
Whereas always has valid values (apart a small region which crosses for QGSJet 01), there are wide energy intervals where is negative. Considering eq. (2.13) one can see that these values occur for energies where the shower fluctuations corresponding to the mean log mass exceed the measured fluctuations. Figure 5 shows that data points are within the allowed physical region only for EPOS 1.99 and Sibyll 2.1. They are partly outside for QGSJet II, and completely outside for QGSJet 01. However the current systematic uncertainties do not allow one to establish stringent tests to the models.
The method presented in this work shows that the Pierre Auger Observatory data can confront hadronic physics models provided that future developments in the shower data analysis reduce systematics. By shrinking the shaded bands in Figure 5 it will be possible to constrain those models.
5 Discussion
The importance of the combined study of the mean values and fluctuations of mass dependent observables has been addressed by several authors [21, 27, 28, 3, 4]. In particular, Linsley [4] showed that a combined analysis of the mean and the variance of can provide a useful representation of the mass transition (if any) to be found in shower profile data. In fact, possible transitions are constrained to a limited region of the (, ) plane. More recently a similar study using the  correlation^{1}^{1}1 In this case the dependence on hadronic models has been accounted for by subtracting the corresponding observables predicted by the models for iron. reached a similar conclusion [29].
Converting data to variables, as described in Sec. 2, one can plot Pierre Auger Observatory data in the (, ) plane. Since this procedure depends on the hadronic model, one gets a plot for each model as shown in Figure 6. Data points are shown as full circles with size increasing in proportion to . The error bars are tilted because of correlations arising from equations (2.10) and (2.13) and represent the principal axes of the statistical error ellipses. The solid lines show the systematic uncertainties. The same figure shows the region allowed for mass compositions. The contour of this region (gray thick line) is generated by mixing neighbouring nuclei in the lower edge and extreme nuclei (protons and iron) in the upper edge. Each of these mixings is an arch shaped line in the (, ) plane.
Figure 6 shows that the Auger data lie outside the allowed boundaries for part of the energy range in some of the models. As noted previously, systematic uncertainties are still large and thus prevent us from more definite conclusions. However the energy evolution is common to all models suggesting that the average mass increases with decreasing log mass dispersion. This behaviour might imply astrophysical consequences.
In fact there are only a few possibilities for extragalactic source models to produce compositions with small log mass dispersion at the Earth. Protons can traverse their path from sources to the Earth without mass dispersion, but this case is excluded by Pierre Auger Observatory data at the highest energies.
Nuclei originating from nearby sources ( 100 Mpc) might be detected with small mass dispersion. For these sources, propagation does not degrade mass and energy so the spectrum and composition reflect closely their values at injection. But, if sources are distributed uniformly, distant sources induce natural mass dispersions. Small dispersions are possible only when there is small observed mass mixing so that, at each energy, only nuclei with a small spread in masses are present. This corresponds to the low edge of the contour of the allowed region in the (, ) plane.
Protons originating by the photodisintegration of nuclei are the main source of mass dispersion because they populate each energy region. The possible end of the injection spectrum based on a rigiditydependent mechanism can reduce the proton component at high energies, thus producing a reduction of the mass dispersion at the highest energies. A complete study of source models under several hypotheses is required to study all the source parameters that limit the mass dispersion in the propagation of extragalactic cosmic rays. Recent studies, see e.g. [30, 31], based on the assumption of a uniform source distribution, have shown that the Auger composition results, when combined with the energy spectrum, require hard injection spectra (i.e. index < 2) with low energy cutoffs and the possible presence of local sources.
6 Conclusions
In this work we presented a method for interpreting and in terms of mass composition. The method is based on an extension of the Heitler model of extensive air showers. The parameterization given in equations (2.6) and (2.7) expresses those two profile observables as a linear combination of the first two moments of the log mass distribution, and , and of the mean shower fluctuations.
We first note that the method provides an effective key to the interpretation of data. The energy dependences of and are sometimes considered as different expressions of the same physical features, e.g. an increase or decrease of the mean log mass. However their different meanings can be easily understood by looking at the dependence on the mass variables. At a fixed energy is only function of ; therefore, it only carries information of the average composition. However, cannot be interpreted as a measure of the average composition since it is also affected by the log mass dispersion. Similarly, the inference of hadronic interaction properties from can be wrong unless the mass dispersion term () is negligible. The parameter represents the dispersion of the masses as they hit the Earth atmosphere. It reflects not only the spread of nuclear masses at the sources but also the modifications that occur during their propagation to the Earth.
The method has been succesfully tested, with the simulation of different mass distributions in the energy interval from 10 to 10 eV showing the robustness of the parameterization. We have applied the method to the Pierre Auger Observatory data to get the first two moments of the distribution. The outcome relies on the choice of a hadronic interaction model to set the parameters and the appropriate shower fluctuations. Four models have been used, EPOS 1.99, Sibyll 2.1, QGSJet 01 and QGSJet II, and the corresponding moments of the log mass distribution have been obtained as a function of energy. Despite the differences in the chosen models, the overall features are quite similar. In particular we find an increasing above 10 eV from light to intermediate masses and a decreasing over the whole energy range, while the mean log mass scale changes with hadronic models.
The results presented in this paper show the capability of the method to infer important features of the mass distribution of the UHECR nuclei. This is a remarkable outcome with respect to the study of source scenarios and propagation. In fact we do not only access the average mass, but also the mass dispersion. While a pure proton beam at the sources is not changed by propagation, nuclei should increase the mass dispersion in their path towards the Earth. The Auger results seem to imply either closeby sources or hard spectral indices, if the energy evolution of the present hadronic interaction models can be trusted.
The proposed method can also be used as a tool to investigate the validity of hadronic interaction models. In particular it has been shown that the intrinsic shower fluctuations are sometimes larger than the measured dispersions. This happens in different energy intervals for the different models. At the highest energies, all models approach the lower boundary, and some of them enter the unphysical region, but the current systematic uncertainties prevent us from confidently rejecting any model. Provided that systematic uncertainties can be reduced in future data analysis, the method can be used to constrain hadronic interaction models. The addition of new measurements, such as the muon content of EAS [32, 33], may allow us to place stronger bounds to the models.
Appendix A Parameterization of the shower mean depth and its fluctuations
The shower code chosen for this work is CONEX [15]. CONEX is a hybrid simulation code that is suited for fast onedimensional simulations of shower profiles, including fluctuations. It combines Monte Carlo simulation of high energy interactions with a fast numerical solution of cascade equations for the resulting distributions of secondary particles. In our CONEX simulation we used the default energy thresholds settings of version v2r3.1^{2}^{2}2 The hadron and muon cutoff (minimum energy) is 1 GeV, the cutoff for electrons, positrons and gammas (e/m particles) is 1 MeV, the threshold energy for solving cascade equations is 0.05, 0.0005 and 0.005 of the primary energy for hadrons, muons and e/m particles respectively. The abovethreshold e/m interaction are simulated with the EGS4 program [34]. The low energy (E 80 GeV) interaction model is GEISHA [35] .
The parameters , , and used in equations (2.6) and (2.7) have been obtained by fitting CONEX showers for four different primaries (H, He, N and Fe) in nine energy bins of width = 0.25 ranging from 10 to 10 eV, and for all the hadronic models used in this work: EPOS 1.99 [11], Sibyll 2.1 [12], QGSJet 01 [13] and QGSJet II [14]. In total, about 25,000 showers have been used for each energy bin and for each hadronic model. The fit procedure always converges with mean (maximum) residuals from the simulated data of about 1 (3) g cm for all the models. The best fit values are reported in Table 1 with their errors.
parameter  EPOS 1.99  Sibyll 2.1  QGSJet 01  QGSJet II 

809.7 0.3  795.1 0.3  774.2 0.3  781.8 0.3  
62.2 0.5  57.7 0.5  49.7 0.5  45.8 0.5  
0.78 0.24  0.04 0.24  0.30 0.24  1.13 0.24  
0.08 0.21  0.04 0.21  1.92 0.21  1.71 0.21  
Shower variances have been fitted using the parameterization given in equations (2.8) and (2.9) and the same simulated data set described above. The mean (maximum) residuals from the simulated data are about 1 (3) g cm for all the models. The best fit parameters are given in Table 2 with their errors.
parameter  EPOS 1.99  Sibyll 2.1  QGSJet 01  QGSJet II 

3279 51  2785 46  3852 55  3163 49  
47 66  364 58  274 70  237 61  
228 108  152 93  169 116  60 100  
0.461 0.006  0.368 0.008  0.451 0.006  0.386 0.007  
0.0041 0.0016  0.0049 0.0023  0.0020 0.0016  0.0006 0.0021  
0.059 0.002  0.039 0.002  0.057 0.001  0.043 0.002 
Acknowledgments
The successful installation, commissioning, and operation of the Pierre Auger Observatory would not have been possible without the strong commitment and effort from the technical and administrative staff in Malargüe.
We are very grateful to the following agencies and organizations for financial support: Comisión Nacional de Energía Atómica, Fundación Antorchas, Gobierno De La Provincia de Mendoza, Municipalidad de Malargüe, NDM Holdings and Valle Las Leñas, in gratitude for their continuing cooperation over land access, Argentina; the Australian Research Council; Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), Financiadora de Estudos e Projetos (FINEP), Fundação de Amparo à Pesquisa do Estado de Rio de Janeiro (FAPERJ), Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP), Ministério de Ciência e Tecnologia (MCT), Brazil; AVCR AV0Z10100502 and AV0Z10100522, GAAV KJB100100904, MSMTCR LA08016, LG11044, MEB111003, MSM0021620859, LA08015, TACR TA01010517 and GA UK 119810, Czech Republic; Centre de Calcul IN2P3/CNRS, Centre National de la Recherche Scientifique (CNRS), Conseil Régional IledeFrance, Département Physique Nucléaire et Corpusculaire (PNCIN2P3/CNRS), Département Sciences de l’Univers (SDUINSU/CNRS), France; Bundesministerium für Bildung und Forschung (BMBF), Deutsche Forschungsgemeinschaft (DFG), Finanzministerium BadenWürttemberg, HelmholtzGemeinschaft Deutscher Forschungszentren (HGF), Ministerium für Wissenschaft und Forschung, NordrheinWestfalen, Ministerium für Wissenschaft, Forschung und Kunst, BadenWürttemberg, Germany; Istituto Nazionale di Fisica Nucleare (INFN), Ministero dell’Istruzione, dell’Università e della Ricerca (MIUR), Italy; Consejo Nacional de Ciencia y Tecnología (CONACYT), Mexico; Ministerie van Onderwijs, Cultuur en Wetenschap, Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO), Stichting voor Fundamenteel Onderzoek der Materie (FOM), Netherlands; Ministry of Science and Higher Education, Grant Nos. N N202 200239 and N N202 207238, Poland; Portuguese national funds and FEDER funds within COMPETE  Programa Operacional Factores de Competitividade through Fundação para a Ciência e a Tecnologia, Portugal; Romanian Authority for Scientific Research ANCS, CNDIUEFISCDI partnership projects nr.20/2012 and nr.194/2012, project nr.1/ASPERA2/2012 ERANET and PNIIRUPD20113014517, Romania; Ministry for Higher Education, Science, and Technology, Slovenian Research Agency, Slovenia; Comunidad de Madrid, FEDER funds, Ministerio de Ciencia e Innovación and ConsoliderIngenio 2010 (CPAN), Xunta de Galicia, Spain; The Leverhulme Foundation, Science and Technology Facilities Council, United Kingdom; Department of Energy, Contract Nos. DEAC0207CH11359, DEFR0204ER41300, DEFG0299ER41107, National Science Foundation, Grant No. 0450696, The Grainger Foundation USA; NAFOSTED, Vietnam; Marie CurieIRSES/EPLANET, European Particle Physics Latin American Network, European Union 7th Framework Program, Grant No. PIRSES2009GA246806; and UNESCO.
References
 [1] For a recent review see e.g. R. Engel, D. Heck and T. Pierog, Extensive Air Showers and Hadronic Interactions at High Energy, Ann. Rev. Nucl. Part. Sci. 61 (2011) 467.
 [2] See e.g. T. K. Gaisser, Cosmic Rays and Particle Physics, Cambridge University Press, Cambridge, 1990.
 [3] J. Linsley, Spectra, anisotropies and composition of cosmic rays above 1000 GeV, rapporteur paper in Proc. 18th International Cosmic Ray Conference (ICRC), Bangalore, India, 12 (1983) 135.
 [4] J. Linsley, Protonair and protonproton cross sections from air shower data, in Proc. 19th International Cosmic Ray Conference (ICRC), San Diego, California, 6 (1985) 1.
 [5] J. Abraham et al. (Pierre Auger Collaboration), Measurement of the Depth of Maximum of Extensive Air Showers above 10 eV, Phys. Rev. Lett. 104 (2010) 091101.
 [6] P. Facal San Luis for the Pierre Auger Collaboration, The distribution of shower maxima of UHECR air showers, in Proc. 32nd International Cosmic Ray Conference (ICRC), Beijing, China, 2 (2011) 225 and arXiv:1107.4804v1 [astroph.HE].
 [7] See e.g. J. Matthews, A Heitler model of extensive air showers, Astropart. Phys. 22 (2005) 387 and references therein.
 [8] J. Linsley, Structure of large air showers at depth 834 g/cm, in Proc. 15th International Cosmic Ray Conference (ICRC), Plovdiv, Bulgaria, 12 (1977), 89.
 [9] T. K. Gaisser et al., Elongation Rate of Air Showers and Implications for 1010 eV Particle Interactions, in Proc. 16th International Cosmic Ray Conference (ICRC), Kyoto, Japan, 9 (1979), 275.
 [10] J. Linsley and A. A. Watson, Validity of scaling to 10 eV and highenergy cosmic ray composition, Phys. Rev. Lett. 46 (1981), 459.
 [11] T. Pierog and K. Werner, Muon Production in Extended Air Shower Simulations, Phys. Rev. Lett., 101 (2008) 171101.
 [12] E. J. Ahn et al., Cosmic ray interaction event generator SIBYLL 2.1, Phys. Rev. D 80 (2009) 094003.
 [13] N. N. Kalmykov et al., Quarkgluon string model and EAS simulation problems at ultrahigh energies, Nucl. Phys. B (Proc. Suppl.) 52 (1997), 17.
 [14] S. Ostapchenko, Nonlinear screening effects in high energy hadronic interactions, Phys. Rev., D 74 (2006) 014026.
 [15] T. Pierog et al., First results of fast onedimensional hybrid simulation of EAS using CONEX, Nucl. Phys. B (Proc. Suppl.) 151 (2006) 159.
 [16] D. J. Bird et al., Evidence for correlated changes in the spectrum and composition of cosmic rays at extremely highenergies, Phys. Rev. Lett. 71 (1993) 3401.
 [17] R. U. Abbasi et al. (HiRes Collaboration), A Study of the composition of ultrahigh energy cosmic rays using the High Resolution Fly’s Eye, Astroph. J. 622 (2005) 910;.
 [18] R. U. Abbasi et al. (HiRes Collaboration), Indications of ProtonDominated Cosmic Ray Composition above 1.6 EeV, Phys. Rev. Lett. 104 (2010) 161101.
 [19] S. P. Knurenko and A. Sabourov, Spectrum and composition of cosmic rays in the energy range 10  10 eV derived from the Yakutsk array data, in Proc. 32nd International Cosmic Ray Conference (ICRC), Beijing, China, 1 (2011) 189.
 [20] C. Jui et al. (Telescope Array Collaboration), Cosmic Ray in the Northern Hemisphere: Results from the Telescope Array Experiment, Proc. APS DPF Meeting, (2011) Providence, RI, USA arXiv:1110.0133.
 [21] R. Walker and A. A. Watson, Measurement of the fluctuations in the depth of maximum of showers produced by primary particles of energy greater than eV, J. Phys. G 8 (1982) 1131.

[22]
D. GarciaPinto, for the Pierre Auger Collaboration,
Measurements of the Longitudinal Development of Air Showers with the Pierre Auger Observatory, in Proc. 32nd International Cosmic Ray Conference (ICRC), Beijing, China, 2 (2011) 87 and
arXiv:1107.4804v1 [astroph.HE].
http://www.auger.org/technical_info/ICRC2011/shower_development_data.txt  [23] M. Unger, EAS Studies of Cosmic Rays above 10 eV, rapporteur paper in Proc. 32nd International Cosmic Ray Conference (ICRC), Beijing, China, 12 (2011) 225.
 [24] L. Cazon and R. Ulrich, The nonlinearity between and induced by the acceptance of fluorescence telescopes, Astropart. Phys., 38 (2012) 41.
 [25] R. Pesce, for the Pierre Auger Collaboration, Energy calibration of data recorded with the surface detectors of the Pierre Auger Observatory: an update, in Proc. 32nd International Cosmic Ray Conference (ICRC), Beijing, China, 2 (2011) 214 and arXiv:1107.4809 [astroph.HE].
 [26] P. Abreu et al. (Pierre Auger Collaboration), Measurement of the protonair crosssection at = 57 TeV with the Pierre Auger Observatory, Phys. Rev. Lett. 109 (2012) 062002.
 [27] A. A. Watson and J. G. Wilson, Fluctuation studies of large air showers: the composition of primary cosmic ray particles of energy E 10 eV, J. Phys. A 7 (1974) 1199.
 [28] R. Walker and A. A. Watson, Measurement of the elongation rate of extensive air showers produced by primary cosmic rays of energy above 2 10 eV, J. Phys., G 7 (1981) 1297.
 [29] K. H. Kampert and M. Unger, Measurements of the Cosmic Ray Composition with Air Shower Experiments, Astropart. Phys. 35 (2012) 660.
 [30] A. M. Taylor, M. Ahlers and F. A. Aharonian, The need for a local source of UHE CR nuclei, Phys.Rev., D 84 (2011) 105007.
 [31] D. Allard, Extragalactic propagation of ultrahigh energy cosmicrays, Astropart. Phys. 3940 (2012) 33
 [32] J. Allen for the Pierre Auger Collaboration, Interpretation of the signals produced by showers from cosmic rays of 10 eV observed in the surface detectors of the Pierre Auger Observatory, in Proc. 32nd International Cosmic Ray Conference (ICRC), Beijing, China, 2 (2011) 83 and arXiv:1107.4804v1 [astroph.HE].
 [33] G. Rodriguez for the Pierre Auger Collaboration, Reconstruction of inclined showers at the Pierre Auger Observatory: implications for the muon contenty, in Proc. 32nd International Cosmic Ray Conference (ICRC), Beijing, China, 2 (2011) 95 and arXiv:1107.4809 [astroph.HE].
 [34] W. Nelson et al., SLAC265, Stanford Linear Accelerator Center (1985).
 [35] GEISHA, H. Fesefeldt, RWTH Aachen report PITHA 85/2 (1985).