Mass Composition of UHECRs from $X_{\rm max}$ Distributions Recorded by the Pierre Auger and Telescope Array Observatories

In this paper we infer the mass composition of the ultra high energy cosmic rays (UHECRs) from measurements of $X_{\rm max}$ distributions recorded at the Pierre Auger (2014) and Telescope Array (TA) (2016) Observatories, by fitting them with all possible combinations of Monte Carlo (MC) templates from a large set of primary species (p, He, C, N, O, Ne, Si and Fe), as predicted by EPOS-LHC, QGSJETII-04 and Sibyll 2.1 hadronic interaction models. We use the individual fractions of nuclei reconstructed from one experiment in each energy interval to build equivalent MC $X_{\rm max}$ distributions, which we compare with the experimental $X_{\rm max}$ distributions of the other experiment, applying different statistical tests of compatibility. The results obtained from both experiments confirm that the mass composition of the UHECRs is dominated ($\gtrsim$$70\%$) by protons and He nuclei {in the energy range investigated $\lg E (\rm eV)$ = [17.8--19.3] (Auger) and $\lg E \rm (eV)$ = [18.2--19.0] (TA).} The indirect comparisons between the $X_{\rm max}$ distributions recorded by the two experiments show that the degree of compatibility of the two datasets is good, even excellent in some high energy intervals, especially above the ankle ($\lg E (\rm eV) \sim 18.7$). However, our study reveals that, at low energies, further effort in data analysis is required in order to harmonize the results of the two experiments.

In this paper we infer the mass composition of the ultra high energy cosmic rays (UHECRs) from measurements of Xmax distributions recorded at the Pierre Auger (2014) and Telescope Array (TA) (2016) Observatories, by fitting them with all possible combinations of Monte Carlo (MC) templates from a large set of primary species (p, He, C, N, O, Ne, Si and Fe), as predicted by EPOS-LHC, QGSJETII-04 and Sibyll 2.1 hadronic interaction models. We use the individual fractions of nuclei reconstructed from one experiment in each energy interval to build equivalent MC Xmax distributions, which we compare with the experimental Xmax distributions of the other experiment, applying different statistical tests of compatibility. The results obtained from both experiments confirm that the mass composition of the UHECRs is dominated ( 70%) by protons and He nuclei in the energy range investigated lg E(eV) = [17.8-19.3] (Auger) and lg E(eV) = [18.2-19.0] (TA). The indirect comparisons between the Xmax distributions recorded by the two experiments show that the degree of compatibility of the two datasets is good, even excellent in some high energy intervals, especially above the ankle (lg E(eV) ∼ 18.7). However, our study reveals that, at low energies, further effort in data analysis is required in order to harmonize the results of the two experiments.

I. INTRODUCTION
The UHECRs (E > 10 18 eV) are the most energetic particles in the universe produced by the most energetic astrophysical objects. Their origin and acceleration mechanisms are still unknown despite the huge efforts in the last decades of the entire astrophysics community [1]. A realistic description of such a complex process of production and transport through the intergalactic medium requires accurate information about the extra-galactic magnetic fields [2,3], interaction with the cosmic microwave background (CMB) and the extragalactic background light (EBL), distribution and type of sources which accelerate these particles [4][5][6][7][8] and, last but not least, accurate information about their energy spectrum [9,10], mass composition [11][12][13][14] and arrival directions [15,16] measured at ground based cosmic ray experiments. These former properties are indirectly reconstructed, making use of different parameters of extensive air showers (EAS) recorded by the cosmic ray experiments such as the Pierre Auger Observatory (Auger) [17] and the Telescope Array (TA) [18].
The mass composition of the primary UHECRs is still a matter of debate even if it was extensively studied by combining different complementary techniques [19][20][21][22][23]. The most reliable parameter from the EAS used to infer the mass composition proved to be X max , the atmospheric depth where the energy deposit profile of the secondary particles reaches its maximum [24]. This parameter is related to the primary particle mass as X max ∝ − ln A.
The measurements of mass composition reported by * nicusorarsene@spacescience.ro the Auger [19,25] and TA [26] experiments, on the basis of the first two moments of the X max distributions ( X max and σ Xmax ) are not in very good agreement on the entire energy range. While the TA results suggest only a light composition above E > 10 18.2 eV, the Auger results clearly indicate a transition to a heavier component starting from E ∼ 10 18.5 eV and becoming increasingly heavier at the highest energies. On the other hand, it was shown that using only the limited information given by the first two moments of a X max distribution may lead, in very particular cases, to a misinterpretation of the mass composition, since different mixes of primary particles can reproduce exactly the same X max and σ Xmax values. To avoid such situations, a method was proposed, which uses the entire shape of each X max distribution by fitting them with Monte Carlo (MC) templates for four fixed primary species (p, He, N and Fe) obtaining in this way information about fractions of individual nuclei [27]. Following this approach, the measurements of X max distributions recorded by Auger and TA experiments were indirectly compared [28] concluding that both measurements are compatible in the limits of statistical and systematical uncertainties.
In a recent paper [29], we showed that by fitting the X max distributions only with four fixed elements (p, He, N and Fe), an artificial worsening of the fit quality can be induced and the reconstructed fractions might be biased as a consequence of a high abundance of some intermediate elements (e.g., Ne/Si) not included in the fitting procedure. We argued that a more appropriate approach is to fit the observed X max distributions with all possible combinations of elements from a larger set of primaries (p, He, C, N, O, Ne, Si and Fe) obtaining the "best combination" of elements that best describe the data.
In this work, we perform an indirect comparison between X max distributions measured by Auger (2014) [19] and TA (2016) [20] experiments following the approach proposed in [29].
In Section II, we describe the simulation procedure for obtaining the MC templates for a large set of primary species (p, He, C, N, O, Ne, Si and Fe), taking into account the detector effects (acceptance and resolution) of both experiments, employing different hadronic interaction models. In Section III, we obtain the fractions of individual nuclei which best describe the X max distributions measured by Auger and TA experiments on the entire energy range lg E(eV) = [17.8-19.3] (Auger) and lg E(eV) = [18.2-19.0] (TA) following two methods. The first method reconstructs the fractions of individual nuclei of the "best combination" of elements which best describe the observed X max distributions. With the second method we extract the average of the reconstructed fractions of each species from all possible combinations of fitting elements whose fit quality was higher than a threshold value, obtaining in this way the evolution of the abundances of all primary species as a function of energy.
In Section IV, we present the compatibility of the measurements of the two experiments on the basis of three statistical tests: p-value as goodness-of-fit, Kolmogorov-Smirnov (KS) and Anderson-Darling (AD). Section IV concludes the paper.

II. MONTE CARLO TEMPLATES
The MC templates (X max distributions or Probability Density Functions (PDFs) of X max ) for 8 primary species (p, He, C, N, O, Ne, Si and Fe) for each energy interval of 0.1 in lg E(eV), employing three hadronic interaction models EPOS-LHC [30], QGSJETII-04 [31] and Sibyll 2.1 [32] were generated using CONEX v4r37 simulation code [33,34]. Due to the experimental limitations, the MC templates for the TA experiment were simulated for only eight energy intervals in the range lg E(eV) = [18.2-19.0] and 15 energy intervals for Auger lg E(eV) = [17.8-19.3]. The zenith angle of the showers were isotropically sampled in the interval θ = [0 • -60 • ] ensuring an isotropic flux on flat surface dN/dcos(θ) ∼ cos(θ). The statistics of each X max distribution are of the order of 10 4 -10 5 events, with a larger number of simulations for proton induced showers in comparison with the heavier nuclei.
A PDF of X max for a primary nuclear species in a given energy interval consists of a binned X max distribution in the range [0-2000] g/cm 2 with a bin width = 20 g/cm 2 for the Auger case, while for the TA experiment the X max distributions are binned in the range [500-1300] g/cm 2 with a bin width = 40 g/cm 2 .
To account for detector effects, the true X max values calculated with CONEX were modified in accordance with the acceptance and resolution of each experiment. In the case of Auger, the true X max values were modified by using Equations (7) and (8) from [19], while for the TA case the true MC templates were modified in accor- dance with the biases (reconstruction + acceptance) and resolutions computed in [20] in Table 1 and Figure 9 of that paper, considering full detector simulations for p, He, N and Fe for the QGSJETII-04 model. For the intermediate elements, the bias and resolution of X max are obtained using a 2nd degree polynomial interpolation. It is worth mentioning that the possible uncertainties on the bias and resolution of the intermediate elements, artificially introduced by this interpolation, would be much smaller than the experimental resolution of the X max parameter. The exact values of the biases and resolutions used to construct the MC templates for the TA MC templates are listed in Table I. For EPOS-LHC and Sibyll 2.1 models we considered the same values of biases and resolutions as computed for the QGSJETII-04 model. Again, we have to stress that this approximation cannot affect the results significantly, since the possible deviations from the true values of bias and resolution for EPOS-LHC and Sibyll 2.1 cannot exceed more than few g/cm 2 . As stated in [20], the bias is larger for deeply penetrating showers, for example, protons. As an example, the difference between the average < X max > induced by protons predicted by Sibyll 2.1 and QGSJETII-04 in the energy interval lg E(eV) = [18.2-18.3] is about 2.4 g/cm 2 , so the bias of the reconstructed X max for protons predicted by these models should be roughly the same. An example of PDFs of X max is given in Figure 1 for proton and iron induced showers in the energy interval lg E(eV) = [18.2-18.3] employing QGSJETII-04 and EPOS-LHC hadronic interaction models and taking into account the specific detector effects of each experiment.
In the next Section we will use these PDFs to fit the experimental X max distributions measured by Auger and TA experiments following a binned maximum-likelihood procedure. The experimental X max distributions in each energy interval recorded by Auger and TA experiments are fitted with all possible combinations of primary elements (p, He, C, N, O, Ne, Si and Fe) following a binned maximumlikelihood procedure. Thus, we obtain the fractions of individual nuclei that best describe the observed X max distributions. In this fitting procedure the minimizing quantity is − ln L, which is defined as: where n i represents the measured counts in the "i"-th bin of the experimental X max distribution and y i is the MC prediction for that bin [35] given by where j = p, He, C, N, O, Ne, Si and Fe, f j (i) is an MC template and c j is the fraction of the j component. We quantify the goodness of fit using the p-value parameter defined as: representing the probability of obtaining a worse fit than the one observed even if the distribution predicted by the fit results is correct. Here, ndf represents the number of degrees of freedom computed as the number of bins (including empty bins) of each X max distribution minus the number of elements considered in the fit, χ 2 is the sum of the square of residuals using the parameters computed by the maximum-likelihood procedure and Γ is the normalized upper incomplete Gamma function.
The results of the fitted fractions of X max distributions measured by Auger and TA are presented following two methods. In method #1 we obtained the "best combination" of elements from all possible combinations that best describe the experimental X max distributions based on the highest p-value max , in each energy interval. The results are displayed in the left plots of Figures 2-4 for EPOS-LHC, QGSJETII-04 and Sibyll 2.1 hadronic interaction models. The error bars of each point were computed using the MINOS technique based on ∆L = 1/2 rule [36]. The reconstructed fractions obtained for both experiments clearly show that the mass composition of primary UHECRs is dominated by light elements (p and He), which present a modulation of their abundances as a function of energy, but keeping the sum (f p + f He ) roughly constant on the entire energy spectrum. This feature is predicted by all hadronic interaction models. An interesting aspect is the presence of Fe nuclei in a quite high abundance (∼20% predicted by EPOS-LHC model) on the entire energy spectrum of TA data which, at least qualitatively, seems to contradict the results of mass composition obtained on the basis of the first two moments of the X max distributions [26].
In method #2 we extract the average of the reconstructed fractions of each species from all possible combinations of fitting elements whose goodness of fit parameter p-value was in the range [0.5 · p-value max , pvalue max ]. For example, if the combination of the elements p + N + Fe and p + He + Si or simply p + He, and so on, give a good fit quality, then the average fraction of protons is the mean of the proton fractions reconstructed from all these combinations. In this case, the error bar associated to one species is computed as the square root of the sum of the squared uncertainties of the species obtained from different combinations of fitting elements whose goodness of fit was in the range mentioned above. In this way, we have an overview of the evolution of all fractions of individual nuclei as a function of primary energy. The results are displayed on the right plots of Figures 2-4 for the three hadronic interaction models. Both methods give similar results on the entire energy range. It is worth mentioning that the shortcoming of method #2 is that the sum of all reconstructed fractions in a given energy interval will be biased (¿1) as a consequence of the limitations of the fitting procedure when using extreme combinations of elements. This bias is evaluated for each energy interval as the sum of all reconstructed fractions, F e i=p f i , for all hadronic interaction models and is represented in Figure 5. As can be seen, the obtained values are not too far from 1 for almost the entire energy spectrum of Auger data, while for TA we observe larger biases in many energy intervals, mainly due to the lower statistics in the experimental X max distributions. Even so, we consider that these bi- ases are not significant compared to the systematic and statistical uncertainties of the experimental X max distributions.
In the next Section, we will perform an indirect comparison between the X max distributions recorded by Auger and TA experiments using the individual fractions of nuclei reconstructed in this Section III.

IV. AUGER VS. TA MASS COMPOSITION COMPATIBILITY
The measurements of X max distributions obtained by both experiments cannot be directly compared because they use different approaches when analyzing the recorded events [37]. This is the reason we produced two sets of MC predictions of PDFs of X max taking into account the specific resolution and acceptance of each experiment.
We test the compatibility of the measurements of Auger data, while in the second approach we translate the individual fractions of each primary species reconstructed from Auger data in equivalent PDFs of X max predicted for TA (PDFs of ) and the comparison is performed between PDFs of X Auger→T A max vs. TA data. In each approach, we use the fitting fractions obtained following method #1 and #2 (see Section III). Each comparison is performed for all hadronic interaction models and is characterized by three statistical tests: p-value as goodness of fit, KS and AD. The p-value parameter is obtained by fitting Auger data with PDFs of X T A→Auger max and TA data with PDFs of X Auger→T A max following the same binned maximum-likelihood procedure. Both KS and AD tests calculate the probability that two distributions come from the same parent distribution. The KS test is one of the most used compatibility procedure with an improved response around the peak of the distribu- tions. On the other hand, the AD test is optimized to probe the compatibility with a higher accuracy on the tails of the distributions which, in our cases, are dominated by protons and He nuclei. In Figure 6 we present the comparison between Auger data and PDFs of X T A→Auger max for the energy interval lg E(eV) = [18.2-18.3] considering the individual fitting fractions obtained from TA data by the two methods (see Section III) for all three hadronic interaction models. The statistics of the X max distribution (N = 1952) as well as the p-value, KS and AD parameters are displayed on the plots. We chose to give this example because it is the energy interval with the highest statistics of events recorded by both experiments. Figure 7 presents the same analysis as in Figure 6 but in the energy bin lg E(eV) = [18. 6-18.7] in which the number of events recorded by Auger is N = 575. In this case, the compatibility of the two datasets is much better, reaching AD = 0.35 for the QGSJetII-04 model. It is important to note that the highest degree of compatibility is obtained with the AD test, which is more sensitive on the tails of the distributions dominated by protons and He nuclei. This is an indication that the measurements of both experiments are more compatible with respect to the light component (p and He) of primary UHECRs.
In Figures 8 and 9 we present the comparison between TA data and the fitted fractions reconstructed from Auger data (second approach) for the energy bins lg E(eV) = [18.2-18.3] and lg E(eV) = [18.6-18.7], respectively.
As in the first approach, the two experiments show that the degree of compatibility of the two datasets is good, reaching excellent agreement in some high energy intervals, especially when using the fractions reconstructed with method #2 (e.g., p-value= 0.83, KS = 0.94 and AD = 0.86) for lg E(eV) = [18.6-18.7] considering QGSJetII-04 model. At lower energies, our study shows that the reconstruction methods need to be refined in order to achieve a better compatibility between the two sets of measurements.
The complete set of probabilities computed with the three statistical tests is listed in Tables III and II for the first approach and second approach, respectively.

V. DISCUSSIONS AND CONCLUSIONS
In this paper, we present an alternative approach to inferring the mass composition of the primary UHECRs using the available X max distributions recorded by Auger (2014) and TA (2016) experiments, by comparisons with MC templates predicted by EPOS-LHC, QGSJETII-04 and Sibyll 2.1 hadronic interaction models.
A common general remark is that the measurements of both experiments suggest that the mass composition of primary UHECRs is dominated ( 70%) by protons and He nuclei, which present a modulation of their abundances as a function of energy but keeping the sum (f p + f He ) roughly constant on the entire energy spectrum. This conclusion holds for all the three hadronic interaction models. An interesting aspect is the presence of Fe nuclei in a quite high abundance (∼20%) in TA data on the entire energy spectrum lg E(eV) = [18.2-19.0] as predicted by the EPOS-LHC hadronic interaction model.
The general conclusion is that, using the approaches proposed in this paper, the degree of compatibility of the two datasets is good, reaching excellent agreement in some high energy intervals above the ankle, as com-puted by the three statistical tests (e.g., p-value= 0.83, KS = 0.94 and AD = 0.86 for lg E(eV) = [18.6-18.7] considering the QGSJetII-04 model). The best agreement is achieved when the indirect comparison is performed using the reconstructed fractions of nuclei from Auger data to build the equivalent distributions, PDFs of X Auger→T A max , which are compared with TA data. However, our study reveals that, at low energies, further effort in terms of data analysis is required in order to harmonize the results of the two experiments, taking into account that the models are not able to accurately describe the data in some energy intervals.
We believe that the current approach could be very useful in future studies on mass composition, especially to crosscheck the measurements of X max between Auger and TA once the statistics of the recorded events will increase. 17