On the measurement of the proton-air cross section using air shower data

The analysis of high-energy air shower data allows one to study the proton-air cross section at energies beyond the reach of fixed target and collider experiments. The mean depth of the first interaction point and its fluctuations are a measure of the proton-air particle production cross section. Since the first interaction point in air cannot be measured directly, various methods have been developed in the past to estimate the depth of the first interaction from air shower observables in combination with simulations. As the simulations depend on assumptions made for hadronic particle production at energies and phase space regions not accessible in accelerator experiments, the derived cross sections are subject to significant systematic uncertainties. The focus of this work is the development of an improved analysis technique that allows a significant reduction of the model dependence of the derived cross section at very high energy. Performing a detailed Monte Carlo study of the potential and the limitations of different measurement methods, we quantify the dependence of the measured cross section on the hadronic interaction model used. Based on these results, a general improvement of the analysis methods is proposed by introducing the actually derived cross section already in the simulation of reference showers. The reduction of the model dependence is demonstrated for one of the measurement methods.


Introduction
The natural beam of cosmic ray particles extends to energies far beyond the reach of any Earth-based particle accelerator.Therefore cosmic ray data provide a unique window to study hadronic interaction phenomena at energies up to several Joules per particle, corresponding to an equivalent center-of-mass energy of up to 450 TeV.
On the other hand, direct observation of the first interaction of ultra-high energy cosmic rays in the upper atmosphere is impossible due to the very low flux of these particles.Only the cascades of secondary particles, called extensive air showers (EAS), can be measured with arrays of particle detectors or optical telescopes.To obtain information on the first interactions in an air shower it is necessary to link the measured air shower characteristics to that of high energy particle production in the shower.This can be done with detailed Monte Carlo simulations of the shower evolution and the corresponding shower observables, but inevitably causes a dependence of the results on hadronic interaction models needed for the shower description.
The total particle production cross section is one of the most fundamental quantities that characterizes hadronic interactions.Considering proton-induced air showers of the same primary energy, the depth of the first interaction point, X 1 , is distributed according to dp dX 1 = 1 λ p−air e −X 1 /λ p−air , where λ p−air is the interaction mean free path of protons in air.The mean depth of the first interaction point and its shower-to-shower fluctuations are directly linked to the size of the proton-air cross section by with the mean target mass of air being m air ≈ 14.45 m p = 24160 mb g/cm 2 [1].It is, therefore, not surprising that there is a long history of attempts to infer this cross section from high-energy cosmic ray data [2][3][4][5][6][7][8][9][10][11][12][13][14][15].A compilation of published proton-air cross section measurements and predictions of ultra-high energy hadronic interaction models is shown in Fig. 1.All data above 100 TeV are based on the analysis of EAS data [7][8][9][10][11][12][13]. Also the ARGO-YBJ measurements around 10 TeV are originating from a high altitude EAS array [14].The data at lower energies stem from unaccompanied hadron analyses [2][3][4][5][6].
All cosmic ray measurements of the proton-air cross section are only sensitive to the particle production cross section [9,22].In addition, interactions with insignificant particle production have no measurable impact on cosmic ray observations.This is the case for air shower based techniques as well as for the unaccompanied hadron method.
modeled as precisely as possible.Sophisticated EAS Monte Carlo simulation programs are available for this task.This results in a highly indirect analysis procedure.
In this article we will perform a detailed Monte Carlo study of different methods to derive the proton-air cross section from air shower data at an energy of ∼ 10 19 eV.We will assume protons as cosmic ray particles in the energy range considered here.Although there exist theoretical models [23][24][25][26][27] and also experimental indications [28][29][30] for this flux being indeed dominated by protons, other elements in the primary flux have also to be taken into account.This will be done in a forthcoming article that will address specifically this topic [31].
Based on the results of the Monte Carlo study of existing analysis methods a general improvement is proposed and explicitly applied to one of the methods.By accounting for the actually measured cross section already in the simulation of showers, the systematic uncertainty due to the limited knowledge of hadronic interactions is significantly reduced.The performance of the improved method is thoroughly tested for the application to high quality data of the depth of the shower maximum.Sources of systematic uncertainties of the resulting cross section are discussed.It is shown that the dependence on the hadronic interaction model, the most important source of systematic uncertainties, can be significantly reduced by incorporating the measured cross section in a consistent way in the shower simulation.
The analyses of this work are done at an energy of 10 19 eV, at which high statistics Definition of variables to characterize EAS longitudinal shower development.Zero slant depth is where the cosmic ray particle enters the atmosphere.The first interaction occurs at X 1 .At X max the shower reaches its maximum particle number N max .After the maximum the shower is attenuated over ∆X 2 before it reaches the particle number N stage at the slant depth of X obs .
data are available or expected from HiRes [32], the Pierre Auger Observatory [33], and Telescope Array [34], but the results are qualitatively also valid at other energies.The expected reduction of the model dependence will be smaller at energies where data from accelerators are available, i.e. 10 15 eV and below.

Relation between air shower observables and depth of the first interaction point
Interactions over many decades in energy are occurring during EAS development.In the startup phase of the shower, relatively few ultra-high energy hadronic interactions are distributing the energy of the primary cosmic ray particle to a quickly growing number of secondary particles.The stochastic nature of these initial interactions is the source of strong fluctuations of EAS initiated by identical cosmic ray particles (shower-toshower fluctuations).A significant part of secondary particles decay to electromagnetic particles and lead to the development of an electromagnetic cascade that, already after a few interactions, contains most of the total EAS energy and constitutes by far the largest fraction of the particles.The interactions of particles in the e.m. cascade are theoretically well understood.Also the large number of these interactions levels out additional large-scale fluctuations .Thus, the development of an air shower can be characterized by two main stages: (i) Startup phase, consisting of few initial hadronic interactions at ultra-high energies.These interaction are the main source of fluctuations.
An exception are electromagnetic showers of E > 10 18 eV that, by chance, happen to start very deep in the atmosphere, for which the Landau-Pomeranchuk-Migdal effects leads to very large shower-toshower fluctuations [35][36][37].
(ii) Cascade phase, in which very numerous interactions of the bulk of the shower particles at intermediate energies take place.No significant large-scale fluctuations are expected from this part.
A clear boundary between the two stages cannot be drawn.The transition is seamless and by itself subject to strong shower-to-shower fluctuations.Only treatment with full air shower Monte Carlo simulation programs can fully account for all fluctuations.
To discuss the relation of the depth of the first interaction point to air shower observables it is useful to introduce a simple longitudinal model of air shower development that connects the proton-air cross section to the observables in a transparent way.The naming conventions for the model are given in Fig. 2. We distinguish between two types of air shower observation: by ground based detector arrays and the direct observation of longitudinal shower profiles by telescope detectors.The relevant air shower observables used in proton-air cross section analyses are the position of the shower maximum, X max , or the combination of the number of electrons N e and muons N µ at the atmospheric slant depth of the detector X obs .
In the following, the correlation of these EAS observables to the characteristics of the ultra-high energy interactions is studied with the one dimensional air shower Monte Carlo program CONEX v2r2 [38].
The particle numbers are defined as the total number of electrons above 1 MeV and muons above 1 GeV.This corresponds to typical quantities observed in air shower arrays, however the exact definition of observables depends very much on the experimental setup and varies strongly from experiment to experiment.By using the total number of particles for our study we are focusing on general air shower properties.
The shower maximum X max is the slant depth of the maximum energy deposit of the shower in the atmosphere.This definition matches the shower maximum derived from the fluorescence light profile of showers and coincides with that of the particle number within ∼ 3 gcm −2 .

Arrays of particle detectors
Using ground based air shower arrays, one can estimate the proton-air cross section by measuring the frequency of air showers of the same energy and stage of their development at different atmospheric depths.By selecting events of the same energy but different directions of incidence, the point of the first interaction has to vary with the angle, in order to observe the shower at the same stage of development.The selection of showers of constant energy and stage can be done only approximately and depends on the particular detector setup, but the typical requirement is a constant (N e , N µ ) at observation level.Examples of measurements of this type are [7][8][9][10]14].By requiring a given number of muons at the detector level does, in first approximation, select EAS of the same primary energy, because the attenuation of muons is small.However, also muons are slowly attenuated in the atmosphere.To correct for that a constant intensity selection can be applied [39].Showers with identical primary energy at the same stage of their shower development are assumed to yield the same number of electrons since, after the shower maximum, the electromagnetic shower attenuation is approximately universal [40][41][42][43][44][45].
The number of showers N sh selected per area dA and time dt, requiring a constant (N e , N µ ) at the atmospheric depth X obs can be written as dN sh dN e dN µ dX obs dA dt = dΩ dE dX where dN CR /(dΩ dE dA dt) is the flux of cosmic ray particles.The distribution P ∆X 1 = dp 1 /d∆X 1 describes shower-to-shower fluctuation of ∆X 1 , see Fig. 3.The twodimensional probability density P ∆X 2 = dp 2 /(d∆X 2 dN e ) is the frequency of showers of given energy E, X 1 , and ∆X 1 to have a shower size N e at a depth distance ∆X 2 from the shower maximum.It is required that the electron number after the shower maximum is attenuated to N e , while X obs = X 1 + ∆X 1 + ∆X 2 .The selection of showers by energy according to their muon number is reflected by the energy dependent probability distribution P Nµ = dp µ /dN µ .
It is important to notice that the distribution P ∆X 1 does not depend on the depth X 1 of the first interaction point.This is shown in Fig. 4 for one hadronic interaction model.At high energy, and for air densities of relevance, almost all charged secondary pions always interact and all neutral pions decay immediately.Furthermore, the electromagnetic cascade initiated by the photons from π 0 decay does not depend on the local air density.This makes the startup phase of a shower to a good approximation independent of X 1 .
The distribution P ∆X 1 depends, however, on the hadronic interaction model used for simulation.This is displayed in Fig. 3: both the mean ∆X 1 and the shape of  the distribution depend strongly on the chosen interaction model.In contrast, the distribution P ∆X 2 exhibits only a small model dependence due to the universality of electromagnetic air showers after shower maximum, see Fig. 5 (left).
An additional source of strong model dependence of Eq. ( 3) is the number of muons predicted for a given shower energy and depth.This can be seen in Fig. 5 (middle) where the frequency of finding a certain muon number for showers at 10 19 eV and different models is shown.The effect of this model dependence is displayed in Fig. 5 (right) by showing the folded distributions In real data analysis it must be taken into account that X obs , N e and N µ are only known with a limited precision due to the detector and shower reconstruction resolution.In the model described by Eq. ( 3) this would furthermore add the corresponding uncertainty distributions P S res and integrations over the three observables.The influence of the measurement resolution will be discussed below but is not included in (3) for sake of clarity.

Optical telescopes
Using fluorescence telescopes, the distribution of X max can be measured directly to obtain a handle on the value of the cross section at the highest energies [46,47].
The direct observation of the position of the shower maximum allows us to simplify Eq. ( 3) by removing the term P ∆X 2 describing the shower development after the shower maximum and the distribution P µ has to be replaced by a corresponding energy estimator based on the longitudinal shower profile.The resulting distribution of X max can be written as Here E em denotes the electromagnetic, i.e. calorimetric energy of the shower that can be obtained from the integration of the energy deposit profile.There is a very good correlation between the primary energy and the calorimetric energy [48].This correlation and the weak energy dependence of the depth of shower maximum allow us to neglect P Eem = dp em /dE em in our numerical studies.

Analysis of cross section measurement methods
The aim of all methods to measure the cross section is the determination of the interaction length λ int from measured distributions that represent the l.h.s. of Eqs.(3,5).

"k-factor" techniques
The approximation of an exponential attenuation of the frequency of air showers after the penetration of large amounts of atmosphere is the basis of the k-factor method [7,46,47], which has been used to analyze both X max and (N e , N µ ) data.The exponential slope of the attenuation is typically denoted by Λ, which is then related to the hadronic interaction length by a so-called k-factor While this method has been applied to data of air shower arrays as well as telescope detectors [7][8][9][10][11][12]14], the definition of the k-factor for a ground array, k S , and for a telescope detector, k X , are not identical.Air shower fluctuations enter differently into k S and k X and the detector resolutions are also very different.This can be understood based on Eq. ( 3) and ( 5), which can both be approximated by exponential distributions  11), and its statistical uncertainty.
for large X max , respectively X obs .For the (N e , N µ ) method this corresponds to while for the X max -tail method it is dN sh This assumes that there are no significant non-exponential contributions to the tails of the distributions, which is not true in general as it will be shown in the following.Since the full distributions are described by Eq. ( 3) and ( 5) it is possible to calculate the tails based on the given convolution integrals.Each of the integrations will give a contribution, so for (N e , N µ ) one gets where λ int is the proton-air interaction length, k ∆X 1 the contribution from the integration of P ∆X 1 , k ∆X 2 the one from P ∆X 2 and k S/X res the part due to the detector resolution.The individual contributions of k ∆X 1 , k ∆X 2 and k S/X res are difficult to compute and not known in most of the existing analyses, except the one in Ref. [12] and even more complete in Ref. [14].
For a given experimental setup and analysis approach it is possible to estimate the values of k ∆X 1 , k ∆X 2 , k X res and k S res as well as their dependence on hadronic interaction models.This is demonstrated here for an artificial experimental setup as described below.For each of the hadronic interaction models QGSJet01c, QGSJetII.3and SIBYLL 2.1 a set of ∼ 400, 000 proton induced air showers is simulated with a primary energy distribution ∝E −3 between 18.7 < log 10 E/eV < 19.4 and a zenith angle distribution of dN ∝ cos θ d cos θ.Two histograms are produced and analyzed for every set of simulations: • Distribution of X max with 18.9 < log 10 E/eV < 19.15.A Gaussian detector resolution with σ(X max ) = 20 gcm −2 and an energy resolution of σ(E)/E = 0.1 is assumed, which corresponds to values reported by the Pierre Auger Collaboration [49].The resulting mean energy of this data sample is 10 19.0 eV.• Distribution of events with 1.26 < N e /10 9 < 3.98 and 0.28 < N µ /10 8 < 1.12 at a vertical atmospheric depth of X vert obs = 860 gcm −2 (corresponding to the Pierre Auger Observatory) versus sec θ.The N e cut selects showers past their maximum (N max e (10 19 eV) ∼ 6.3•10 9 ) as done, for example, in the pioneering Akeno analysis [7] ¶.The detection resolution of N e and N µ are assumed to be Log-Normal with a resolution of σ(log 10 N e )/ log 10 N e = 0.05 and σ(log 10 N µ )/ log 10 N µ = 0.1 as well as the zenith angle uncertainty of σ(θ) = 1.0 • .Due to significant modeldifferences in the prediction of muon numbers the selected data samples have slightly differing energies (from 10 18.96 eV for QGSJet01 up to 10 19.04 eV for SIBYLL).
The histograms are fitted with an exponential function by a log-likelihood fit, which allows to correctly include empty bins in the analysis.The fit range is chosen as follows: While the end of the fit, X fit end , is always 50 gcm −2 beyond the last non-zero entry in the histogram, the start of the fit, X fit start , is varied.The resulting dependence of the exponential slopes on the start of the fitting range is then parameterized by where the asymptotic slope Λ 0 for X fit start → ∞ is taken to be the fit range-independent value of the exponential slope.In Fig. 6 this procedure is displayed for histograms obtained with the QGSJetII interaction model.It is found that the choice of the fit range has a major impact on the outcome of the slope.This is consistent with the results of [39].For the (N e ,N µ )-method, generally no plateau is found and the resulting value for Λ 0 is highly unstable.Thus, it is hardly possible to define a meaningful slope of the tail for this case.The situation for X max is somewhat better.While a stable plateau is still difficult to identify, Λ 0 can be estimated reliably.With this analysis approach the resulting asymptotic values, Λ 0 , are slightly lower than the values obtained for Λ directly from fits to the distributions.So neither the slopes nor the k-factors can be directly compared to results from previous analyses.Nevertheless, the found model dependence and other methodical problems will apply in a similar way to previous analyses.The strong dependence of the slope on the chosen fitting range is a severe methodical problem of the k-factor technique.Without a given method to infer a meaningful slope from such distributions, as for example the one proposed here, the k-factor technique cannot produce reliable results.The obtained slopes Λ 0 are then compared to the interaction mean free path length of the primary particle in the atmosphere, λ p−air , of the interaction model used in the CONEX simulations.In Tab. 1 the resulting k-factors are listed together with their propagated statistical uncertainties.The total model induced uncertainties on the k S/Xfactors are ∼ 7 % for the X max -tail and ∼ 28 % for the (N e , N µ )-method.Looking to the k-factor components one can see that a part of the model induced uncertainty enters already in the shower development up to the shower maximum (k ∆X 1 ), while the shower development after the shower maximum (k ∆X 2 ) adds the more significant amount of model-dependence.Another large contribution comes from k S res , which is mostly related to very differing model predictions on the number of muons.The factor k X res is only marginally model-dependent, which is originating from the very small model differences in the energy dependence of X max .
While this study mostly serves illustrative purposes, and the chosen experimental and analysis parameters are arbitrary, it nevertheless demonstrates the impact of model dependence on k-factor techniques.The found model dependence is rooted in the underlying air shower physics as well as in typical detector characteristics.

Unfolding of the X max -distribution
An improvement of the cross section measurement techniques is achieved by explicitly accounting for air shower fluctuations [13].This allows one to use not only the slope but also the shape of the X max -distribution.The measured X max -distribution, Eq. ( 5), is unfolded using a given P ∆X 1 -distribution to retrieve the original X 1 -distribution.Of course, the ∆X 1 -distribution needs to be inferred from simulations, which ultimately introduces a comparable model dependence as in the k-factor techniques [50].The model dependence of the most important part of the kernel function, P ∆X 1 , can be seen in Fig. 3.
In the unfolding technique, a larger range of the X max -distribution is used.If the primary particles are known to be protons only, a cross section analysis can be done already with very limited shower statistics.On the other hand, the results of this method are more sensitive to a possible fraction of primary particles that are not protons.

An improved method to derive the proton-air cross section
One of the shortcomings of cross section analysis methods applied so far is the missing connection between the cross section of the first interaction, that is to be measured, and the cross sections used in the calculation of the various probability distributions in Eqs. ( 3) and ( 5).
In the following we consider the method of unfolding the X max distribution, applicable to fluorescence telescope measurements, and improve it by consistently accounting for the "measured" high-energy cross section.We will calculate the full shape of the distribution of X max , depending only on properties of hadronic interactions at ultra-high energies.The impact of changing features of hadronic interactions at ultrahigh energies on the longitudinal air shower development, i.e.P ∆X 1 , is parameterized and used within the calculation.It is straightforward to refine the modeling by incorporating detector acceptance as well as the energy distribution of analyzed data.

Description of the X max -distribution
The description of the distribution of observed X rec max in terms of σ p−air is based on Eq. ( 5), now written with the term for the experimental X max resolution The parameter X shift allows us to shift the P ∆X 1 -distribution by replacing ∆X 1 with ∆X 1 + X shift .The introduction of the X shift parameter is necessary in order to reduce the model dependence of the analysis [50].This can be easily understood by looking at the P ∆X 1 -distributions shown in Fig. 3, which appear to be shifted between different models by up to ∼ 60 gcm −2 .It is found that the X shift parameter is highly model dependent and comprises many additional effects of the characteristics of high energy hadronic interactions on the X max -distribution.It can be demonstrated that differences between the inelasticity or the secondary multiplicity of the models may act as a cause for such an additional shift [51].The only important assumption related to the role of X shift in the cross section analysis is that any additional and unknown changing characteristics of hadronic interactions at extreme energies contributes mainly to a global shift of X max and thus ∆X 1 , while leaving the shape of the distributions unchanged.
The cross section analysis proceeds then as follows.The X rec max -distribution is calculated from Eq. ( 12) for a given interaction length λ int and shift parameter X shift and compared to the measured distribution.By performing a log-likelihood fit with the interaction length and an overall shift in depth as free parameters, the cross section and the 1σ uncertainty band is found [50].
One major improvement with respect to previous cross section analysis approaches is the consideration of the impact of a changing cross section on the resulting shower development described by P ∆X 1 .It is assumed that the cross section is reasonably well known at the energy corresponding to that of the Tevatron collider.Starting from this energy of ∼ 2 × 10 15 eV, cross sections of the model used for the calculation of P ∆X 1 are rescaled by a factor f (E) that increases logarithmically with energy and ensures that the modified model cross section matches the corresponding value of the fit parameter λ int at the considered shower energy.Here we will use 10 19 eV as reference for the scaling parameter being f 19 at this energy.Then the scaling factor reads ln(E/10 15 eV) ln (10 19 eV/10 15 eV) , for E > 10 15 eV and f (E) = 1 otherwise.The idea of the rescaling of the cross section is shown in Fig. 7, where the proton-air cross section of SIBYLL is shown together with cross sections that are scaled up and down by 20 % at 10 19 eV.So far the dependence of P ∆X 1 on the cross section has been neglected in air shower based cross section measurement.In fact, any attempt to determine the proton-air cross section without taking this dependence into account is overestimating the impact of σ p−air on the analyzed observables since part of the measured effect must be attributed to a modified development of the air shower and not to the fluctuations of the first interaction point.
To implement the idea of modified cross sections in the data analysis one needs to parameterize the ∆X 1 -distribution and determine its energy dependence and the modification of this distribution as function of the cross section scaling parameter.This is done using a Moyal distribution extended by one parameter and described in detail in Appendix A. In Fig. 8

Fitting range and stability
To investigate the sensitivity of the method to the choice of the X max -range used for the fit, we generate sets of 3000 simulated data showers and analyzes them with Eq. ( 12).The Gaussian detector resolution for X max is chosen to be 20 gcm −2 , which is the value reported by the Pierre Auger Collaboration [49].
The interval of the X max -distribution used to fit the model is important for the resulting statistical as well as systematical uncertainties.The first aspect is obvious since a reduction of the dataset is clearly leading to a reduced statistical power of the reconstruction.The latter one is mostly because of the possible contamination of the X max -distribution by cosmic ray primaries other than protons.All primary nuclei heavier than protons are producing shallower X max compared to protons.Primary photons, on the other hand, are deeply penetrating and have a larger X max than protons.A restricted fitting range in X max can thus be used to enrich the considered fraction of protons and reduce a possible contamination by other cosmic ray primaries.
The resulting impact of the chosen fitting range on the performance of the reconstruction is shown in Figure 9.The position of the peak of the X max -distribution, X peak , is used as a reference to define the fitting range.The starting point as well as the ending point of the fit are expressed only relative to X peak .
Evidently, it is beneficial to chose the fitting range as large as possible to get the smallest resulting statistical uncertainties.It is possible to shift the beginning of the fitting range relatively close to X peak without loosing too much statistical power (cf. Figure 9 (a)).It is found that, with the used parameterizations, the beginning of the fitting range can be set to 50 gcm −2 in front of X peak .In the following this is the adopted default choice.
The choice of the end of the fitting range has a similar impact on the reconstruction.In Figure 9 (b) it is shown how the reconstruction is degrading, while choosing a shorter (b) End of the fitting range, X stop max (relative to the peak of the X max -distribution) Figure 9. Impact of the chosen fitting range in X rec max on the resulting cross section (left) and X shift (right) at 10 19 eV.Each point corresponds to the mean of 100 reconstructions and the error bars denote the resulting RMS.Both, the reconstructed X max as well as the ∆X 1 -distribution are produced independently with the specified interaction model.The lines are just to guide the eye.fitting range.Since the photon fraction of ultra-high energy cosmic ray primaries is already strongly constrained [52] there is no special need to restrict the upper end of the fitting range.In the following the upper end of the fitting range for the log-likelihood fit of the X max -distribution is set to the value of the maximum X max plus 40 gcm −2 .
The cross section dependent parameterizations of the ∆X 1 -distributions are introducing a slight systematic overestimation of the reconstructed cross sections of 10 − 20 mb, corresponding to < 5%.This can be taken into account within the determination of the total systematic uncertainties of the measurement.The statistical resolution of the reconstruction is around 20 mb for 3000 events.

Comparison of the performance of different analysis methods
In order to test the ability of a reconstruction method to recover a changing input cross section, air shower simulations with a modified cross section (cf.Appendix A.2) are performed.The simulated air shower data are then reconstructed and the found cross  sections σ rec are compared to the modified input cross sections σ modified = f 19 σ model .The result of this analysis for all cross section reconstruction methods is summarized in Fig. 10.The central shaded area together with the solid line demonstrates the principal self-consistency of the analysis technique.The outer shaded area can be interpreted as the maximum possible model-dependence of the reconstruction method.
With the (N e , N µ )-method it is very difficult to reconstruct a modified cross section.Additionally the model dependence of the result is generally large.It is difficult to quantify, but can easily be of the order of hundreds of mb.
A very much better performance is achieved by the X max -tail analysis.The deviation of the reconstruction cross section from the input value is not getting larger than 50 mb and the model-dependence is between 100 − 200 mb.
The X max -unfolding technique can find cross sections that are not deviating much from the original model value.However, it shows a clear systematic trend of underestimating small cross sections and strongly overestimating large cross sections.Also the model dependence is not negligible, ranging from ∼ 75 mb at very low cross sections, ∼ 200 mb at intermediate cross sections to many hundreds of mb at large cross sections.
Finally, the X max -model demonstrates its unique ability to retrieve very consistently modified cross sections over a wide range.The model dependence is smaller than for the X max -tail method and is < 50 mb at low cross section while it grows to ∼ 100 mb at large cross sections.
The generally better sensitivity to reconstruct smaller compared to larger cross section is resulting from the increasing importance of the fluctuations of X 1 for the final X max -distribution.This facilitates the measurement of a small cross section, while for large cross sections X max -distributions are mostly shaped by fluctuations of the shower development P ∆X 1 and the detector resolution, making a measurement of the cross section generally more difficult.
Due to the combined impact of the fluctuation of both electron as well as muon numbers, the remaining sensitivity of the (N e , N µ )-method to the proton-air cross section is very much reduced compared to methods using the observable X max , see Fig. 11.
It is not straightforward to quantify the systematic effect caused by interaction models on the results of the proton-air cross section analyses.While the magnitude of the outer uncertainty band in Fig. 10 suggests relatively large model dependence of the results, this just reflects the existing incompatibility of the underlying models.Thus, it cannot be used to make any conclusion on the comparison of real data to the different models.For the analysis of real data, model dependence should be evaluated by the independent reconstruction of the data based on all available interaction models, and subsequent quoting of the mean result.The interval spanned by the reconstructions can then be used as an estimation of the systematic uncertainty induced by the models.This results in a smaller dependence on interaction models compared to Fig. 10.

Summary
It is demonstrated in methodical studies that all techniques to determine the proton-air cross section from air shower data are in fact based on the same fundamental formulation that describes the longitudinal air shower development.The requirement of using air shower Monte Carlo simulations to interpret EAS observables inevitably introduces a dependence on the hadronic interaction models used for the simulations.
It is found that the magnitude of the model dependence is similar for all cross section reconstruction methods; Only the (N e , N µ )-technique exhibits a significantly larger model dependence, which is due to the additional strong model dependence on the prediction of muon numbers.In addition, since the prediction of electron and muon numbers are both depending on interaction models, it is not possible to define model independent cuts for the selection of air shower events.Thus, already on the level of event selection, a model dependence is introduced in (N e , N µ )-analyses.
All k-factor based techniques are suffering from non-exponential contributions to the slope of the X max -respectively (N e , N µ ) frequency-distribution at large depths.Thus the exponential approximation used to define k-factors, is only of limited accuracy.Any non-exponential contribution creates a strong dependence of the exponential slope on the chosen fitting range (see also Ref. [53]).
Generally, k-factors are depending on the resolution of the experiment and can therefore not be simply transferred from one experiment to another, as done in some recent analysis [54].In particular, k X -factors are inherently different from k S -factors and can therefore not be transferred from an X max -tail analysis to that of ground based frequency attenuation or vice versa.
An X max -based analysis technique, as proposed in Ref. [13], was further improved to reduce methodical and model induced systematic uncertainties on the cross section result as much as possible.The presented improved method is a model of the X maxdistribution based on simple considerations on longitudinal air shower development combined with the parameterized impact of a changing cross section on the resulting air shower development.Previous reconstruction techniques assume a static ∆X 1distribution or k-factor, independent of the modified cross section.As, in general, the reconstructed cross section deviates from the predicted value of the hadronic interaction model used to generate the ∆X 1 -distribution or k-factor, a discrepancy of the data with the air shower simulations is inevitable.The data can not be described by the models used for the air shower simulations and it is not sufficient to adopt the measured cross section just for the first interaction point, while keeping the rest of the shower simulations unchanged.
For the case of a proton-dominated cosmic ray data sample around 10 19 eV, which was studied in this article, it is possible to determine the proton-air cross section with good precision.The statistical uncertainties of the analysis of samples with 3000 events are about 20 mb (5 %) and the systematic uncertainties introduced by the parameterizations of the X max -model are of the same order of magnitude.The modelinduced systematic uncertainty is between 50 and 100 mb, corresponding to relative uncertainties of < 10 % to ∼ 20 %.A further reduction of this model-related uncertainty is possible, but needs improvements in the understanding of the underlying hadronic interaction physics.Furthermore, the second free parameter of the new analysis technique, X shift , is a further handle that is sensitive to the physics of hadronic interactions at ultra-high energies.A non-zero result of X shift indicates the deficiency of the underlying hadronic interaction model to describe the measured data distribution.It can be interpreted in terms of a modified characteristics of secondary particle production in hadronic interactions [51].A2.Furthermore, it is convenient to introduce the relative change as Appendix A.3.Parameterization of the ∆X 1 -distribution Combining the found dependence of ∆X 1 -distributions on the energy E and the cross section modifier f 19 , it is now possible to construct a parameterization of P ∆X 1 in terms of E and σ(E).Firstly, the principle energy dependence of P ∆X 1 is evaluated by calculating α e (E), β e (E) and γ e (E) using Eq.(A.3) and the results from Tab. A1.Secondly, the effect of the changed cross section is added.To achieve this, the first step is to calculate the corresponding deviation of the cross section σ(E) extrapolated to 10 19 eV assuming Eq. ( 13).The resulting factor f 19 is then used to evaluate ∆α σ (f 19 ), ∆β σ (f 19 ) and ∆γ σ (f 19 ) using Eq.(A.5) and the results listed in Tab.A2.In accordance to Eq. ( 13) the energy dependence is assumed to be logarithmic and vanishing below 10 Figure1.Compilation of proton-air production cross sections from cosmic ray measurements[2][3][4][5][6][7][8][9][10][11][12][13][14].The data are compared to model predictions[16][17][18][19][20][21].

Figure 2 .
Figure 2.Definition of variables to characterize EAS longitudinal shower development.Zero slant depth is where the cosmic ray particle enters the atmosphere.The first interaction occurs at X 1 .At X max the shower reaches its maximum particle number N max .After the maximum the shower is attenuated over ∆X 2 before it reaches the particle number N stage at the slant depth of X obs .

1 Figure 3 .
Figure 3. Distribution of the depth distance from the point of the first interaction up to the shower maximum for proton primaries at 10 19 eV.The predictions of different models are compared.

Figure 4 .
Figure 4. Independence of P ∆X1 of X 1 , demonstrated for primary protons at 10 19 eV with QGSJetII.

Figure 5 .
Figure 5. Model dependence of longitudinal air shower development for proton primaries.Left panel: Slant depth distance from the shower maximum up to the depth where showers are attenuated to N e (E e > 1 MeV) = 10 9 electrons.Middle panel: Number of muons above 1 GeV after 1000 gcm −2 of shower development.Right panel: Same as left panel, but for showers selected according to 7.6 < log 10 N µ (E µ > 1 GeV) < 7.7.While the plots on the left and in the middle are obtained from fixed energy simulations, the right plot is based on an energy spectrum ∝ E −3 .

Figure 6 .
Figure 6.Left panel: (N e , N µ )-method.Right panel: X max -tail method.The fitted open symbols show the dependence of the resulting exponential slope on the choice of the begin of the fit range.The arrow indicates the asymptotic behavior of the fit, deduced from Eq. (11), and its statistical uncertainty.

Figure 7 .
Figure 7. Modified extrapolation of the proton-air cross section.Shown are the original prediction from SIBYLL together with a 20 % increased respectively decreased cross section at 10 19 eV.
example ∆X 1 -distributions are shown together with fits of the extended Moyal distribution and also the final parameterizations.

1 Figure 8 .
Figure 8. Simulated ∆X 1 -distributions at 10 19 eV for two hadronic interaction models.Included are the fits of the extended Moyal distribution Eq.(A.1) as well as the final parameterization.
Start of the fitting range, X start max (relative to the peak of the X max -distribution)

Figure 10 .
Figure 10.Resulting systematics and dependence on models of cross section reconstruction techniques.Every analysis was performed 100 times on 3000 simulated EAS; Shown are the mean results.The central shaded area shows the ability of a method to reconstruct the input cross section of simulated air shower data with the same model that was used for the simulations.The thick line marks the mean result for all used models, which are QGSJet01c, QGSJetII.3and SIBYLL2.1.The width of the central band is caused by statistical fluctuations.The outer shaded area denotes the maximum deviation of the cross section reconstructed with a model that differs from the one used to generated the simulated dataset.

Figure 11 .
Figure 11.Shower-to-shower fluctuations of X max (left) and log 10 N e (right) for primary protons.The solid line is the original model prediction of QGSJetII, while the dashed line is for showers having their first interaction at the fixed slant depth corresponding to the predicted proton-air interaction length of QGSJetII.The dotted line is for simulated air showers with increased cross section by a factor of two, f 19 = 2, in addition to the fixed first interaction point.

Figure A2 .
Figure A2.Parameterization of the f 19 -dependence of the ∆X 1 -parameterizations with respect to a changing cross section.The lines denote the interpolations as discussed in the text.

Table 1 .
Resulting k-factors from CONEX study for primary protons at 10 19 eV.In each column the maximal and minimal k-factors are highlighted.The difference between those maxima is denoted by ∆ and is a measure of model-dependence.