Depth of maximum of air-shower profiles: testing the compatibility of measurements performed at the Pierre Auger Observatory and the Telescope Array experiment

At the Pierre Auger Observatory and the Telescope Array, the measurements of depths of maximum of air-shower profiles, $X_\mathrm{max}$, are performed using direct observations of the longitudinal development of showers with the help of the fluorescence telescopes. Though the same detection technique is used at both installations, the straightforward comparison of the characteristics of the measured $X_\mathrm{max}$ distributions is not possible due to the different approaches to the analysis of the recorded events. In this work, the Auger-Telescope Array composition working group presents a technique to compare the $X_\mathrm{max}$ measurements from the Auger Observatory and the Telescope Array. Applying this technique the compatibility of the first two moments of the measured $X_\mathrm{max}$ distributions is qualitatively tested for energies $10^{18.2}~{\rm eV}<E<10^{19.0}~{\rm eV}$ using the recently published Telescope Array data from the Black Rock Mesa and Long Ridge fluorescence detector stations. For a quantitative comparison, simulations of air showers with EPOS-LHC, folded with effects of the Telescope Array detector, are required along with the inclusion in the analysis of the systematic uncertainties in the measurements of $X_\mathrm{max}$ and the energies of the events.


Introduction
The large amount of data on the ultra-high-energy cosmic rays, collected by the Pierre Auger Observatory [1] operating in the Southern Hemisphere since 2004 and the Telescope Array (TA) [2] operating in the Northern Hemisphere since 2008, opens unprecedented possibilities for performing a comparison of the spectra, arrival directions and mass compositions of the cosmic rays coming from the complementary regions of the sky. Any differences in the characteristics of the primary radiation found in such type of analysis would have important astrophysical implications and would help to clarify the origin of the most energetic particles in the Universe. For an unequivocal attribution of such differences to factors having an astrophysical nature, the comparison of the data of Auger and TA should rely on a deep understanding of the systematic effects present in the measurements of the two experiments which employ different kinds of the surface detectors (SD), different atmospheric monitoring equipment and programs, different designs of the fluorescence detectors (FD) and eventually different approaches to the data analysis. To address these issues a close collaboration of the people from various experiments is required, and to this scope, several working groups were created in preparation for the ultra-high-energy cosmic ray symposium that took place in Geneva in 2012 [3].
The activities of the mass composition working group in the recent years have been focused on the comparison of the Auger and TA data on the depth of the maximum of the extensive air showers measured with the use of the fluorescence telescopes. Until the present time, a good agreement between the two data sets was found as regarding the evolution of X max with energy [4], so regarding the compatibility of the shapes of the X max distributions [5] (the latter analysis was done for energies below 10 EeV). In 2018 the TA Collaboration has published new data [6] including results on the first two moments of the X max distributions with the larger statistics than ever before. In this paper we present a preliminary qualitative analysis of the compatibility of these TA measurements with the data of the Pierre Auger Observatory and outline the next steps required for finalizing the comparison.

Pierre Auger Observatory
At the Pierre Auger Observatory the longitudinal development of air showers is measured with the FD consisting of 24 fluorescence telescopes covering each 30 • in azimuth and 1.5 • − 30 • in elevation. The telescopes are grouped in units of six at four sites around the SD array of the area of 3000 km 2 . There is an additional site with three high elevation telescopes for detection of air showers with energies E 10 18 eV, but events with such energies will not be discussed in this paper.
An extensive program of monitoring of atmospheric parameters, important for a proper reconstruction of X max and energy of air showers, is run at the Observatory (see [1] for more details). Pressure, humidity and temperature at different atmospheric depths are determined on a 3-hour basis using data from the Global Data Assimilation System. Vertical aerosol optical depth (VAOD) is measured hourly with the FD using light profiles of the laser shots from two central laser facilities. The presence of clouds in the field of view of the FD telescopes is controlled each 15 minutes with the help of cloud cameras and elastic lidars, additional information on cloud coverage is obtained using laser shots from the laser facilities and data from the Geostationary Operational Environmental Satellites. For events of the extremely high energies or having unusual longitudinal profiles, a rapid monitoring of atmospheric conditions is performed using the lidars and the Photometric Robotic Atmospheric Monitor.
The reconstruction of the recorded events is performed using the Auger Offline software framework [7]. Only hybrid events, i.e. the events having at least one triggered SD station, that satisfy a number of quality selection criteria, are included in the X max analysis. Further, the high quality data pass an additional fiducial field-of-view selection to guarantee an unbiased acceptance of the showers almost independently on their X max (see Fig. 1). Finally, the data are corrected for the reconstruction, residual acceptance biases and resolution effects, and thus the resulting X max and σ(X max ) can be compared directly to the predictions from the ideal (without any detector) MC simulations (see [8] for full details). The most recent Auger data [9] come from the period 12/2004−12/2015 and contain about 10550 events in the energy range E > 10 18.2 eV discussed further.
In Fig. 2 one can see the evolution of the first two X max moments measured with the Auger Observatory. Starting from lower energies up to about 2 EeV, X max in the data is approaching the MC predictions for protons, thus the primary mass composition is becoming lighter. For energies above 2 EeV, the behavior of both X max and σ(X max )  is consistent with the increase of the average mass in the primary beam [10].

Telescope Array experiment
The TA data discussed in this paper were collected by the fluorescence telescopes installed at the Black Rock Mesa (BR) and Long Ridge (LR) sites [6]. There are 12 fluorescence telescopes at each of the sites covering the total field-of-view of 108 • in azimuth and 3 • − 33 • in elevation. These FD telescopes are placed in the southern part of the TA SD of area of 700 km 2 . The properties of the atmospheric profiles at the TA site are obtained from the data of the Global Data Assimilation System and, for estimating the systematic uncertainties in the measurements of X max , from the data collected with the NOAA National Weather Service radiosondes. The VAOD is measured each 30 minutes using laser shots from the central laser facility. The reconstruction of X max is performed with an average VAOD = 0.04 for all nights and the VAOD measurements are then used for estimation of the systematic uncertainties of X max and σ(X max ). The cloud coverage and cloud thickness are judged by eye and logged hourly by operators in the field.
The TA events used in the X max analysis are hybrid, as in the case of Auger. For an event to be accepted it is required the triggering of both the FD and of three SD counters adjacent to each other. The details on the reconstruction and selection of the events can be found elsewhere [6]. Differently from the Auger analysis, no fiducial field-of-view selection or corrections for the reconstruction and acceptance biases are applied to the TA data. Consequently, a larger number of events is kept in the data set, but the comparison of the data to the predictions of   the interaction models can be done only using the MC simulations processed through the same analysis chain as the data and including thus the analysis biases and the effects of the TA detector. The period of the data taking is 05/2008 − 11/2016, the data set consists of 3330 events with E > 10 18.2 eV. From the comparison of X max measured by the TA BR/LR to the predictions of the QGSJetII-04 model (Fig. 3) and from the analysis of the shapes of X max distributions [6,11], it follows that the data are compatible within the statistical and systematic errors to pure protons for all energies E > 10 18.2 eV. Due to the lower detector acceptance for deep high energy showers and due to the low number of the detected events, for E > 10 19.0 eV the shapes of the measured X max distributions become compatible also to the MC predictions for other pure primary compositions (helium, nitrogen) [6,11]. uncertainty 14.1 g cm as a conservative estimate. Other sources are added in quadrature, and we nd the total systematic uncertainty in max to be 17.4 g cm . The results are summarized in Table  The systematic uncertainties of ( ) max from the sources discussed above are also evaluated and given in Table  Adding in quadrature, we obtain 21.1 g cm As seen in Figure 14, within systematic uncertainties, max of the data is in agreement with QGSJet II-04 protons and helium for nearly all energy bins. There is clear separation between the region of systematic uncertainty and heavier elements such as nitrogen and iron. In the last two energy bins there is some overlap between the systematic uncertainty region of the data and the nitrogen, but statistics in the data there are very poor. Care must be taken in interpreting Figure 14, since max by itself is not a robust enough measure to fully draw conclusions about UHECR composition. When comparing max of data to Monte Carlo, in addition to detector resolution and systematic uncertainties in the data that may hinder resolving the different elements with relatively similar masses, the issue of systematic uncertainties in the hadronic model used to generate the Monte Carlo must also be recognized. This will be discussed in Section . Referring back to Figures 12 and 13, we can see that although the max of the data in Figure 14 lies close to QGSJet II-04 helium, the ( ) max of the data is larger than the helium model allows for energy bins with good data statistics. For this reason, we will test the agreement of data and Monte Carlo not just by comparing max and ( ) max , but by using the entire distributions. The elongation rate of the data in Figure 14 found by performing a t to the data is found to be 56.8 5.3 g cm decade . The dof of this t is 10.67 9. Table summarizes the observed rst and second moments of TA s observed max for all energy bins.

Method
If one wishes to draw conclusions about agreement between the data and the models, we should employ a test that measures

Remarks on the Auger and TA X max measurements
For energies above ∼ 10 18.5 eV the first two moments, the shape of X max distributions and the elongation rate of X max measured by Auger are incompatible to the MC predictions for protons [8,9] (the same conclusions have been obtained from the measurements involving the information from the SD [12,13]). The interpretation stating that the TA X max measurements are compatible to the predictions for QGSJetII-04 protons [6,11] does not necessarily contradict the Auger results. In case of the TA BR/LR data, only the compatibility to the pure beams was tested, thus it can not be excluded that the TA data can be described well by the mixed compositions of QGSJetII-04 (as was already shown in [5]) and of other interaction models. It should be noted as well that the statistical significance of such a comparison will be lower than in case of Auger due to the difference in the sizes of the data sets (Fig. 4) reflecting the difference in hybrid exposures (sizes of the SD arrays, data taking periods and detection/selection efficiencies) of the Auger and TA experiments.
Instead of comparing the interpretations, the primary aim of the composition working group is the comparison of the X max measurements themselves, including the check of the compatibility of the first two X max moments, of the shapes of the X max distributions and of the X max elongation rates. Further we present the technique for doing such a comparison, perform a preliminary qualitative comparison of X max and σ(X max ) and outline the next steps needed for a completion of the full quantitative evaluation of the compatibility of the Auger and TA X max measurements.  . Numbers of events (top) and their ratio (bottom) in the X max data sets of Auger [9] and TA BR/LR [6] for the common energy range E > 10 18.2 eV. The energy binning is the same as used by the TA [6]. Data are binned using the energy scales of each of the experiments.

A method of comparison of the Auger and TA X max measurements
As explained in the previous section, the Auger and TA X max measurements cannot be compared directly because of the use of different approaches to the analysis of the recorded events. In the case of Auger, due to the application of the fiducial field-of-view selection and consequent correction of the biases, the X max moments are free from the detector effects and can be compared directly to the ideal MC simulations. The data of the TA on the other hand include all detector and analysis biases.
To compare the Auger and TA data one should convert Auger X max values into the values folded with the TA detector effects (Auger X max ⊗ TA). To do this, it was proposed [14] to find the MC compositions giving the best description of the shapes of the Auger X max distributions and to process these compositions through the full simulations and analysis chain used for the TA BR/LR data. The MC compositions describing the Auger X max distributions were found in [9,15] and will be referred hereafter as AugerMix.
In Fig. 5 an example of AugerMix compositions for QGSJetII-04 and EPOS-LHC in comparison to the Auger data is shown for the energy bin 10 19.0 eV < E < 10 19.1 eV. One can see that the MC template for EPOS-LHC provides a better description of the data compared to the template for QGSJetII-04. As demonstrated in [9], the p-values characterizing the agreement between Auger-Mix of QGSJetII-04 and the Auger data are ∼ 0.01 in all energy bins between 1 and 10 EeV. The reason for a non-satisfactory description of the data might be that the AugerMix of QGSJetII-04 is composed only of protons and helium and the widths of the X max distributions for such a mix of light elements are too wide compared to the widths found in data [9]   This fact is illustrated in Fig. 6 where one can see that σ(X max ) for AugerMix of QGSJetII-04 is significantly larger compared to the values found in the data. Since the AugerMix compositions of QGSJetII-04 do not reproduce correctly the shapes of the X max distributions measured by Auger, they are not the optimal choice for performing the comparison of the Auger and TA measurements. For the proper comparison, one can use the mass composition templates found using EPOS-LHC: as shown in Fig. 6, with this interaction model both measured X max and σ(X max ) are described well, as also indicated by the good p-values found in [9] testing the agreement between the X max distributions.
At the moment of writing of this paper the only available simulations for the TA BR/LR detector were done with QGSJetII-04 while EPOS-LHC simulations have not been finished yet 2 . Thus for a preliminary comparison of the X max moments, one can only resort to a re-weighting procedure first introduced and thoroughly tested in [4]. The re-weighting allows one to transform X max distributions for the AugerMix of QGSJetII-04 into the X max distributions of the AugerMix of EPOS-LHC taking all re-  construction, selection biases and TA detector effects into account. Technically this is achieved by the re-weighting of each reconstructed event of QGSJetII-04 by a factor equal to the ratio of the X max p.d.f.'s AugerMix(EPOS-LHC)/AugerMix(QGSJetII-04) estimated at the true (generated) X max and energy of this event.
For the comparison presented here the re-weighting was performed using the AugerMix compositions from [8] for the energy range 10 18.2 eV < E < 10 19.0 eV. From Fig. 7 one can see that X max ⊗ TA are ≈ 5 g cm −2 smaller than the generated values. This bias can be caused by the relatively smaller detector acceptance for the showers in the deeper tails of the X max distributions. Similarly, σ(X max ) ⊗ TA is also a few g cm −2 smaller than in the ideal MC.
In Fig. 8 the comparison of AugerMix ⊗ TA for EPOS-LHC (representing the Auger data transferred into the TA detector) to the TA BR/LR measurements is shown. One can see that X max of the two experiments agree within the statistical and systematic errors with mostly shallower X max values (heavier mass) for TA compared to Auger and that the agreement is getting better for E > 10 18.5 eV. In the past (e.g. [16]) the agreement of the Auger data on X max to the Middle Drum TA measurements was tested, and it was found that both measurements are compatible with an average difference of ∆ = (2.9 ± 2.7 (stat.) ± 18 (syst.)) g cm −2 and no apparent dependence on energy (see Fig. 9). The Middle Drum has different hardware and is located 4 − 5 km farther away from the border of the SD compared to the BR/LR stations and might have worse efficiency for low energy (E < 10 18.5 eV) showers [6].
The energy dependence in the comparison of the Auger and TA data can be naturally tested using the anal-  ysis of the measured X max elongation rates. The elongation rate of Auger is measured to be (80±1) g cm −2 /decade between 10 17.2 and 10 18.3 eV, above 10 18.3 eV it flattens to 18 Figure 9. Comparison of the Auger [8] and TA Middle Drum [17] measurements of X max . The plot is taken from [16].
(26 ± 2) g cm −2 /decade [9] indicating the increase of the primary mass. As noted previously, the TA elongation rate is folded with the detector effects and thus, due to the detector acceptance, it can be flatter than expected from the ideal MC. Given the acceptance biases, the statistical and systematic uncertainties on X max , it is premature to make interpretations claiming an evolution in composition from the TA BR/LR data. If X max is fit for all data observed, the slope is found to be 53 ± 3 g cm −2 /decade, but with χ 2 /ndf of 24/9 (p-value = 4 × 10 −3 or 2.7σ significance) (see also [18]). If the lower energy of the fit is adjusted to E > 10 18.4 eV, the slope is found to be 40±5 g cm −2 /decade with an acceptable χ 2 /ndf of 5.4/7. Model predictions of the elongation rate, folded with the BR/LR detector effects, for a single species is ≈ 50 g cm −2 /decade. Regarding the comparison of σ(X max ) (Fig. 8) up to 10 18.7 eV, there is a good agreement between the Auger and TA data. Above this energy, the statistical fluctuations in the TA data become quite large making it more difficult to draw any firm conclusions on the compatibility of the Auger and TA measurements.

Discussion
In this paper the Auger -TA working group on the mass composition has presented a method for comparing X max data of the two experiments and the preliminary qualitative checks of the compatibility of X max and σ(X max ) for energies 10 18.2 eV < E < 10 19.0 eV using the recently published data from the TA BR/LR stations. While we could not identify any discrepancies beyond the statistical and systematic errors, more definite conclusions can be obtained after the shower simulations with EPOS-LHC for the TA BR/LR will be available. Other factors which will need to be taken into account are the difference between the Auger and TA energy scales [19] and systematic errors on the fractions of nuclei in the AugerMix MC templates used by Auger to fit the measured X max distributions. The eventual quantitative comparison will include the tests of the compatibility of the first two X max moments, of the X max distributions and of the X max elongation rates for the whole energy range E > 10 18.2 eV.