Measurements of UltraHigh Energy Cosmic Rays UHECR suggest a complex composition with significant contributions from heavy nuclei at the highest energies. We systematically explore how the selection and number of primary nuclei included in the analysis impact the inferred UHECR mass composition. Introducing a distance measure in the space of Xrm max distribution moments we demonstrate that limiting the analysis to a few primaries can introduce significant biases particularly as observational data improves. We provide lists of primaries approximately equidistant in the new measure which guaranty unbiased results at given statistical confidence. Additionally we explore consistent inclusion of nuclei heavier than iron and up to plutonium deriving first observational upper bounds on their contributions to UHECR with the Pierre Auger Open Data.
['Astrophysics - High Energy Astrophysical Phenomena', 'High Energy Physics - Phenomenology']
Resolution of Heavy Primaries in Ultra High Energy Cosmic Rays
{'Resolution of (Heavy) Primaries in Ultra High Energy Cosmic Rays': 'Blaž Bortolato, 1, ∗ Jernej F. Kamenik, 1, † and Michele Tammaro 2, ‡  \n1 Jožef Stefan Institute, Jamova 39, 1000 Ljubljana, Slovenia  \nFaculty of Mathematics and Physics, University of Ljubljana, Jadranska 19, 1000 Ljubljana, Slovenia 2 INFN Sezione di Firenze, Via G. Sansone 1, I-50019 Sesto Fiorentino, Italy  \n(Dated: September 12, 2024)  \nMeasurements of Ultra-High Energy Cosmic Rays (UHECR) suggest a complex composition with significant contributions from heavy nuclei at the highest energies. We systematically explore how the selection and number of primary nuclei included in the analysis impact the inferred UHECR mass composition. Introducing a distance measure in the space of X max distribution moments, we demonstrate that limiting the analysis to a few primaries can introduce significant biases, particularly as observational data improves. We provide lists of primaries approximately equidistant in the new measure, which guaranty unbiased results at given statistical confidence. Additionally, we explore consistent inclusion of nuclei heavier than iron and up to plutonium, deriving first observational upper bounds on their contributions to UHECR with the Pierre Auger Open Data.', 'I. INTRODUCTION': "Measurements of Extensive Air Showers (EAS) at Pierre Auger Observatory (PAO) [1] and at Telescope Array [2] indicate that the mass spectrum of Ultra-High Energy Cosmic Rays (UHECR) is best explained by a mixture of different primaries. Their inferred composition also depends on the cosmic ray energy: at lower energies ( E ∼ 10 18 eV), it seems to be dominated by protons and possibly other light elements, while measurements at higher energies ( E ∼ 10 19 -10 20 eV) point towards a significant fraction of heavy nuclei [3-6].  \nThe UHECR composition is closely related to other open questions in the field, regarding the sources of UHECRs, their acceleration mechanisms and their propagation in the galactic and intergalactic medium.  \nTypically, the nuclei assumed to be part of UHECRs (apart from protons) span from stable isotopes of helium (He, A = 4 ) to iron (Fe, A = 56 ), as Fe has the largest binding energy per nucleon and it is likely the end product of stellar evolution. While heavier elements ( A > 56 ) are rarer in cosmic radiation, they can be produced at UHE, e.g. in GRB [7], and have a similar propagation phenomenology as lighter elements [8-11]. Thus it is reasonable to explore scenarios where these super heavy nuclei make up a fraction of the UHECR and to study their observational features at observatories like the PAO. Recent studies indicate the presence of Z ∼ 53 primaries as a possible origin of the 'second knee' in the UHECR energy spectrum [12, 13].  \nAt terrestrial observatories, the composition of UHECRs can be inferred by studying the longitudinal profile of an EAS, that is the intensity of fluorescent light emitted by nuclei in the atmosphere, typically nitrogen, excited by the passage of charged particles, and measured  \nas a function of the s.c. slant depth ( X ) of the shower. In general, the longitudinal profile has a clear peak, at X max , corresponding to the point of maximum population of electrons and positrons in the shower evolution. Comparing observed data on X max distributions at a particular energy against simulations for a set of primaries, one can infer the allowed fraction of each considered element in the observed collection of UHECRs.  \nPast studies of UHECR composition have considered either a binned distribution of X max or its global features such as its mean ⟨ X max ⟩ and width σ ( X max ) [3, 5]. Typically the measured values of these were compared to simulations of a few chosen primaries. Recently [14], we have have shown how a systematic description of the X max distribution in terms of its central moments together with efficient inference techniques based on bootstrapping and nested sampling can be used to infer the full composition of all stable primaries up to iron. We have also shown that while the fractions of individual primaries in this case can be highly correlated and thus degenerate with each other, one can still derive robust confidence intervals for fractions of all primaries heavier than any chosen nuclear mass or atomic number.  \nIn this work we extend our previous studies of UHECR composition by addressing two main questions concerning the choice of considered primaries when inferring composition, namely how does the (1) range of primaries and (2) their total number affect the inferred composition confidence intervals. This is particularly important for both present and future studies. Firstly, the current results are based on the (strong) assumption that a mixture of only four primaries, typically (p,He,C/N/O,Fe) [15], describes well the UHECR mass spectrum. Secondly, it is warranted in order to to derive bounds on possible fractions of primaries beyond iron, because a brute-force approach of including all stable enough primaries up to uranium ( Z = 92 ) or plutonium ( Z = 94 ) is computationally extremely expensive using existing techniques. To this end we first introduce a natural metric in the space of X max distribution moments. Because simulated X max distributions of primaries have intrinsic uncertain-  \nties, this allows on one hand to identify a-priori practically (given finite event statistics) indistinguishable primaries, and on the other allows us to include only distinguishable primaries in the inference procedure. Conversely, we study in detail using mock data, how omission of distinguishable primaries can potentially lead to biased inference of primary UHECR composition. In particular, we show that mass composition results of existing studies using the full PAO statistics but with limited number of primaries below Z = 26 are potentially biased. Finally we explore the effects of including primaries beyond Z = 26 in the composition and in particular derive first upper bounds on fractions of UHECR primaries up to plutonium using the PAO open dataset.  \nThe remainder of the paper is structured as follows: in Sec. II we introduce our methodology and data, in particular the PAO open dataset, the CORSIKA EAS simulation framework, the relevant observables - moments of X max distributions, and our statistical inference procedure. In Sec. III we study the implications of inferring composition using different numbers of primaries. We first define a distance measure in the space of X max distribution moments and show how it can be used to define distinguishable primaries. We apply this measure to select primaries to be included in the inferred mixture. We test this procedure both using mock data as well as the PAO open dataset, where we also infer composition including primaries beyond iron and up to plutonium. Our main conclusions are summarised in Sec. IV, while the Appendices A and B contain some auxiliary computations, as well as additional tables and plots, respectively.", 'II. METHODOLOGY': "a. Open Data : we use publicly available data from the 2021 Pierre Auger Open Data (PAOD) release [16]. This dataset consists of 10% of the full recorded data by the PAO. In particular, we focus on so called hybrid showers: events recorded by both the Surface Detectors (SD) and the Fluorescent Detectors (FD). There are 3348 hybrid showers, distributed in the reconstructed primary energy as ∼ E -2 . 6 , from E min = 0 . 65 EeV to E max = 63 EeV. We select the energy bin E ∈ [0 . 65 , 1] EeV, containing 934 events, and use it in all our remaining discussion. In this way we work with a fairly large statistical sample within a small energy bin. For completeness, we include the respective results in the energy bin E ∈ [1 , 2] EeV, encompassing 1234 events, in Appendix B.  \nb. Simulations : we simulate EAS with CORSIKA v7.7550 [17]. To achieve fast and reliable simulations of longitudinal shower profiles, we employ the CONEX option. In conjunction with the EPOS hadron interaction model, this option allows to include primary nuclei up to A = 295 . Note however that composition results can strongly depend on the choice of the hadronic model [14]. Thus, our results cannot explore the potential effects of different parametrizations of UHE hadronic cross sec-  \ntions. We simulate 10 4 showers per primary, with the energy uniformly distributed in the chosen energy bin. The atomic number A of the primary is taken as the one of the most abundant stable isotope; we thus unambiguously indicate each primary nucleus by its atomic number Z . We simulate all primaries from Z = 1 (proton) to Z = 26 (iron); in addition, we simulate showers for heavier primaries with Z = 27, 28, 29, 30, 34, 39, 40, 42, 43, 44, 47, 50, 57, 58, 64, 72, 81, 82, 91, 92, 94. We chose to simulate a set of primaries with atomic number immediately beyond iron, and additionally a set of primaries that cover the interval up to Plutonium, Z = 94 , which is the heaviest that can be simulated with CONEX. Finally, we extract the value of X max from each simulated shower by fitting the longitudinal profile to a GaisserHillas curve [18]. The set of X max per primary furnishes a distribution p ( X max | Z ) .  \nc. Observables : we describe the resulting X max distribution by its moments, defined as  \n⟨ X n max ⟩ = ∫ p ( X max | Z ) X n max d X max ∫ p ( X max | Z )d X max . (1)  \nFrom this expression we can obtain the central moments as  \nz n = ⟨ ( X max -⟨ X max ⟩ ) n ⟩ = n ∑ k =0 ( n k ) ⟨ X n -k max ⟩ ( -1) k ⟨ X max ⟩ k . (2)  \nThe number n of moments needed to infer the composition depends on the interplay of statistical and systematic uncertainties, which we discuss in the next paragraph. We have shown in Ref. [14] that with the data available in the Auger Open Data, n = 3 is sufficient and additional moments do not improve the results.  \nd. Uncertainties : We consider both systematic and statistical uncertainties on X max . The former originate from the finite resolution of the fluorescent detector [19]; we include these by reweighting the distribution p ( X max | Z ) by the detector smearing [14]. The latter are a consequence of the limited statistics available, either in the simulation or in the observed data. We take these into account by means of bootstrapping [20], that is repeatedly sampling sets of X max from p ( X max | Z ) and computing moments of these sets; the spread of the moment distributions obtained will depend on the total number of events available.  \nThe next observational run of Auger Prime [21] is expected to deliver a larger and more precise dataset. In the following we perform projections to larger statistics by assuming that the central values of the measured X max moments remain unchanged, while the variances of moment distributions scale with the statistical multiplier f as 1 /f , when increasing the total available statistics by a factor of f . 1 While such scaling is formally only valid  \nfor statistical uncertainties, at the current and even projected PAO statistics, also the dominant systematic uncertainties are expected to effectively scale with statistics. Namely, the detector resolution is calibrated using (calibration) data and the resulting uncertainties are statistics dominated [19].  \ne. Inference procedure : we follow the steps of our work, Ref. [14], to infer the composition of UHECR, w . For a given mixture of primaries, the log-likelihood is  \nln L ( w ) = ∫ ln [ p ( z | µ w , Σ w ) ] p ( z | ˜ µ, ˜ Σ ) d z , (3)  \nwhere p ( z | ˜ µ, ˜ Σ ) is the moment distribution of the Auger data, while p ( z | µ w , Σ w ) is the distribution of moments for simulated momenta, weighted by the composition w . The number of free parameters depends on the set of Z allowed to be in the mixture, which is worked out in Sections III A and III. The best fit composition can be obtained by maximizing the log-likelihood in Eq. (3). However, the typical number of free parameters and thus the dimensionality of the likelihood manifold is large. In order to be able to effectively extract the confidence intervals of our best fits, that is the shape of ln L around the maximum, we use Nested Sampling (NS) techniques [23, 24]. The latter allow us to sample highdimensional likelihoods, with up to 100 free parameters, thanks to a large reduction in computational time with respect to traditional sampling methods.  \nWe present the results as a cumulative composition : for each Z included in the fit, we show the allowed fraction of primaries heavier than Z . This method yields a reliable presentation of the results, including well defined confidence intervals that account for correlations among individual primaries' fractions. It also allows to compare composition results with different numbers of primaries since, contrary to individual primary fractions, the cumulative fraction values for any Z are always well defined and have the same statistical meaning irrespective of the actual primaries included in the fit (as long as protons are included).", 'III. COMPOSITION WITH DIFFERENT NUMBERS OF PRIMARIES': 'The procedure of inferring the composition of UHECR, introduced in Ref. [14] and summarized in Sec. II raises an important question: what is the minimal number of composition components required to obtain an unbiased fit to data? Namely, defining the composition vector w = ( w Z 1 , w Z 2 , . . . , w Z k ) , we wish to quantify the minimal length k and find a suitable list of Z k primaries, such  \nthat the inference procedure yields (cumulative) composition fractions consistent with their truth level values. The set of primaries (p,He,C/N/O,Fe), that is N L = 4 , has been largely assumed in the literature [4, 15] to represent the UHECR mass spectrum. Previous works, see Refs. [5, 25], have however shown that the goodness of fit for the UHECR composition is strongly affected by the number of primaries chosen, and that it generally improves when including more intermediate elements, up to N L = 8 . Here we extend these studies in a systematic way, employing a distance measure in the space of X max moments to provide a quantitative criterion for the selection of primaries.', 'A. Primary Distance in Space of X max Moments': "The decomposition of the X max distributions in terms of a series of central moments allows us to quantify the separation between different primaries Z . Given two primaries Z 1 and Z 2 , we define the distance  \nd 2 n ( Z 1 , Z 2 ) ≡ ( µ Z 1 -µ Z 2 ) T (Σ Z 1 +Σ Z 2 ) -1 ( µ Z 1 -µ Z 2 ) , (4)  \nwhere n indicates the number of moments used. Here µ Z 1 , µ Z 2 are the n -vectors of moment means and Σ Z 1 , Σ Z 2 the respective n × n covariance matrices 2 . This definition builds upon the property of the bootstrapped distributions of moments, that they quickly approach normal distributions. The same feature allows for a straightforward statistical interpretation of d 2 n : the distance behaves as a chi-squared distribution with n degrees of freedom, d 2 n ∼ χ 2 n , thus its value can be converted to a confidence level of distinguishing among two sets of moments. For example, if we fix n = 1 moments and the distance between two primaries is d 2 1 ∼ χ 2 1 ≃ 2 . 7 , then we can interpret this as 90% confidence that the compositions consisting of the two primaries are in principle distinguishable using measurements of only the first moment.  \nSince X max distributions of individual primaries cannot be inferred directly from data, distances d 2 n need to be estimated using simulations. The main systematics in calculating the moments and in particular their means is then expected to come from the parametrization of UHE hadronic interactions, i.e. the choice of the hadronic model. In addition however, the corresponding covariance matrices also scale with the simulation statistics. This means that the separability of pairs of primaries crucially depends on the assumed (simulated) data statistics. In particular, for N simulated or measured events, the distances d n scale as √ N . The scaling remains true even in the presence of systematic uncertainties on the  \nX max measurements. As detailed in Ref. [19] and discussed in Sec. II, the main systematics namely come from the limited statistics of the observed data.  \nIn the top plot of Fig. 1 we show the n = 1 distances d 1 for the full list of our simulated primaries in the chosen energy bin, assuming the same statistics as PAOD. As expected, the distance for a fixed Z 1 grows with the difference | Z 2 -Z 1 | as the distributions (and consequently their moments) of X max drift apart for more distant primaries. Somewhat surprisingly however, we find that the distance d n between any two primaries is to a very good approximation independent of the number of moments considered n . This observation can be traced back to the covariance matrices Σ Z , which exhibit a high degree of correlation between the different moments, as first observed already in Ref. [14]. This property allows us to disregard higher moments in the context of distances, and consider only the first moment z 1 for all the simulated primaries. In addition, the monotonic nature of z 1 allows us to employ an analytic interpolation formula for the mean and variance of z 1 as functions of Z (see Appendix A for details) in order to compute d 1 for all pairs of primaries (including those not explicitly simulated) up to Z = 94 , as shown in the bottom plot of Fig. 1.  \nLeveraging on the statistical interpretation of the distance between primaries and making use of these interpolating functions, we can obtain lists of approximately equidistant primaries. We fix the first element to always be proton, Z = 1 , and choose a step distance d 0 . We then iteratively find the next element of the list as the one whose distance from all previous is at least d 1 ≥ d 0 . Each choice of d 0 thus directly corresponds to a set of Z and defines the level of separation between neighbouring primaries in the shower composition.  \nWe show some examples in Table I, assuming 10 3 events per primary, where we also indicate the number of primaries with Z ≤ 26 and Z > 26 as N L and N H respectively. Note that a set with N L = 4 , as widely employed in the literature, would correspond to a relatively large distance step, d 0 = 16 . 6 at these statistics, while asking for at most d 0 < 2 separation leads to N L = 20 . Due to the scaling of d n with assumed primary statistics, these lists can be used at any given statistics if we also rescale the threshold value d 0 by the corresponding statistical multiplier √ f . For example, by increasing the dataset tenfold, f = 10 , to 10 4 events per primary, the same lists are obtained by imposing the new thresholds d ' 0 = √ fd 0 ∼ 3 d 0 . For completeness, we show the behaviour of N L and N H as functions of d 0 assuming 10 3 events per primary in Fig. 2.  \nWe note that the distributions of X max moments change with CR energy and thus both N L,H as well as the lists of d 0 -equidistant primaries can change for the same d 0 in different energy bins, see e.g. Table II, so these should be recomputed for any specific CR energy range considered.  \nFinally, in any measured CR event sample of unknown composition, only the total number of registered showers  \nFIG. 1: Top: distance d 1 for each pair of simulated primaries. Bottom: distance d 1 for each pair of primaries, using the interpolating functions in Eq. (A1). Note that the scale is linear for d 1 ∈ [0 , 1] while it is logarithmic above this interval. See text for details.  \n<!-- image -->  \n<!-- image -->  \n1  \nis known a-priori, and not the statistics of individual primary components. Thus, in practice one can only compute an upper bound on d 2 n ( Z 1 , Z 2 ) by assuming that the dataset set is composed of only two primaries ( Z 1 and Z 2 ) of equal fractions. The actual separability within the dataset will then always be lower and depend on the (unknown) true composition. Our d 0 values and lists in Tables I and II can thus be considered as very conservative if applied to event samples of given total statistics. They nonetheless provide meaningful benchmarks, as we demonstrate using simulated mock data examples of known composition in the next section.", 'B. Composition Inference on Mock Data': 'We demonstrate the implications of different d 0 based w length choices for the mixture fit results using specific mock compositions of UHECR based on simulated data. We first select a set of primaries and their respective frac-  \nd 0 N L ( N H ) List of atomic numbers ZTABLE I: Lists of approximately equidistant primaries in the chosen energy bin. The first column shows the chosen step d 0 assuming 10 3 events per primary, while the second indicates the values of N L and N H obtained. The last column lists all selected Z .  \nFIG. 2: Stacked N L (gray) and N H (black) values as a functions of the distance threshold d 0 , assuming 10 3 events per primary in the energy bin [0 . 65 , 1] EeV.  \n<!-- image -->  \ntions to represent the truth-level composition, w true . We compute the moments of the truth-level X max distribution by randomly sampling 1000 events, apportioned following w true . In order to ensure consistent truth level  \nOn the other hand, we observe that this behavior is not universal and that the relative goodness of fit does  \ncompositions at large values of f , we do not separate these events from the rest of the simulations. Otherwise, the computed values of the distribution moments would be slightly different from the ones computed using all events; thus at large f , this difference could lead to significant deviations from w true .  \nWe study two example compositions:  \nEx1 w 1 = w 10 = 1 / 2 , that is we assume the mass spectrum consists of equal fractions of only two primaries, Z = 1 and Z = 10 ;  \nEx2 w 1 = w 2 = · · · = w 10 = 1 / 10 , that is we take equal fractions of the lightest ten primaries.  \nFor each example, we infer the composition using the lists of light primaries 3 with N L = 4 , 8 , 12 , 16 , 20 , 24 , 26 , as shown in Table I, and for several values of the statistical multiplier f . The results are presented in Fig. 3 for Ex1 (top) and Ex2 (bottom), respectively. On the left side we show the truth level cumulative composition w true (in black lines), 4 and compare it to the 95% confidence intervals for three choices of N L (note that we only show the cumulative fractions at the Z values included in w for each choice); on the right side we plot the differences of maximum ln L values (corresponding to best fit compositions) with respect to the value obtained with N L = 26 , ln L 26 , that is by considering all light primaries in the composition.  \nIn both examples and for all choices of N L , cumulative composition confidence intervals shrink with increasing statistics ( f ), as expected. However, for too small N L (depending on composition and statistics) the 95% confidence intervals fail to cover the truth-level values w true thus indicating biased inferred compositions.  \nAn extreme example is given by the N L = 4 case in Ex2. In the left plot, the red bars show an increasing precision on the cumulative fractions with increasing statistics. However, the intervals are systematically pulled away from the truth-level values indicated by the black line. A significant deviation is present already at f = 1 , that is assuming 10 3 measured events in total, and it quickly grows with larger statistics. The N L = 4 case thus leads to an underestimation of the heavier primary fractions and to biased results in general, even for event statistics comparable to PAOD. At the same time the relative goodness of fit degrades significantly with increasing f , as seen in the right plot.  \nnot necessarily degrade significantly with increasing f . It can also remain comparable to the full N L = 26 case even as the composition results become biased. This can be most clearly seen for the N L = 8 case in Ex1, where the log likelihood difference remains almost constant even as the fraction of the heavier element in the composition is significantly underestimated. Interestingly, in Ex1, the same happens even for N L = 16 , where at f = 5 the fit  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFIG. 3: Cumulative composition confidence intervals and relative goodness of fit for Ex1 ( top ) and Ex2 ( bottom ) mock UHECR compositions. Left : the black line indicates the truth-level composition, while colored bars indicate 95% confidence level intervals for the respective primaries. Red, green, blue and gray bars refer to the inference with N L = 4 , 8 , 16 , 26 primaries respectively, while shades of the same color indicate the multiplier f , as specified by the legend on the right. Right : Difference of ln L values of best fit compositions, with N L = 4 , 8 , 12 , 16 , 20 , 24 , with respect to the case N L = 26 . See text for details.  \n<!-- image -->  \nstarts to underestimate the heavier primary fraction and thus actually yields biased composition results. Note that this choice of N L corresponds to d 0 ≃ 6 assuming 5 × 10 3 events per primary. Although the true composition in Ex1 might be considered extreme, it nonetheless demonstrates that existing mass composition studies based on the full PAO dataset which all consider N L < 16 , potentially yield 95% confidence intervals of primary fractions which do not cover their true values.', 'C. Application to Pierre Auger Open Data': 'Building upon the mock data examples studied in the previous Section, we next infer the composition of UHECR from PAOD in the chosen energy bin. We show our results for the cumulative fractions in Fig. 4, for different choices of N L and at various projected statistics factors f (actual PAOD statistics correspond to f = 1 ). Note again that for the moment we are restricting the procedure to lists of primaries up to iron, Z ≤ 26 (and so N H = 0 ). Contrary to mock simulated data, it is of course impossible to know a priori the true mass composition of the UHECR spectrum. However, we can use the results of the full N L = 26 fit, in particular its maximum  \nFIG. 4: Cumulative fraction of primaries, with 95% CL, obtained with N L = 4 , 8 , 16 , 26 (red, green, blue, gray), using PAOD. The solid black line indicates the best fit obtained with N L = 26 .  \n<!-- image -->  \nlikelihood composition w best , as reference against which we compare other N L choices. We observe that a familiar pattern emerges: including only N L = 4 primaries, the fit tends to underestimate the cumulative fractions obtained with the full 26 primaries, especially the total fraction of nuclei with Z > 1 . Even more pronounced however is the underestimation of the 95% confidence intervals, very similar to what is observed in mock data, where it signals a biased composition. The N L = 16 fit results instead are already well compatible with the reference N L = 26 case including the total coverage of the 95% CL intervals, at least at PAOD statistics (for f = 1 ).', 'D. Bounds on Admixtures of Nuclei Beyond Iron': "Finally, we study the effect of including possible 'heavy' primaries, that is nuclei with Z > 26 . We have seen that, in order to have an unbiased estimation of the cumulative fraction of primaries in PAOD, we require at least N L = 16 . We can then use the same procedure to extend our lists to include approximately d 0 -equidistant primaries up to Z ≤ 94 , as shown in Table I.  \nTo portray the effect of adding N H heavy primaries on", 'IV. CONCLUSIONS': "In studies of the UHECR mass spectrum, the vector of primary composition fractions, w , fundamentally drives the inference procedure: its size defines the dimensionality of the problem, and as a consequence the computa-  \nthe likelihood, we first infer the composition after adding two consecutive d 0 -equidistant heavy primaries at a time, following the list for N L = 16 in Table I; the results are shown in the left-hand plot of Fig. 5. The bounds on cumulative fractions do not change significantly, as it can be seen by comparing the black bands ( N H = 0 ) to the colored ones; the inferred CL intervals increase only slightly, due to the addition of new free parameters with the same statistical sample size. As a direct result, we find that the fraction of primaries with Z > 26 in the energy interval E ∈ [0 . 65 , 1] EeV is within 95% CL bounded by  \nw ( Z > 26 , E ∈ [0 . 65 , 1] EeV) ≤ 24% , (5)  \nwhile for E ∈ [1 , 2] EeV we get  \nw ( Z > 26 , E ∈ [1 , 2] EeV) ≤ 18% . (6)  \nBy construction, the upper bounds on heavy primaries are decreasing monotonically with increasing Z . Consequently, bounds on so-called super heavy nuclei, such as uranium and plutonium can be much tighter than the above values. As an illustration, we thus also make a conservative estimation on fractions of such primaries as  \nw ( Z > 81 , E ∈ [0 . 65 , 1] EeV) ≤ 10% , w ( Z > 85 , E ∈ [1 , 2] EeV) ≤ 6% , (7)  \nboth at 95% CL. These values are obtained using the full d 0 -equidistant range of N L = 16 and N H = 10 primaries on PAOD (see Tables I and II for the lists).  \nA different outcome would be obtained if instead we only added super heavy primaries, e.g. Z = 81 , 91 , to the original 16 primaries with Z ≤ 26 . This is shown in the right-hand plot of Fig. 5. In particular, we observe a substantial underestimation of the iron fraction compared to the results with d 0 -equidistant primaries: both for N H = 10 , as well as with N H = 2 , where the consecutive heavy primaries Z = 30 and Z = 34 are included in the fit. We thus conclude that at least when applied to PAOD, the method of including consecutive d 0 -equidistant primaries with Z > 26 yields reliable upper bounds on cumulative fractions of heavy elements in UHECRs even when the full range of heavy nuclei is not included in the fit. Conversely, including only ultraheavy primaries in the fit in addition to Z ≤ 26 can lead to biased composition results.  \ntional power required to address it; the specific primaries that form w potentially impose a strong prior on the inferred composition, which in turn can lead to biases in the fit results. It is thus important for the experimental collaborations determining the mass spectrum to base their procedures upon a solid and quantitative framework for  \n<!-- image -->  \nFIG. 5: Cumulative fraction of primaries in PAOD at 95% CL. Right : Results with N L = 16 (black), and adding N H = 2 , 4 , 6 , 8 , 10 (red,yellow,dark blue, light blue, cyan) sequential d 0 -equidistant heavy primaries. Left : Results with N L = 16 and N H = 2 ultra-heavy primaries ( Z = 81 and Z = 91 ) (black), and N H = 10 d 0 -equidistant heavy primaries (gray). See text for details.  \n<!-- image -->  \nselecting optimal w .  \nIn order to provide such framework, we can leverage the discriminating power of the X max distributions, to define a measure of distance between primaries in the space of X max moments. By construction, the distance gives the statistical significance to distinguish between two primaries in a composition, as function of the statistical sample size. In other words, for normally distributed moments, this distance is equivalent to the significance, in units of σ , that the two sets of moments are distinct from each other. Furthermore, at least for the set of simulated primaries, this distance is well approximated by including only one (i.e. the first) moment. Consequently, we can build lists of approximately equidistant primaries for the vector w ; for convenience, we separate these primaries in 'light', with Z ≤ 26 , and 'heavy', with Z > 26 , and indicate their numbers in the lists as N L and N H . Taking 10 3 events per primary as an example relevant for the statistics available in the PAOD, a minimum distance of d 0 ∼ 3 between primaries implies w to be at least of length N L = 16 . Given the typical scaling with the statistical sample size, requiring the same d 0 with the full PA statistics would lead to even larger dimensionality for w .  \nWe first tested this framework on two sets of simulated mock compositions. By fixing the underlying true w , we can observe the effect of imposing different N L on the inferred composition. The examples demonstrate that too low values of N L for given sample size (or equivalently too high d 0 ) in general lead to biased composition fraction estimates with underestimated uncertainties, see Fig. 3. As noted previously [5, 25], the relative goodness of fit, computed as the maximum value of the log-likelihood, in some cases tracks this increasing discrepancy of low N L  \nresults with respect to sufficiently large N L values. However we find that this behavior is not universal and that in general a high goodness of fit score does not guarantee unbiased composition results.  \nA similar pattern is observed when applying the same strategy to PAOD, shown in Fig. 4; here the N L = 4 best fits seem to underestimate the fraction of primaries other than protons, with respect to the N L = 16 or N L = 26 results. The effect is expected to be statistically significant for the existing full PA dataset. We thus recommend a minimum of N L = 16 w components for the inference of the UHECR mass composition, assuming that no elements heavier than iron are present.  \nBy lifting the assumption that Z ≤ 26 , we can repeat the procedure including so-called 'heavy' primaries. We explored possible ways to add N H heavy primaries to w . The results shown in Fig. 5 indicate that adding very distant primaries; in the example presented, the introduction of only Z = 91 and Z = 94 in addition to Z ≤ 26 strongly biases the fit. Fortunately, the bias can be reliably avoided by instead including sequential equidistant heavier primaries starting from a given N L set.  \nAt cosmic ray energies considered in this work, no heavy elements are actually expected, as the data favors a mixture of protons and lighter nuclei. Using PAOD we nonetheless put first observational bounds on the presence of primaries heavier than iron, as well as more specifically on elements heavier than uranium and plutonium in UHECRs in the energy range E ∈ [0 . 65 , 2] EeV. The same analysis cannot be carried out at higher energies, as the PAOD does not contain a meaningful number of events above ∼ 2 -5 EeV.  \nOur work could be extended in several directions. Besides applying our methodology to the full PA (and Tele-  \nscope Array) datasets, the determination of an optimal composition vector w for a given CR dataset of unknown composition remains an open problem. We have argued that, while necessarily conservative, interpreting d 0 values computed with individual primary statistics comparable to the full dataset as resolution significances (measured in normal distribution sigmas) gives meaningful benchmark values to determine both the size and the relevant components of w . However, the potential reduction of computational cost as well as inferred composition  \n- [1] Pierre Auger, A. Aab et al., Nucl. Instrum. Meth. A 798 , 172 (2015), 1502.01323.\n- [2] H. Kawai et al., Nuclear Physics B - Proceedings Supplements 175-176 , 221 (2008), Proceedings of the XIV International Symposium on Very High Energy Cosmic Ray Interactions.\n- [3] Pierre Auger, A. Aab et al., Phys. Rev. D 90 , 122006 (2014), 1409.5083.\n- [4] P. Lipari, Phys. Rev. D 103 , 103009 (2021), 2012.06861.\n- [5] N. Arsene and O. Sima, Eur. Phys. J. C 80 , 48 (2020), 2001.02667.\n- [6] Pierre Auger, A. Aab et al., Phys. Rev. D 96 , 122003 (2017), 1710.07249.\n- [7] B. D. Metzger, D. Giannios, and S. Horiuchi, Mon. Not. Roy. Astron. Soc. 415 , 2495 (2011), 1101.4019.\n- [8] L. N. Epele and E. Roulet, JHEP 10 , 009 (1998), astroph/9808104.\n- [9] G. Bertone, C. Isola, M. Lemoine, and G. Sigl, Phys. Rev. D 66 , 103003 (2002), astro-ph/0209192.\n- [10] D. Allard, Astropart. Phys. 39-40 , 33 (2012), 1111.3290.\n- [11] B. T. Zhang, K. Murase, N. Ekanger, M. Bhattacharya, and S. Horiuchi, (2024), 2405.17409.\n- [12] T. K. Gaisser, T. Stanev, and S. Tilav, Front. Phys. (Beijing) 8 , 748 (2013), 1303.3565.\n- [13] X.-J. Lv et al., (2024), 2403.11832.\n- [14] B. Bortolato, J. F. Kamenik, and M. Tammaro, Phys. Rev. D 108 , 022004 (2023), 2212.04760.\n- [15] Pierre Auger, E. W. Mayotte et al.\n- [16] T. P. A. Collaboration, Pierre auger observatory 2021 open data, 2021.\n- [17] D. Heck, J. Knapp, J. N. Capdevielle, G. Schatz, and T. Thouw, CORSIKA: a Monte Carlo code to simulate extensive air showers. (, 1998).\n- [18] T. K. Gaisser and A. M. Hillas, Reliability of the Method of Constant Intensity Cuts for Reconstructing the Average Development of Vertical Showers, in International Cosmic Ray Conference, , International Cosmic Ray Conference Vol. 8, p. 353, 1977.\n- [19] Pierre Auger, A. Aab et al., Phys. Rev. D 90 , 122005 (2014), 1409.4809.\n- [20] B. Efron and R. J. Tibshirani, An Introduction to the Bootstrap, Monographs on Statistics and Applied Probability No. 57 (Chapman & Hall/CRC, Boca Raton, Florida, USA, 1993).\n- [21] Pierre Auger, J. Stasielak, AugerPrime -The upgrade of the Pierre Auger Observatory, in 9th International Conference on New Frontiers in Physics,  \nuncertainties due to a lower number of free parameters when working with a shorter composition vector w motivate further optimization, which we leave for future work.", 'ACKNOWLEDGMENTS': 'B.B. and J.F.K. acknowledge the financial support from the Slovenian Research Agency (grant No. J1-3013 and research core funding No. P1-0035).  \n- 2021, 2110.09487.\n- [22] B. Bortolato, J. F. Kamenik, and M. Tammaro, Phys. Rev. D 109 , 043023 (2024), 2304.11197.\n- [23] J. Skilling, AIP Conference Proceedings 735 , 395 (2004), https://aip.scitation.org/doi/pdf/10.1063/1.1835238.\n- [24] J. Buchner, UltraNest github repository, "https:// github.com/JohannesBuchner/UltraNest , 2019.\n- [25] N. Arsene, (2021), 2109.03626.', 'Appendix A: Interpolation': 'In Fig. 6 we show the values of the mean (top) and variance (bottom) of the z 1 distributions as function of Z , for the [0 . 65 , 1] (left) and [1 , 2] (right) energy bins, respectively. The black dots, indicating the values extracted from simulations, can be fitted by two simple function of Z ,  \nh mean ( x ) = a m x + b m x 2 + c m + d m · log( x ) , h var ( x ) = a v x + b v √ x + c v + d v · log( x ) . (A1)  \nAs a consequence, we can extract the values of z 1 , and thus of the distances, for all primaries between Z = 1 and Z = 94 , without the need to simulate the longitudinal profile of the respective showers. The results of the fit are shown as red lines in Fig. 6.  \nFIG. 6: Fits of mean and variance of z 1 for [0 . 65 , 1] EeV (left) and [1 , 2] EeV (right) energy bins. The black dots represent simulated values, while the red line indicates the result of the fit.  \n<!-- image -->  \nAppendix B: [1,2] bin plots  \nTABLE II: Same as Table I, for the energy bin [1 , 2] EeV.  \n<!-- image -->  \nFIG. 7: Same as Fig. 5, for the energy bin [1 , 2] EeV.  \n<!-- image -->  \nFIG. 8: Same as Fig. 4, for the energy bin [1 , 2] EeV.  \n<!-- image -->'}