{'fraction_non_alphanumeric': 0.03727714748784441, 'fraction_numerical': 0.01469475958941113, 'mean_word_length': 4.172115115954177, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On the influence of statistics on the determination of the mean value of the depth of shower maximum for ultra high energy cosmic ray showers

The chemical composition of ultra high energy cosmic rays is still uncertain. The latest results obtained by the Pierre Auger Observatory and the HiRes Collaboration, concerning the measurement of the mean value and the fluctuations of the atmospheric depth at which the showers reach the maximum development, Xmax, are inconsistent. From comparison with air shower simulations it can be seen that, while the Auger data may be interpreted as a gradual transition to heavy nuclei for energies larger than ~ 2-3x10^18 eV, the HiRes data are consistent with a composition dominated by protons. In Ref. [1] it is suggested that a possible explanation of the observed deviation of the mean value of Xmax from the proton expectation, observed by Auger, could originate in a statistical bias arising from the approximated exponential shape of the Xmax distribution, combined with the decrease of the number of events as a function of primary energy. In this paper we consider a better description of the Xmax distribution and show that the possible bias in the Auger data is at least one order of magnitude smaller than the one obtained when assuming an exponential distribution. Therefore, we conclude that the deviation of the Auger data from the proton expectation is unlikely explained by such statistical effect.


Introduction
The nature of the primary cosmic rays is intimately related to the astrophysical objects capable of accelerating these particles to such high energies. Also, propagation in the intergalactic medium depends on the composition, which affects the resulting spectral distribution of the flux observed at Earth. A knowledge of the composition is also very important for primary energy reconstruction and for anisotropy studies.
One of the most important limitations of composition analyses comes from the lack of knowledge of the hadronic interactions at the highest energies. Composition studies are based on the comparison of experimental data with Monte Carlo simulations of atmospheric cosmic rays showers, which makes use of hadronic interaction models which extrapolate the available low energy accelerator data to the energies of the cosmic rays.
One of the most sensitive parameters to the mass of the primary cosmic ray is the atmospheric depth at which the showers reach their maximum development. Lighter primaries generate showers that are more penetrating, producing larger values of X max . Also, the fluctuations of this parameter are smaller for heavier nuclei. The Pierre Auger Observatory and the HiRes experiment are able to observe directly the longitudinal development of the showers by means of fluorescence telescopes. Therefore, in both experiments, the X max parameter of each observed shower can be reconstructed from the data taken by the telescopes.
The mean value and the standard deviation of X max , as a function of primary energy, obtained by Auger [2] and HiRes [3] appear to be inconsistent. From the comparison with simulations, the Auger data suggest a transition to heavier nuclei starting at energies of order of 2 − 3 × 10 18 eV, whereas, the HiRes data are consistent with protons in the same energy range. In Ref. [1] a new parameter, the difference between the mean value and the standard deviation of X max , was introduced in order to reconcile the Auger and HiRes results. This new parameter has the advantage of being much less sensitive to the first interaction point than the mean value and the standard deviation separately. From a comparison of the experimental values of this parameter, obtained by Auger and HiRes, with simulated data, they infer that the composition of the cosmic rays is dominated by protons. They say that the energy dependence of the distribution of X max , observed by Auger, seems to be caused by an unexpected change in the depth of the first interaction point, which can be explained by a rapid increase of the cross section and/or increase of the inelasticity. Both possibilities require an abrupt onset of new physics in this energy range, which makes them questionable. They also suggest that the deviation of the distribution of X max from the proton expectation, present in the Auger data, could be originated in the statistical techniques used to analyze the data. In particular, they suggest that the deviation of the mean value of X max from the proton expectation could be explained by a bias originated from the exponential nature of the X max distribution and the decreasing number of events as a function of primary energy.
In this work we show that, considering a better description of the X max distribution, the bias in the determination of the mean value of X max become more than one order of magnitude smaller than the one obtained for the exponential distribution. We find that the value of the bias in the last energy bin (the one with the smallest number of events) of the Auger data, published in Ref. [2], is 1.5 g cm −2 , which is much smaller than the systematic errors on the determination of the mean value of X max estimated in Ref. [2].

Numerical approach
Following Ref. [1] let us introduce the parameter, where X max is the mean value of the X max distribution, is the sample mean corresponding to samples of size N and mode[X N max ] is the value ofX N max that occurs most frequently, i.e. the maximum of the distribution function of X N max . Therefore, the bias on the determination of X max appears when a particular realization of the sample mean is equal to the mode of the sample mean distribution function. Note that the sample mean (Eq. (2)) is an unbiased estimator of the mean of the exponential distribution, i.e. E[X N max ] = X max . In Ref. [1] it is shown that approximating the X max distribution by an Exponential function the parameter ξ(N) is given by: ξ E (N) = 1/N.
In order to better describe the distribution of X max two different types of functions are considered. They are chosen in such a way that the distribution ofX N max can be obtained, at least, in a semi-analytical way. The first function considered is a shifted-Gamma distribution [4], where k = 5 and the other two parameters can be obtained from the mean value and the standard deviation of X max , The second function under consideration is the convolution between an exponential function and a Gaussian (Exp-Gauss), where α, β and λ are fitting parameters and A library of simulated showers was generated by using the program CONEX (v2r2.3) [5]. Monochromatic samples of 10 4 proton showers were generated from log(E/eV ) = 18 to log(E/eV ) = 19.5 in steps of ∆ log(E/eV ) = 0.1. The arrival directions of the showers follow an isotropic distribution, such that the zenith angle is in the interval [0 • , 60 • ]. The hadronic interaction models considered are QGSJET-II [6] and EPOS 1.99 [7].
The mean value and the standard deviation (needed for the description of the X max distribution using the shifted-Gamma function) were fitted with a quadratic function and a linear function of log(E), respectively, i.e., Figure 1 shows the simulated data as well as the fits, for both hadronic interaction models considered. The values of the parameters corresponding to Eqs. (8) and (9)   given in table 1. Table 1. Parameters corresponding to the quadratic and linear fits of X max and σ[X max ], respectively (see Eqs. (8) and (9)), obtained from simulations for QGSJET-II and EPOS 1.99. The distribution functions of X max , for every energy and hadronic interaction model considered, were fitted with the Exp-Gauss function, Eq. (6). The parameters α, β and λ were fitted with linear functions of log(E), in order to obtain the Exp-Gauss representation of the X max distribution for every value of energy in the interval [10 18 , 10 19.5 ] eV, see Appendix A for details. Figure 2 shows the distributions of X max , obtained by using CONEX with QJSJET-II, for log(E/eV ) = 19 and log(E/eV ) = 19.5. Red solid lines correspond to the fits of the simulated data with the Exp-Gauss function. The blue dashed lines correspond to the shifted-Gamma function, Eq. (3), for which the parameters X 0 and τ X are obtained by using the expressions of X max and σ[X max ] in Eqs. (8) and (9) to calculate X 0 and τ X from Eqs. (4) and (5), respectively. From the figure it can be seen that the Exp-Gauss function is a better fit to the simulated data than the shifted-Gamma function. It can also be seen that the tail to larger values of X max is slightly overestimated by the Exp-Gauss distribution and underestimated by the shifted-Gamma function. Therefore, the distribution function of the universe (samples with N → ∞) should fall between this two functions. The distribution ofX N max can be calculated by means of the characteristic function, which is defined as the expectation value of exp(itX max ), i.e. φ Xmax (t) = E[exp(itX max )]. It is straightforward to show that the characteristic function ofX N max is given by The characteristic function of the shifted-Gamma distribution is φ G Xmax (t) = exp(iX 0 t) (1 − itτ X ) −k and then the characteristic function ofX N max is given by which corresponds also to a shifted-Gamma distribution. Therefore, the distribution function ofX N max is given by, By using Eq. (10) it is easy to show that, In this case, ξ is also proportional to 1/N but it is suppressed by the ratio between the standard deviation and the mean value of X max . A similar expression is obtained when the distribution function of X max is described by a truncated exponential function, see Appendix B for details. The blue solid line on the left panel of Fig. 3 corresponds to ξ G as a function of N for log(E/eV ) = 19.5, approximately the mean value of the energy (weighted by the spectrum) for the last bin considered in Ref. [2]. Note that, the number of events in this bin is 34. From the figure, it can be seen that ξ G is more than one order of magnitude smaller than the function 1/N. The distribution function of X max is affected by the presence of fluctuations introduced by the detectors. The distribution function ofX N max , including a Gaussian uncertainty on the determination of X max is given by, where σ Rec is the standard deviation of such uncertainty. The mode of this distribution is calculated numerically. Dashed and dashed-dotted lines on the left panel of Fig.  3 correspond to parameter ξ G (N) obtained for σ Rec = 20 g cm −2 and σ Rec = 40 g cm −2 , respectively. When a symmetric uncertainty on the determination of X max is included, the parameter ξ becomes still smaller and decreases for increasing values of the uncertainty. This is due to the fact that ξ is larger for asymmetric distributions, like the exponential, and the convolution of the pure X max distribution with a Gaussian is more symmetric than the original one.
The characteristic function of the Exp-Gauss distribution is the product of the characteristic function of the exponential distribution, (1 − iλt) −1 , with the one corresponding to a Gaussian, exp(iαt − β 2 t 2 /2). Then, the characteristic function of X N max is then given by, which corresponds to the convolution of a Gamma distribution with a Gaussian, Last integral is calculated numerically in order to obtain the mode of the resultant distribution. The solid red line in the right panel of Fig. 3 shows ξ EG as a function of the sample size for log(E/eV ) = 19.5. Note that ξ G is smaller than ξ EG , this is due to the more extended tail to larger values of the Exp-Gauss distribution compared with the corresponding one to the shifted-Gamma distribution. In any case, ξ EG is still about one order of magnitude smaller than 1/N. As for the case of the Gamma distribution, dashed and dashed-dotted red lines correspond to σ Rec = 20 g cm −2 and σ Rec = 40 g cm −2 , respectively. In this case the effect of the uncertainty on the determination of X max is included inP EG just by replacing the parameter β byβ = β 2 + σ 2 Rec . As expected, the curves that include the uncertainty on the determination of X max fall bellow the one corresponding to the ideal case. The left panel of Fig. 4 shows the parameter ξ as a function of energy corresponding to the number of events in each energy bin taken from Ref. [2], for the case in which there is no uncertainty on the determination of X max (which gives larger value of ξ, as shown before). The energy assigned to the ith bin, used to calculate ξ, corresponds to the mean value of the energy in the bin weighted by the broken power law fit of the cosmic rays energy spectrum, J(E), of Ref. [9], where E L i and E U i are the lower and upper limits of the ith bin. It can be seen that the values of ξ, obtained by using the Exp-Gauss distribution and the shifted-Gamma distribution, are more than one order of magnitude smaller than the corresponding one for the exponential distribution, in the whole energy range and for both hadronic interaction models considered. As in the previous calculation, ξ G results are smaller than ξ EG . In fact, the ξ curve corresponding to the true distribution of X max should fall between the curves corresponding to the Exp-Gauss and the shifted-Gamma representation of the X max distribution.  Figure 4. ξ (left panel) and ∆X max (right panel) as a function of log(E/eV ), for the statistics of the Auger data of Ref. [2]. Solid lines correspond to QGSJET-II and dashed lines correspond EPOS 1.99.
The right panel of Fig. 4 shows the parameter ∆X max = X max ξ which gives the grammage of the shift suffered by X max ifX N max takes the value of the mode of its distribution. It can be seen, that for the last energy bin, the one with 34 events, ∆X max is 1.5 g cm −2 , which is much smaller than the systematic uncertainties on the determination of X max estimated in Ref. [2].
The energy bins considered in the analysis of Ref. [2] have a width of ∆ log(E/eV) = 0.1 in the energy range from E = 10 18 eV to E = 10 19 eV. Between E = 10 19 eV and E = 10 19.4 eV, ∆ log(E/eV) changes to 0.2 and the last bin corresponds to E ≥ 10 19.4 eV. Therefore, the number of events per bin decreases in the energy range from E = 10 18 eV to E = 10 19 eV, it increases from 96 in the bin [10 18.9 , 10 19 ] eV to 138 in the bin [10 19 , 10 19.2 ] eV and then, it decreases for the last two bins. This change in the bin width generates the structure around E ∼ = 10 19.1 eV seen on the curves of Fig. 4.
Note that ξ EG calculated by using EPOS 1.99 is larger than the corresponding one for QGSJET-II, this is due to the fact that the X max distributions obtained with EPOS 1.99 are more asymmetric (increase faster, coming from small values of X max , and have a more extended tail) than the corresponding ones to QGSJET-II.
Concerning iron showers, it can be seen that ξ takes smaller values than the ones for protons. This is due to the large suppression of fluctuations in iron showers, the ratio of the standard deviation to the mean value of X max is smaller than for protons, producing smaller values of ξ (see Eq. (11)). In particular, ξ f e G = ξ pr G /K where K increases from ∼ 2.3 at E = 10 18 eV to ∼ 2.4 at E = 10 19.5 eV for QGSJET-II.

Conclusions
In this work we studied in detail statistical bias in the determination of the mean value of X max , suggested in Ref. [1], as a possible explanation of the deviation of Auger data from the proton expectation. We used two different functions to fit the X max distribution obtained from simulations: (i) the convolution of an Exponential distribution with a Gaussian and (ii) a shifted-Gamma distribution. We find that the bias obtained by using these two functions is more than one order of magnitude smaller than the corresponding one of the Exponential distribution, the one used in Ref. [1]. We find that the values of the bias, obtained for the convolution of the Exponential function with the Gaussian, are larger because it presents a more extended tail to larger values of X max than the shifted-Gamma distribution. We also find that the bias diminishes when a Gaussian (symmetric) uncertainty on the determination of X max is included.
We also calculated the expected bias, as a function of primary energy, using the actual number of events in each energy bin of the Auger data, published in Ref. [2], for both hadronic interaction models considered in this work, QGSJET-II and EPOS 1.99. We find that the largest value of the bias, corresponding to the bin with the smallest number of events, is smaller than 1.5 g cm −2 , much less than the systematic errors on the determination of X max estimated in Ref. [2].

Appendix A. Parameters for the Exp-Gauss fits
The parameters α, β and λ, obtained from the fits of the X max distributions with the Exp-Gauss function (see Eq. (6)), are fitted with linear functions of log(E/eV ) as shown in figure A1. They can be written in the following way, where the coefficients C i , i = 1...6, are given in  Figure A1. Parameters α, β and λ corresponding to the fits of the X max distribution with the Exp-Gauss function for QGSJET-II and EPOS 1.99. The straight lines correspond to the linear fits of the points. It is possible to describe the X max distribution function by a truncated exponential distribution, which is given by, where Λ is a parameter that describe the tail of the X max distribution and X c is the truncation value. The characteristic function of this distribution is, φ T E Xmax (t) = exp(itX c ) (1 − itΛ) −1 and then, the characteristic function of the sample mean is given by, φ T Ē X N max (t) = exp(itX c ) (1 − itΛ/N) −N , which corresponds to a shifted Gamma distribution. Therefore, the distribution function of the sample mean is given by, By using that X max = X c + Λ, it is easy to show that, Note that, it can be seen, form Eq. (B.4), that ξ T E takes a very similar form to the one obtained for the shifted-Gamma function, see Eq. (11). Typical values of the parameters, obtained experimentally, are X c ∼ = 700 g cm −2 and Λ ∼ = 56 g cm −2 [10] (note that these parameters depend on primary energy and the ones used here, obtained from Ref. [10], correspond to the energy interval [10 18 , 10 18.5 ] eV, in any case, they are just used to roughly estimate the suppression factor of the bias). Therefore, ξ T E (N) ∼ = 0.125/N, which is suppressed by a factor 0.125 with respect to the corresponding one to the exponential distribution.