{'fraction_non_alphanumeric': 0.03444795279258761, 'fraction_numerical': 0.018116879755680937, 'mean_word_length': 4.309743025972241, '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': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A New Numerical Technique to Determine Primary Cosmic Ray Composition in the Ankle Region

In this paper we introduce a new multiparametric technique that attempts to tackle simultaneously the problems of composition determination and hadronic interaction uncertainty. Employing simulations of a real world detector under its planned operational conditions, and disregarding systematics, we can asses that the present technique should be able to determine the composition of a binary mixture of p and Fe with a statistical confidence of few percent, in a way that is independent of the assumed hadronic interaction model. Moreover, the combination of real data with the tools developed and presented here should give an indication of the reliability of the various hadronic interaction models in current use in the area. We center our study in the region of the ankle, where composition carries critical astrophysical information, and use two main parameters: the number of muons at 600 m from the shower axis and the depth of the shower maximum obtained from the hybrid operation of the planned muon counters and high elevation fluorescence telescopes of the AMIGA and HEAT Auger enhancements.


Introduction
The cosmic ray energy spectrum, in the high energy region above 10 17 eV, presents two main features observed by several experiments, the second knee, observed at 4 × 10 17 eV [1,2,3,4] and the ankle. There is evidence of a fourth feature situated at the highest observed energies, consistent with the so-called GZK suppression [5,6], which would be caused by the interaction of the ultrahigh energy protons with the photons of the cosmic microwave background radiation (CMBR) [7,8]. For the case of heavier nuclei a similar effect is expected because of their interaction with photons from the infrared and microwave backgrounds [9].
The origin of the second knee is still unclear. It has been interpreted as the end of the efficiency of the acceleration in Galactic supernova remnant shock waves [10], a change in the diffusion regime in our galaxy [11,12] or even the transition between the Galactic and extragalactic components of the cosmic rays [13].
The ankle is a broader feature. It has been observed by Fly's Eye [2], Haverah Park [14], Yakutsk [3], HiRes [4] and Auger [5] in Hybrid mode at approximately the same energy, ∼ 3×10 18 eV. The origin of the ankle is also unknown. It can be interpreted as the transition between the Galactic and extragalactic components [15] or the result of pair production by extragalactic protons after the interaction with photons of the CMBR during propagation [16,17].
There are three main models of the Galactic-extragalactic transition. The first is the mixed composition model [15], in which extragalactic sources inject a spectrum of masses similar to lower energy Galactic cosmic rays and for which the transition takes place at the ankle. The second is the ankle model [18], a two-component transition from Galactic iron nuclei to extragalactic protons at the ankle energy. The third is the dip model [13], in which the ankle is due to pair production of extragalactic protons that interact with the photons of the CMBR (in this scenario the transition occurs at the second knee).
In order to rule out or substantiate any of those models, additional information is necessary besides the energy spectrum shape and absolute intensity. Detailed measurements of the composition as a function of energy, while not sufficient, would be extremely valuable to break the present degeneracy among competing models for the Galactic-extragalactic transition [19,20]. Furthermore, this kind of information could help to determine what the highest energy accelerators in the Galaxy are and provide indicators of the kind and level of magnetohydrodynamic turbulence present in the intergalactic medium traversed by the lowest energy cosmic ray particles [10].
Several experiments have measured the cosmic ray composition in the region where the transition takes place. Nevertheless, large discrepancies exist between different experiments and experimental techniques [21]. One of the main reasons behind the plurality of sometimes contradicting results is that composition is determined by comparing experimental data with numerical shower simulations. These simulations include models for the relevant hadronic interactions which are extrapolations, over several orders of magnitude in center of mass, of accelerator data to cosmic ray energies. This is a source of considerable uncertainty which is confirmed, to a certain extent, by the fact that there is experimental evidence of a deficit in muon content of simulated showers with respect to real data [22].
In this paper we introduce a new statistical method to test the compatibility of the high energy hadronic interaction models and real data. In fact, the hadronic interaction model is assumed to be the main systematic uncertainty in the present analysis. Throughout our analysis, we consider QGSJET-II [23,24] and Sibyll 2.1 [25]. Our method allows us to verify whether the experimental data are compatible with the hadronic models under consideration and, if so, to estimate the composition.
Although this new technique is of general applicability, we study its potential in the context of AMIGA (Auger Muons and Infill for the Ground Array) [26] and HEAT (High Elevation Auger Telescopes) [27], the lower energy extensions of the Southern Pierre Auger Observatory. These two enhancements will extend the energy range down to 10 17 eV, encompassing the second knee and ankle region where the Galactic-extragalactic transition takes place.
The most relevant mass sensitive parameters that will be obtained from AMIGA and HEAT are the number of muons at 600 m from the shower core, N µ (600), and the atmospheric depth of maximum shower development, X max , respectively. Consequently, we develop our statistical technique by using mainly this pair of parameters. Combinations of N µ (600) with parameters from the Cherenkov detectors like the slope of the lateral distribution function, rise-time of the signals and curvature radius are also studied.
Note that our technique should also be applicable to Telescope Array and its low energy extension [28,29], which have hybrid capabilities in the region of the ankle and also plan to include muon detectors [29,30]. The quoted expected error for X max is ∼ 20 g cm −2 [30], which is comparable to the value estimated for HEAT.
for S(r 0 ), the lateral distribution function value at a fixed distance r 0 ∼ 600 m, as measured by the AMIGA surface array of water Cherenkov detectors. The same procedure should be, in principle, applicable to other detectors with hybrid capability like Telescope Array. This is so because N µ (600) and energy are correlated through S 600 and the latter also receives an important contribution from the electromagnetic lateral distribution function, while N µ (600) is directly used as a composition parameter. Therefore, at energies of few 10 18 eV or smaller, under the assumption of a binary mixture of proton and iron, the separation in N µ (600) due to composition alone at fixed energy is larger than the separation in S 600 due to composition at the same energy. This can be clearly seen, for example, in figures 12.a and 12.b of Ref. [31], where the merit factors (or discrimination power) for different composition indicators are shown as a function of energy for two zenith angles. The merit factor is, basically, the separation between the distribution functions of each parameter at a given energy normalized by the combined dispersion of both distributions. It can be seen that, at a fixed energy, the separation between the distributions of N µ (600) for proton and iron, the main composition indicator, is much larger than the separation between the corresponding distributions of S 600 , the energy estimator. Consequently, the fact that N µ (600) is correlated simultaneously to energy and composition does not inhibit its use as a composition parameter. Section 3.3 shows our main results. In particular, it is shown there that, by working on a plane defined by two parameters derived from X max and N µ data, an estimation of the cosmic ray composition can be obtained that is reasonably independent of the uncertainties related to the underlying hadronic interaction model. Furthermore, we demonstrate that rather small samples of events at a given reconstructed energy, compatible with the level of statistic expected from detectors currently under construction that will operate in the ankle region, are enough to this end. Additionally to the determination of composition, given two possible interaction models and an observed data set, the technique can be used to assess the compatibility of these models with the experimental data.

Abundance Estimator
In order to develop a statistical method to infer the composition of the cosmic rays (i.e. the abundance of a given primary type), let us consider two possible types of primaries, A = a, b, and samples of size N = N a + N b , where N a and N b are the number of events corresponding to type a and b, respectively. From each event of an individual sample it is possible to extract several observable parameters sensitive to the primary mass. Therefore, for a given mass sensitive parameter q we define, where q A i are N A random variables distributed as f A (q) and which is the probability that an event is of type a for a given q, assuming no prior knowledge of the primary type. Note that we restrict our analysis to the case in which the cosmic rays are the superposition of two components.
The new statistic, so defined, is an estimator of the abundance of the primary of type a. ξ q is a random variable because it is a function of N random variables, ξ q = ξ q (q a 1 . . . q a Na , q b 1 . . . q b N b ). Therefore, its mean value as a function of the composition is given by, where c a = N a /(N a +N b ) is the composition (or the abundance) corresponding to the primary a for the sample of N events. Noting that f A (q) are probability density functions and, therefore, that dq f a (q) ≡ 1, the dependence of ξ q on c a can be explicitly seen after integrating Eq. 3, where, Equation (4) shows that the mean value of ξ q increases linearly with the composition of the samples and from equations (5,6) we see that as the overlap between the distributions f a (q) and f b (q) decreases, ξ q (c a ) tends to the identity function. We also see that ξ q (1/2) = 1/2, independent of the shape of the distributions f a (q) and f b (q).
The variance of ξ q is given by, where, Therefore, the variance of ξ q also increases linearly with c a . Note that it is proportional to N −1 , i.e., it approaches zero as the size of the sample approaches infinity.
The composition estimator ξ q is defined for each parameter q, and therefore a covariance can be calculated for any two of such estimators, say q µ and q ν . Let f A (q µ , q ν ) be the distribution function of these variables and .
The covariance between ξ qµ and ξ qν is given by, From Eq. (13) we can see that Cov(ξ qµ , ξ qν ) is also a linear function of c a , and a sufficient condition to be zero is that the variables q A µ and q A ν are independent.
The parameter ξ q is the sum of N random variables, therefore, for large enough values of N, this variable follows a Gaussian distribution because of the central limit theorem. parameter η = | q 1 − q 2 | /[σ 2 (q 1 ) + σ 2 (q 2 )] 1/2 , which measures the discrimination power of q, can be easily calculated and takes the values η = x/ √ 2. Fig. 1 shows the corresponding slope and the intercept of ξ (c a ) as a function of η. We see that, as η increases, the intercept approaches zero and the slope approaches one, i.e., as expected, ξ (c a ) tends to the identity function.

Shower and detector simulations
In order to perform composition analyses around E = 1 EeV (EeV = 10 18 eV), including the effect of the energy uncertainty, we generated a library of atmospheric showers by using the AIRES package version 2.8.2 [32]. We used a relative thinning of 0.6 and statistical weight factor 0.2 (see Ref. [32] for details). For simplicity, and to a good approximation, we assume that the showers follow a power law energy spectrum with spectral index 2 γ = −2.7 in the energy interval [0.6, 2] EeV. We generated eight sets of 166 showers each corresponding to two kinds of primaries, protons and iron nuclei, two values of zenith angle, θ = 30 • and θ = 45 • , and two different models of high energy hadronic interactions, QGSJET-II and Sibyll 2.1.
Each of these showers was used to generate events of the AMIGA-HEAT detector system. The simulation of the AMIGA detectors was performed with a dedicated package described in Ref. [31]. Each shower was used 40 times by distributing impact points uniformly in the 750 m-array. We assumed 30 m 2 muon detectors segmented in 192 cells and buried 2.5 m underground [31]. We used a time binning of 20 ns and an efficiency of each segment equal to one. The shower arrival directions and core positions were reconstructed with the standard CDAS package Er-v3r4 [33] specifically developed to reconstruct the Cherenkov detector information in Auger. To reconstruct the muon lateral distribution function, we used the method introduced in Ref. [31].
The effects of the response of the HEAT telescopes and the reconstruction procedure for the longitudinal profile in the determination of X max were included using the approach of Ref. [31]. For each simulated AMIGA event, corresponding to a given shower, we obtained a reconstructed X max by taking a random value from a Gaussian distribution of mean value equal to the one calculated internally in AIRES (obtained by fitting the longitudinal profile with a Gaisser-Hillas function) and σ, for the corresponding energy, obtained from the interpolation of the simulated data given in Ref. [27]. Note that, despite the fact that the fluctuations in X max have been included in a realistic way, the same cannot be said of the possible correlations between the fluctuations in fluorescence and surface parameters for hybrid events. Nevertheless, we expect that such correlations could be dealt with by proper quality cuts without affecting our main conclusions.
In short, for each simulated event, we obtained all the parameters with the corresponding fluctuations of the Cherenkov, muon and fluorescence detectors. We obtained 8 sets of N 0 ∼ = 166 × 40 = 6640 events each (depending on the reconstruction efficiency). To refer to each set, we will use the notation S(θ, A, h) where θ = 30 • , 45 • , A =Proton, Iron and h =QGSJET-II, Sibyll 2.1.

Probability density function estimation
In order to calculate ξ (see Eq. (1)), we need the distribution functions, f A (q), of the different parameters sensitive to the primary mass considered, including the effects of the detectors and reconstruction methods. We use the nonparametric method of kernel superposition [34,35,36,37] as an estimate of these probability density functions from the simulated data.
We assume a Gaussian error of 25% in order to include the uncertainty in the determination of the primary energy, which has an important effect in the discrimination power of N µ (600). Therefore, we consider samples of events obtained in the following way: for each simulated event, belonging to a given set S(θ, A, h), of real energy E, we estimate the reconstructed energy by taking a random value from a Gaussian distribution of mean E and σ = 0.25 × E. Since the bias introduced by the rapid fall of the spectrum shifts the center of the bin from 1 EeV to 1.14 EeV (see Appendix B), if the "reconstructed" energy falls in an energy interval of width 0.5 EeV centered at 1.14 EeV, the event is added to a new sample, S Πr (θ, A, h). We repeat this procedure 10 times for each set S(θ, A, h). Therefore, we obtain 10 samples S i The parameters used in the analysis are defined as follows: (i) N µ (600) is the number of muons at 600 m from the shower core, which is estimated from a fit to the lateral distribution of the number of muons measured by the AMIGA muon counters [31], (ii) R is the curvature radius of the shower front, (iii) β is the slope of the lateral distribution function of the signal in the water Cherenkov detectors and (iv) t 1/2 is a parameter constructed with the rise-time of the signal in a selected subset of the triggered water Cherenkov detectors, where N T is the number of stations with signal greater than 10 VEM (Vertical Equivalent Muon 3 ), t i 10 and t i 50 are the times at which 10% and 50% of the total signal is collected, respectively, and r i is the distance from the i-th station to the shower axis. Only stations at a distance to the shower axis greater than 400 m are included in Eq. (14).
Better estimates of the density functions are obtained using the adaptive bandwidth method introduced by B. Silverman [34]. We perform a first estimation of each density function by using a Gaussian kernel with a fixed smoothing parameter, where x is a two-dimensional vector of one of the pairs of the parameters considered, N is the size of the sample, V is the covariance matrix of the data sample and h 0 = 1.06 × N −1/6 is the smoothing parameter corresponding to Gaussian samples. The latter is used very often in the literature [39] because it gives very good estimates even for non Gaussian samples. From the density estimate obtained by using Eq. (15) we calculate the parameters, and then we obtain the density estimate from, where h i = h 0 λ i . Fig. 2 shows the average over the 110 density estimates corresponding to the parameters X max and N µ (600) for protons, θ = 45 • and QGSJET-II as the hadronic interaction model.
Every estimatef( x) has an associated uncertainty because it is constructed from a finite set of N elements. To take into account this uncertainty in the composition determination we use the smoothed bootstrap technique [34]. For eachf ( x) we obtain 10 different samples, of the same size used to obtainf( x), by selecting random two-dimensional vectors from it. To every bootstrap sample we perform the same procedure done to the the original sample used to obtainf ( x). Therefore, for every primary type, zenith angle, hadronic model and pair of parameters we obtain 110 estimates of the corresponding distribution function. Note that the smoothed bootstrap technique just allows us to estimate the variance of eachf ( x) but not the bias ( which we assume is of negligible importance because we are using adaptive smoothed parameters. interaction models. From this figure we see that, on average, QGSJET-II produce more muons than Sibyll 2.1 and that the difference between the X max of protons and iron nuclei is smaller for QGSJET-II.

Composition determination
The method presented here is statistical. In principle, we want to infer from a data sample the average composition of the cosmic rays assuming a mixture of proton and iron nuclei, i.e., 0 ≤ c p ≤ 1. We consider samples of N = 100 and N = 1000 events which are the number of hybrid and surface detector events in the energy interval considered, respectively, expected for the 750 m-array in two years of data taking [40]. N = 1000 is also the sample size of the hybrid events that fall in the energy interval under consideration for the life time of Auger, ∼ 20 years.
For each pair of mass sensitive parameters, zenith angle, hadronic interaction model and for each value of c p from 0 to 1 in steps of ∆c p = 0.1, we sample the average proton and iron distributions to generate 1000 independent samples of N = 100 events and 300 independent samples of N = 1000 events.
For each of these samples, corresponding to a given pair of parameters (q 1 , q 2 ), we calculate ξ q 1 and ξ q 2 from, where ν = 1, 2 and Heref k p,ν (q ν ) is the kth marginal distribution function for protons corresponding to the parameter q ν calculated from the kth estimatef k p (q 1 , q 2 ).f l f e,ν (q ν ) is the lth marginal distribution but for iron nuclei. We take QGSJET-II as the "true" high energy hadronic interaction model, therefore, we use just the QGSJET-II estimates to calculate P kl p,ν (q ν ).
For each sample corresponding to a given pair of mass sensitive parameters, hadronic interaction model, zenith angle and proton abundance we obtain 110 × 110 values of ξ qν which correspond to all possible combinations of the marginal distribution functions of proton and iron nuclei. Finally, for each pair of parameters, hadronic interaction model, zenith angle and c p we calculate the mean values ξ q 1 , ξ q 2 , the standard deviations σ(ξ q 1 ), σ(ξ q 2 ) and the covariance cov(ξ q 1 , ξ q 2 ). generated with QGSJET-II and Sibyll 2.1. It also shows the one and two sigma regions of ξ q . As expected, σ is much smaller for N = 1000, because, as explained in section 2, the variance of ξ is inversely proportional to the sample size.
From Fig. 4 we see that the mean value of ξ increases linearly with the proton abundance of the samples, but with slope smaller than one, which corresponds to the ideal case where the iron and proton distributions do not overlap. The slopes of ξ µ and ξ Xmax for the QGSJET-II samples are comparable, meaning that the discrimination power of N µ (600) and X max is similar. This is consistent with the results of Ref. [31], where we showed that the discrimination power of N µ (600) is considerably larger than for X max . But when we take the energy uncertainty into account, they become comparable since showers' muon content depends almost linearly on primary energy.
For samples generated with Sibyll 2.1 (as mentioned, we take QGSJET-II as the reference model) we see that the behavior of ξ q is quite different, in particular for c p = 1/2 the value of ξ q is no longer 1/2. Moreover, the mean value of ξ µ , for a given value of c p , is larger for Sibyll 2.1 since, on average, QGSJET-II produces more muons than Sibyll 2.1 (see Fig. 3). Focusing on X max , we see that for smaller values of c p (larger values of iron abundances), the mean value of ξ Xmax is smaller for Sibyll 2.1. This is due to the fact that X max for iron nuclei is, on average, smaller for Sibyll 2.1 (see Fig. 3). For larger values of c p , ξ Xmax for Sibyll 2.1 is of order or even larger than the corresponding parameter in QGSJET-II. Again, this happens because the mean value of X max for Sibyll 2.1 corresponding to protons is slightly larger than for QGSJET-II. As a final result, the slope of the straight line resulting from Sibyll 2.1 is larger than the slope resulting from QGSJET-II.
We found similar results for θ = 30 • . In particular, in this case, the slope of ξ Xmax is greater than that corresponding to ξ µ , which means that the discriminating power of ξ Xmax is greater than that corresponding to N µ (600). Although the behavior of ξ µ and ξ Xmax for θ = 30 • and θ = 45 • are qualitatively similar, the dependence on the zenith angle is not negligible.
As mentioned in section 2, the distribution functions of the variables ξ q1 and ξ q2 are Gaussian. The mean value and the covariance matrix depend on the proton abundance of the samples and therefore so will the ellipses that enclose regions of a given value of probability. Fig. 5 shows the ellipses corresponding to 68% and 95% probability for the parameters ξ µ and ξ Xmax for our case study: θ = 30 • and 45 • , N = 100 events, for samples built using QGSJET-II and Sibyll 2.1. The evolution of the abundance on the ξ µ −ξ Xmax plane, and the shape and size of the associated ellipses, allow for a smooth estimation of the composition in a way that is reasonably independent of the assumed hadronic interaction model. Furthermore, given two possible interaction models and an observed data set, a figure like the ones depicted in figure 5 can be used to assess simultaneously the compatibility of these models and the experimental data.
Having an experimental sample, we can obtain a point on the ξ µ − ξ Xmax plane by using the density estimates obtained from simulations. This point, together with a diagram like the one in Fig. 5 allow for a quick evaluation of the compatibility between the experimental data and the hadronic interaction models under consideration. Moreover, if the position of the experimental data point on the diagram is compatible with any of the hadronic interaction models, one can also obtain a rough estimation of the composition by simple inspection of the nearest ellipses. As an example, let us consider the top panel of that is located too far from the curves is inconclusive from the point of view of composition or hadronic interaction model, but is a strong indicative of large systematic errors in the detector. more stringent test of the compatibility of experimental data with the hadronic interaction models.
If the experimental point falls close to the region of the ellipses, i.e., the hadronic interaction models are compatible with the data, the composition can be estimated in a more formal way. For this purpose we perform a linear interpolation of the mean values and the elements of the covariance matrices (see Eqs. (4,13)), as a function of c p for each zenith angle, sample size and high energy hadronic interaction model considered. We then assumed a linear dependence in the ξ µ −ξ Xmax plane that links the same value of composition, c p , of the two hadronic models considered ("composition isolines") and choose λ as a variable in the [0, 1] interval which takes the values 0 and 1 for QGSJET-II and for Sibyll 2.1, respectively. This approach is based on the assumption that a continuous and smooth parametric variation between different hadronic models is possible on the ξ µ − ξ Xmax plane.
In this way we obtain, for each zenith angle and sample size, the functions µ(c p , λ) = ( ξ µ (c p , λ), ξ Xmax (c p , λ)) and V(c p , λ). The intermediate values of λ correspond to hadronic models for which the values of the parameters ξ µ and ξ Xmax fall in between those corresponding to QGSJET-II and Sibyll 2.1. Therefore, λ parametrizes the composition isolines for hadronic models that would yield points between those corresponding to QGSJET-II and Sibyll 2.1.
Let us suppose that for a given experimental sample with zenith angle θ and number of events N we obtain ξ exp = (ξ exp µ , ξ exp Xmax ). We can estimate the proton abundance,ĉ p , by solving the equation µ(ĉ p , λ) = ξ exp . The regions in the (c p , λ) plane compatible with ξ exp at a given confidence level are the solutions of the inequality, where R 2 (α) = −2 ln(1 − α) and α is the confidence level (for instance, α = 0.68, α = 0.95, etc.).
In order to obtain the c p intervals compatible, at 95% confidence level, with the considered points, we have to project the regions of Fig. 7 onto the x axis, corresponding to proton abundance. Table 1 shows the inferred composition and its uncertainty at 95% confidence level. Note that until now we just consider the composition of a given sample of cosmic rays, not the one corresponding to the universe (sample of infinite number of events), to see how to obtain the proton abundance of the cosmic rays from the composition of a sample see Appendix B. Table 1 Inferred proton abundance of the samples and its uncertainty at 95% confidence level for the points P i and Q The Auger Observatory measures other composition-sensitive parameters apart from X max and N µ (600) by means of the water Cherenkov detectors. Although these parameters strongly depend on X max and N µ (600), Cherenkov detectors (and muon counters) have 100% duty cycle, as opposed to fluorescence telescopes, which have only 10% duty cycle. This fact highlights the statistical value of surface parameters. For this reason we also studied the applicability of our method to the combination of N µ (600) with other parameters obtained from the Cherenkov detectors: t 1/2 , β and R (see subsection 3.2). Fig. 8 shows the ellipses corresponding to 68% and 95% probability for the combination of ξ µ with ξ t 1/2 , ξ β and ξ R . We consider samples of 1000 events which is the number of events expected in two years of data taking for the 750 m-array of AMIGA. Fig. 8 shows that, also in these cases, the position of the ellipses allows us to test the high energy hadronic models as well as to obtain the proton abundance of a sample, when the experimental data falls in the proximity of the ellipses. Again, the composition uncertainty increases in the region close to c p = 1 especially for θ = 45 • .
A note must be made regarding the potential effects of systematic errors on our technique. Reconstruction biases are rather well known and, very likely, included in a fairly acceptable way in the reconstruction packages used. Furthermore, all the distributions use here are based on reconstructed simulated data. The energy bias, on the other hand, is a fundamental problem for high energy cosmic ray shower measurement at present and it is, certainly, a potential problem for any technique that attempts to infer cosmic ray composition from extensive showers at high energies. Our technique is not an exception in this respect, specially because we use muon number as one of our main parameters which is particularly sensitive to energy. Nevertheless, there are experimental and theoretical facts that can attenuate the problem, or even give it a potentially interesting twist: (i) the bias in energy is bounded to around 30% or, at most, 50% and (ii) the sign of the bias is known, i.e., Monte Carlo simulations overestimate the energy and not the other way around, (iii) the difference in muon content between hadronic interaction models in the ankle region is about 5%, (iv) X max has a logarithmic dependence on energy. This means that, in the ξ Xmax − ξ µ space, curves move in a predictable direction due to the bias and by an amount that is bounded but larger than the expected one due to uncertainties in the hadronic interaction model. Furthermore, the displacement is mostly parallel to the ξ µ axis, without changing appreciably the distance between points of equal composition along curves of constant hadronic interaction model. This implies that, even in the presence of an energy bias, the position of an experimental point on this plane still carries significant physical information. In fact, while one would lose the capacity of discriminating between hadronic interaction models for large systematic biases in energy, one could gain the capacity of experimentally constraining the bias with an uncertainty of the order of 5%.
In any case, the existence of systematic errors is an important and by no means trivial problem, which deserves further studies and the search for more appropriate strategies [41], like the modification or refinement of the parameters actually used in the context of the technique proposed here. Therefore, one alternative line of action, might be, for example, the redefinition of a muon content parameter in the way proposed by Hillas [42].

Conclusions
Detailed composition studies at the energies of the second knee and ankle of the cosmic ray spectrum will be crucial to weight different astrophysical models of the Galactic-extragalactic cosmic ray flux transition. Experiments like Auger and Telescope Array in the near future will be instrumental on the later.
In this paper we present a new statistical method to perform composition studies in a two dimensional space. The method was designed having in mind the enhancements AMIGA and HEAT presently under construction by the Pierre Auger Observatory, but its applicability extends to other detectors, like TALE, which also have hybrid capability.
A main advantage of the method is that it minimizes the effects of the present uncertainty associated with the hadronic interaction models, used to simulate cosmic ray showers, on the inferred composition. Furthermore, in the case in which systematic errors in energy are known or smaller than 5%, the method allows, besides the determination of the composition, an independent verification of the compatibility between real shower data and hadronic interaction models. In the case of larger systematic errors (e.g., ∼ 30%), the technique still allows to make reliable composition estimation and, additionally, to set an upper limit on the size of the systematic error in energy.
The novelty of the Auger enhancements is the combination, at energies between ∼ 10 17 and ∼ 10 18.5 eV, of hybrid measurement with additional simultaneous, and independent muon number information. We exploit this wealth of data by working on a two-dimensional space ξ 1 -ξ 2 which encodes, e.g., X max and N µ information. This space has the property of clearly separating composition and hadronic interaction model dependencies. Confidence levels are calculated in the form ellipses that enclose regions of 68% and 95% probability. It is also shown the way in which the ξ distribution functions decrease as the size of the sample increases in a way compatible with the exposure time scale appropriate for AMIGA-HEAT. We show that, for the case in which the data is not dominated by systematic errors in energy, the constraints imposed on the hadronic models become stronger as the exposure grows while the composition error diminishes to an unprecedented accuracy for astrophysical applications over the lifetime of the experiment. Therefore, besides being able to constrain the high energy hadronic models, it should be possible to determine the composition as a function of energy in the ankle region with errors varying from ∼ 20% to ∼ 5%, at the 95% confidence level, as the data taking progresses from 2 to 20 yr.

Acknowledgments
The authors have greatly benefited from discussions with several colleagues from the Pierre Auger Collaboration, of which they are members. We also want to acknowledge Greg Snow for carefully reading the manuscript and for his valuable comments. GMT acknowledges the support of DGAPA-UNAM through grant IN115707.

A Optimum energy bin
The primary energy is obtained by fitting a lateral distribution function to the total signal in each surface station of the water Cherenkov array. This allows us to interpolate the shower signal at a fixed distance from the core which, in turn, is used as an energy estimator. This reference distance is such that the shower fluctuations go through a minimum in its vicinity, and its exact value depends on the geometry of the array; for Auger, the reference distance is 1000 m for the 1500 m baseline spacing and 600 m for the AMIGA infill of 750 m spacing.
The signal at the reference distance is calibrated with the telescopes via hybrid events. The corresponding energy uncertainty for the 1500 m-array of Auger is ∼ 20% [5]. Guided by this experimental result, we assume in this work a 25% Gaussian energy uncertainty.
In order to study the composition at energies of order E 0 = 1 EeV, we first have to determine the range of reconstructed energies, centered at E r0 , with the form Π r = [(1 − δ)E r0 , (1 + δ)E r0 ] (δ = 0.25 for 25% of energy uncertainty) such that the fraction of events in the interval Π 0 = [(1 − δ)E 0 , (1 + δ)E 0 ] is maximum. The intervals Π r and Π 0 are different because of the spectrum. The contamination of events of real energies smaller than 1 EeV which are outside Π 0 is greater than the one corresponding to energies, also outside Π 0 , but above 1 EeV.
We assume that the number of cosmic ray showers with real energies between E and E + dE is well represented by the following power law spectrum, where, for every set of simulated events we have E 1 = 0.6 EeV, E 2 = 2 EeV, N 0 ∼ = 6640 (depending on the reconstruction efficiency) and γ = 2.7 (see subsection 3.1).
The number of events whose real energy belongs to Π 0 such that the reconstructed energy falls in Π r is given by, On the other hand, the number of events with real energy outside Π 0 whose reconstructed energy falls in Π r is given by, Ratio of the number of events whose real energy does not belong to Π 0 but the reconstructed energy falls in Π r and the number of events whose real energy belongs to Π 0 but the reconstructed energy falls in Π r , F (E r0 ). The minimum is located at E r0 ∼ = 1.14 EeV.
Therefore, the value of E r0 for which the fraction of events belonging to Π 0 that fall in Π r is maximum is obtained by minimizing the function F (E r0 ) = f B (E r0 )/f A (E r0 ). Fig. A.1 shows F (E r0 ) for δ = 0.25 from which we see that it has a minimum but at an energy greater than 1 EeV. The minimum is located at E r0 ∼ = 1.14 EeV, then, Π r = [0.86, 1.43] EeV. where C p is the proton abundance of the cosmic rays.
If we have a sample of N events corresponding to a binomial distribution with n 0 positive trials, an estimator of the parameter p of the binomial formula (C p in the Eq. (B.1)) is given byp = n 0 /N and the upper, p max , and lower, p min , limits of the interval which contains the real value of the parameter p with probability α (confidence level), are the solutions of, for n 0 = 0 and n 0 = N. When n 0 = 0 and n 0 = N the estimators of p are zero and one, respectively, and we can obtain, at a given confidence level α, an upper limit p up for the case n 0 = 0 and a lower limit p low for n 0 = N, (B.5) By using the method described in subsection 3.3 we can find an interval, [c 1 , c 2 ], corresponding to a given confidence level α, for the composition of a sample. Therefore, the number of protons of the sample, at a confidence level α, is contained in [n 1 , n 2 ] = [Nc 1 , Nc 2 ]. For the case where n 1 = 0 and n 2 = N we can obtain the lower limit for the cosmic rays composition, C min p , by solving Eq. (B.3) with n 0 = n 1 and the upper limit, C max p , by solving Eq. (B.2) with n 0 = n 2 . For the case in which n 1 = 0, we just have to calculate the upper limit of the interval by solving Eq. (B.2) with n 0 = n 2 . The same happens for the case where n 2 = N. We just have to calculate the lower limit of the interval by solving Eq. (B.3) with n 0 = n 1 . The criterion adopted to obtain the central value for the cosmic rays composition is to take it equal to the one of the sample,Ĉ p =ĉ p . Table B.1 shows the central values and their uncertainty at 95% confidence level for the composition of the cosmic rays inferred from the composition of the samples obtained for the points P i and Q i with i = 1, 2, 3 shown in table 1. As expected, the uncertainty due to the finite size of the samples is more important for the case of N = 100 events.