Testing the Compatibility of the Depth of the Shower Maximum Measurements performed at Telescope Array and the Pierre Auger Observatory Auger-TA Mass Composition Working Group Report

. The Telescope Array and the Pierre Auger Observatory estimate the composition of ultra-high-energy cosmic rays by observing the distribution of depths of air-shower maxima, X max . Both experiments directly observe the longitudinal development of air showers using ﬂuorescence telescopes with surface particle detectors used in conjunction to provide precision in determining air-shower geometry. The two experiments di ﬀ er in the details of the analysis of events, so a direct comparison of X max distributions is not possible. The Auger – Telescope Array Composition Working Group presents their results from a technique to compare X max measurements from Auger with those of Telescope Array. In particular, the compatibility of the ﬁrst two moments of the X max distributions of Auger with the data from the Black Rock Mesa and Long Ridge detectors of the Telescope Array is tested for energies above 10 18 . 2 eV. Quantitative comparisons are obtained using air-shower simulations of four representative species made using the Sibyll 2.3d high-energy interaction model. These are weighted to ﬁt the fractional composition seen in Auger data and reconstructed using the Telescope Array detector response and analysis methods.


Context and History
The Telescope Array (TA) observatory [1] and the Pierre Auger Observatory [2] (Auger) both measure the nuclear composition of ultra-high-energy cosmic rays (UHECRs) by observing the distribution of the depths of shower maxima (X max ) of extensive air showers reconstructed in hybrid mode, combining surface detector (SD) and fluorescence detector (FD) measurements. The conclusions and interpretations made from these measurements have differed between the two observatories, which has lead to confusion amongst outside observers. One source of this confusion may be that TA and Auger employ different strategies in selecting the data sets used in the measurements. Auger selects events with analysis cuts designed to minimize bias in X max acceptance and reconstruction, and then corrects the resulting X max moments for remaining biases. TA on the other hand selects all well-reconstructed events and then models the biases using Monte Carlo (MC) * e-mail: bergman@physics.utah.edu * * e-mail: auger_spokespersons@fnal.gov simulations. This difference in event selection strategy makes direct comparison of X max distributions problematic.
Beginning with UHECR2012 [3], the Auger-TA Mass Composition Working Group has tried to assess the degree of agreement between Auger and TA X max measurements. The procedure for comparison involves making a representative breakdown of the Auger X max distributions into several single species fractions using a given high-energy interaction model, then using the representative fractions and high-energy interaction model in the TA MC simulation. The result is interpreted as the X max distributions TA would have measured given the X max distributions actually measured at Auger. These distributions, specifically their first and second moments, are compared to the actual TA measurements to assess the degree of agreement. In practice, the four representative nuclear species are taken to be hydrogen (H), helium (He), nitrogen (N), and iron (Fe). The relative fractions of these species in energy bins found to represent the Auger data is called the AugerMix.
At UHECR2014 [4], a comparison of the AugerMix to TA Middle Drum hybrid data was made, just comparing X max . This showed agreement within systematic uncertainties. At UHECR2016 [5] and UHECR2018 [6], the AugerMixes using both QGSJetII-04 [7] and EPOS-LHC [8] high-energy interaction models were compared to TA Black Rock and Long Ridge hybrid data using both X max and σ(X max ). The comparison showed agreement within systematic uncertainties for X max but tension at energies above 10 18.8 eV in σ(X max ). Due to problems in the implementation of EPOS-LHC within CORSIKA [9] at that time, the EPOS-LHC result was actually performed by re-weighting QGSJetII-04 distributions.
While QGSJetII-04 works well in simulating the response of the TA SD (and FD) to air showers, it cannot describe Auger measurements because it predicts broader distributions than observed. It is thus useful to use another high-energy interaction model. We show here the result of using the Sibyll 2.3d [10] high-energy interaction model in both making the AugerMix and in simulating the result of that mix in the TA detector.

The AugerMix with Sibyll 2.3d
Auger hybrid X max data [11] were binned by energy into X max histograms. The energy bins were chosen to match those used by TA, which are constrained by the TA statistics. The bin edges in log 10 4, 19.9). In each of the energy bins, X max distributions from Conex [12] simulations of single species (H, He, N, Fe) using Sibyll 2.3d were created. These single-species distributions are then modified by the X max acceptance and resolution following [13] and used as templates in fitting the observed Auger distribution into a set of fractions. While there are correlations in the fit fractions, only the best-fit fractions are used in this analysis.
The fit AugerMix fractions using Sibyll 2.3d as a function of energy are shown in Figure 1. These fractions describe the Auger measured data as the weighted sum of the templates can reproduce both X max and σ(X max ) as shown in Figure 2.

Generating the TA distributions
The AugerMix fractions were used to generate X max distributions as they would have been seen in the TA detector. Single-species MC sets for each of H, He, N, and Fe were generated using CORSIKA v7.7402 with Sibyll 2.3d as the high-energy interaction model. In total, 250 events for each species in each tenth-of-a-decade bin from 10 18 eV to 10 20.6 eV were generated using CORSIKA/Sibyll with optimal thinning [14]. The events were dethinned using the Stokes method [15] and the response of TA SDs placed every 6 m in a grid around the shower core was calculated using GEANT [16] to create a "tile file" for each simulated event. These events were then sampled at random locations and with random azimuthal orientations about the TA detector area with sufficient sampling to produce   the measured flux. The response of the TA FD was also simulated at this time from the stored longitudinal shower information from CORSIKA. This produces a spectral hybrid composition data set for each individual species. This data is in the same format as the collected TA hybrid data, and the same program of analysis cuts is applied to the simulated data that is used with the measured data.
The X max distributions resulting from this analysis are combined using the AugerMix fractions as weights. As each single-species set is a spectral set using the same flux, this is the same as selecting events according to the Auger-Mix at the "thrown" level of MC.

Comparison of the AugerMix to TA Data
The results of the comparison between TA and the Auger⊗TA results using the Sibyll 2.3d AugerMix are shown in Figure 3. The X max comparison shows agreements within systematic uncertainties throughout the compared energy range, with some slight tension at the lowest energies. However, σ(X max ) shows some disagreement between TA and Auger⊗TA in the upper half of the 10 18 eV decade; the TA results show a roughly constant width throughout the decade, while the Auger⊗TA results show a continual narrowing of the distribution.

EPJ Web of Conferences
, 02008 (2023) https://doi.org/10.1051/epjconf/202328302008 283 UHECR 2022 Figure 3. Comparison of X max (left) and σ(X max ) (right) between measured TA data distributions (blue squares) and the Auger⊗TA results (red circles) using Sibyll 2.3d.  To compare the shape of the distributions directly, we perform an Anderson-Darling test [17] on the two distributions in each energy bin. The test is made after removing the difference in means between the two distributions. In Figure 4 we plot the p-values of the test as a function of the energy for both the Auger⊗TA comparison to TA, and for the AugerMix comparison to the Auger data. Examples of distributions that match and that do not match are shown in Figure 5  One systematic effect that was not considered in these comparisons, was the variation in atmospheric aerosol levels on a night-by-night basis. The MC simulations use a constant vertical aerosol optical depth (VAOD) value of 0.04, whereas the data was collected with varying VAODs values, which have a mean of 0.04. The systematic effect of this variation may add up to 18.9 g/cm 2 to the measured width of the X max distributions. If we smear the Auger⊗TA distributions by and additional 18.9 g/cm 2 and redo the Anderson-Darling comparisons as before, we get a set of p-values consistent with both distributions coming from the same parent distribution. The values from this comparison are shown in Figure 6. It should be noted that this additional smearing of the X max distributions does not significantly ameliorate the differences in widths seen in Figure 3. The additional 18.9 g/cm 2 will be added in quadrature to the significantly larger intrinsic widths of the Auger⊗TA data.

Discussion
We have made a comparison of the cosmic ray composition measured by Auger as it would have been seen at TA. This was performed using a set of representative fractions of single species X max distributions generated using Sibyll 2.3d to simulate the response of the TA detector, also using Sibyll 2.3d, and compare the resulting distribution to the actual data distributions measured in TA. The result shows that the Auger⊗TA X max measurements agree within systematic uncertainties. However, there is disagreement in the width σ(X max ) at energies between 10 18.5 and 10 19 eV. This difference in widths is similar to that seen previously in Auger⊗TA comparisons preformed using both QGSJetII-04 and EPOS-LHC models. The difference in widths can be slightly ameliorated, but not removed, by smearing the Auger⊗TA distributions by an additional factor of 18.9 g/cm 2 to account for variations in atmospheric conditions which effect TA data but are not present in these simulations.