Tests of hadronic interactions with measurements by Pierre Auger Observatory

The hybrid design of the Pierre Auger Observatory allows for the measurement of a number of properties of extensive air showers initiated by ultra-high energy cosmic rays. By comparing these measurements to predictions from air shower simulations, it is possible to both infer the cosmic ray's mass composition and test hadronic interactions beyond the energies reached by accelerators. In this paper, we will present a compilation of results of air shower measurements by Pierre Auger Observatory which are sensitive to the properties of hadronic interactions and can be used to constrain the hadronic interaction models. The inconsistencies found between the interpretation of different observables with regard to primary composition and between their measurements and simulations show that none of the currently used hadronic interaction models can provide a proper description of air showers and, in particular, of the muon production.


Introduction
The Pierre Auger Observatory [1] is designed to measure ultra-high energy cosmic rays (E > 10 17 eV) through the detection of extensive air showers (EAS). The two main parts of the observatory are the Fluorescence Detector (FD), composed of 27 telescopes grouped in 4 sites and the Surface Detector (SD), which consists of an array with 1660 water-Cherenkov detectors (WCDs) distributed over 3000 km 2 . The FD measurements allow the reconstruction of the longitudinal profile of EAS, whereas the SD can measure the spatial distribution of particles on the ground and their arrival times. The information obtained from both detectors can be used to reconstruct a number of EAS observables. By comparing the measurements of these observables to predictions from Monte Carlo simulations using hadronic interaction models, one can either infer the cosmic-ray's mass composition or study the properties of hadronic interactions by testing the models. For each one of these goals, a different set of observables turns out to be more suitable.
The most reliable observable for mass composition inference is X max , the depth of maximum energy deposit. Lighter primaries penetrate deeper into the atmosphere than heavier ones, producing on average showers with larger X max values. It is well known that the X max is mostly driven by the properties of only the first interaction, which starts the dominant electromagnetic cascade producing most of the electron and positron contribution to the energy deposit profile. Although the exact relation between X max , the primary energy, and mass depends on * e-mail: raul.prado@desy.de * * Full author list: http://www.auger.org/archive/authors_2018_05.html the properties of hadronic interactions, the theoretical uncertainties on the predictions of X max by the most recent hadronic models are relatively small when compared to other observables (e.g. number of muons). The difference between the model predictions of X max is ≈ ±20 g/cm 2 (which translates into a difference in ln A of ≈ ±0.8) and is constant over energy [2]. Because of this, the cosmic-ray composition derived from the X max measurements is commonly used as a reference, allowing the hadronic models to be tested by comparing the composition interpretation from other observables with that using X max . The X max measurements by Auger can be found in Ref. [3,4].
The muonic component of EASs is particularly interesting in the context of testing hadronic interactions. Muons are mostly produced by the decay of charged mesons when their energies are low enough such that their decay length is comparable to the interaction length. Before reaching such low energies, a number of generations of hadron-air interactions is required. As a consequence, the muonic component is sensitive to the chain of hadronic interactions and, in particular, to their particle production properties. It is well known that the current hadronic models cannot provide a satisfactory description of the muon production in EAS. The most popular manifestation of this fact is the so-called muon deficit problem, i.e. the number of muons (N µ ) predicted by EAS simulations is significantly smaller than that in measurements. This deficit has been observed by several experiments [5][6][7][8][9][10]. Therefore, experimentally accessing the muonic component of EAS, through the measurements of N µ as well as further observables sensitive to muons, is of extreme value in constraining hadronic models and possibly infer properties of hadronic interactions.   Figure 1. Results of the top-down analysis [8]. See text for details.
In this paper, we will present a selected compilation of measurements by the Pierre Auger Observatory which can be used to test hadronic interaction models. In Sec. 2, we present the results of the so-called Top-Down analysis, in which a set of hybrid events are compared in detail with simulations in a way that it is possible to disentangle the contributions from the electromagnetic and muonic components of EAS. In Secs. 3 and 4, we present the results of two analyses based on EAS observables which are purely muonic: the number of muons measured in highly inclined events (Sec. 3) and the muon production depth (Sec. 4). In Secs. 5 and 6, two analyses based on the parameter risetime (t 1/2 ) are presented: the measurements of the azimuthal asymmetry of risetime(Sec. 5) and the so-called ∆ s -method (Sec. 6). Finally, in Sec. 7 we summarize the results.

Top-Down analysis
By using the experimental setup of Auger, the muon content of EAS can be directly measured in highly inclined events (θ > 60 • ) by using the atmospheric attenuation to eliminate the electromagnetic component (see Sec. 3). However, it is not possible to isolate the muonic from the electromagnetic component at the detector level for vertical events (θ < 60 • ). To overcome this difficulty, an original analysis procedure was developed aiming to evaluate the muon content in hybrid events by means of a detailed comparison with Monte Carlo simulations. The full description of the analysis and the results can be found in Ref. [8].
For a set of 411 high quality hybrid events with 10 18.8 < E/eV < 10 19.2 , the first analysis step was to find simulated events in which the longitudinal profile matches those measured by the FD. By construction, these simulated events have an X max and energy compatible with the corresponding measured event. The above procedure allows us to compare the SD signals from measured and simulated events without having to account for differences in the longitudinal development of the showers. The first step of this comparison can be seen in Fig. 1 (left), where we show the average ratio R between the shower size parameter S 1000 obtained from data and simulations as a function of the shower zenith angle sec(θ). Apart from two different hadronic models, Epos LHC [11] and QGSJet II-04 [12], we also show two different composition assumptions: proton only and a mixed composition scenario ob-tained from the interpretation of the X max measurements by Auger. It can be seen that for any model and composition scenario, the ratio R is greater than ≈ 1.2, which shows that the ground signal of measured events is at least 20% larger than the simulated ones with the same longitudinal development. Furthermore, an evolution of R with zenith angle is observed, which will be important for the next step of the analysis.
In a second step, the particle distributions in simulated events were rescaled in a way to make the data and simulations compatible in terms of the shower size S 1000 . Given that the different components of the shower evolve differently with zenith angle and a large zenith angle range is covered by the present analysis, it becomes possible to separately rescale the electromagnetic component, represented by R E , and the hadronic component, represented by R had . While R E is directly related to the shower energy scale, R had is responsible for scaling the muon content, meaning that by increasing R had the N µ in the shower would also increase by the same factor. A model derived from simulations was used to perform the rescaling and the results can be seen in Fig. 1 (right). To make the ground signal of simulations compatible with those observed in the data, an increase from 30% to 70% of R had is required, considering all the possible combinations of hadronic models and composition assumptions. On the other hand, it is observed that the energy scale factor is never required to increase by more than 10% and the case without rescaling (R E ) is covered by the systematic uncertainties in all cases.
The results described above represent a manifestation of the muon deficit problem. The best case scenario of the mixed composition assumption with Epos LHC as hadronic model still requires an increase of 30% on the N µ predicted by simulations at energies ≈ 10 19 eV. For the first time, the muon deficit was evaluated by disentangling the contribution of the muonic component from the energy scale.

Number of muons in highly inclined showers
An efficient way of measuring muons using the SD without contamination from the electromagnetic component is by using the atmosphere as a shield. Above a certain atmospheric depth (∼ 2000 g/cm 2 ), the electromagnetic component is strongly attenuated while most muons still penetrate down to the ground. Such large atmospheric depths are reached by highly inclined events. In this analysis, events with zenith angle 62 • < θ < 80 • measured by both SD and FD are selected and their muon content is reconstructed. The full analysis description and results can be found in Ref. [7]. Because of the complicated dependencies of the muon density (ρ µ ) with the zenith angle (θ) and azimuthal angle (φ) of an event and the axis distance of a station (r), a template fit method was adopted to describe the distribution of muons on the ground of measured events. The templates were built using Monte Carlo simulations of proton primaries at 10 19 eV using QGSJet II-03 [12] as the hadronic model. The measured SD signals were then fitted by using the templates in which the normalization parameter (N 19 ) was left free. By applying the same fit to simulated events, the bias on N 19 was estimated. After correcting for the bias, an unbiased estimator for the template normalization R µ was obtained. For any given zenith angle, N µ can be recovered from R µ . For example, R µ = 1 corresponds to 1.2 × 10 7 muons at the ground level above 0.3 GeV (muon energy threshold in the WCDs) for θ = 70 • .
The average R µ over energy bins as a function of energy is shown in Fig. 2 (left), where the energy was reconstructed independently of R µ from the FD measurements. The predictions for proton and iron primaries using two hadronic models are also shown for comparison. The systematic uncertainties, depicted with square brackets, are dominated by the uncertainties on the energy scale.
By comparing our measurements with predictions of simulations, we can observe that, even considering the lowest values of R µ allowed by the systematic uncertainties, the muon content of inclined events would imply a very heavy composition interpretation. The compatibility of this interpretation with the composition scenario expected from the X max measurements is tested in Fig. 2 (right), where we show ln R µ versus X max at 10 19 eV. One can conclude that, by assuming the X max composition, the measured value of R µ is significantly larger than that expected from the simulation, regardless of the hadronic model. This fact shows another manifestation of the muon-deficit problem, which is complementary to the result presented in Sec. 2 in terms of the event zenith angle.

Muon production depth
The time structure as measured by the SD contains a lot of information about the EAS development. The analyses presented in the next three sections make use of the time traces recorded by the WCDs to reconstruct EAS properties.
Since a muon propagates nearly linearly from its production point down to the ground, its arrival time can be mapped into their production depths by using simple geometrical considerations. In practice, this procedure requires SD stations with a large muon purity and a precise knowledge of the muon time delay. The former requirement is satisfied by selecting stations far from the shower axis in highly inclined showers, while a model based on simulations is used to correct for the muon time-delay.
After reconstructing the muon production depth for each event, the maximum of the profile (X µ max ) can be determined, the energy evolution of its moments can be evaluated, and compared to predictions by simulations. The X µ max carries the information about the depth of the first interaction, folded with the longitudinal development of the hadronic component of the EAS which is responsible for the muon production. It has been shown that X µ max is very sensitive to the properties of pion-air interactions [14,15].
A first version of the muon production depth analysis was published in Ref. [16], and recently, an updated version has been released in Ref. [13]. Among other improvements, the updated analysis extends the zenith angle and the lateral distance range of the selected stations, increasing significantly the number of available events. With the resulting larger statistical power, it is possible to derive not only the X µ max , but also the σ[X µ max ]. In Fig. 3, we show both X µ max and σ[X µ max ] as a function of energy in comparison with the prediction from simulations of proton and iron primaries. Similarly to what can be seen in the results published in Ref. [16], the composition interpretation of the X µ max implies the presence of primaries as heavy as iron for QGSJet II-04 [17] and even heavier for Epos LHC. This interpretation is clearly inconsistent with that obtained from the X max measurements, which means that the hadronic models cannot properly describe the longitudinal development of muon production yet. December 2016 have been used in this analysis. Considering the applicability ranges of this work and the selection criteria described in Sec. 2, the number of UHECR events here analyzed is 2227. Data have been studied as a function of the primary energy. A bin width ∆log 10 (E/eV)=0.1 is chosen for energies log 10 (E/eV) between 19.2 and 19.8. Not having enough statistics to keep the same binning, data are integrated in one bin for log 10 (E/eV) in the range [19.8-20]. For each energy bin, the first two moments of the X * µ max distribution are evaluated on data and are compared directly to the expectations obtained from Monte Carlo simulations after the reconstruction procedure (Sec. 2). We note that data and Monte Carlo are both equally biased by the reconstruction, so the relative distance to the reference lines does not vary in X * µ max (see below conversion to the mean logarithmic mass ln A ) and no systematics are associated to these effects. On the contrary, the physical X µ max would display the mass and model spread as systematics, as discussed previously. The overall systematic error on the first two moments of the X * µ max distribution turns out to be around 4 g/cm 2 and 3 g/cm 2 respectively, and due to two sources: the small dependence of the selection efficiency of the quality cuts on the primary mass ( 1 g/cm 2 ) and the time variability of data. An additional systematic error of 7.5 g/cm 2 can be associated with the event selection and procedure to fit the MPD profiles and needs to be taken into account in the determination of ln A (see below). The results on X * µ max are shown in Fig. 3 (left) by black squares, with their statistical (black line) and systematic uncertainties (gray band). For each energy bin, the number of events is indicated. From the comparison with the predictions, the inconsistency among models and data is evident. In the case of EPOS-LHC, data are at odds with predictions for all reasonable masses, in the whole energy range. Considering instead QGSJetII-04 and in particular iron expectations, a mild incompatibility arises at the highest energies. We have also checked that when converting X * µ max to X µ max by using the reconstruction bias (Fig. 2 left) averaged on mass/models, we obtain a good agreement with the results shown in Fig. 3 and with results presented in [8]. The inconsistencies outlined here make it difficult to draw firm conclusions on composition with our measurement of X * µ max : we see that the predictions of X * µ max from the two hadronic models are significantly different in absolute value (≈ 35 g/cm 2 ). However, we can note that the muonic elongation rate, i.e. the rate of change of X * µ max with primary energy, is predicted to be about 6 Figure 3. Results of the muon production depth analysis [13].

Azimuthal asymmetry of the risetime
Both analyses presented in this and in the next section are based on the risetime parameter (t 1/2 ). It is defined as the time of increase from 10% to 50% of the integrated signal of each WCD station. Given that muons on average arrive earlier at the detector than electrons, t 1/2 turns out to be sensitive to the relationship between these two components. This implies that, for example, increasing the relative number of muons arriving in a given station would result in a faster increase in the WCD signal and thus a smaller t 1/2 . Because of this, the t 1/2 parameter is suitable for composition inferences and tests of hadronic models.
Because t 1/2 is a station quantity with a very complicated dependence on parameters like shower zenith-angle, station azimuthal-angle and station lateral-distance, it is impractical to define one unique t 1/2 -based parameter for each event. Instead, the two analyses presented in this paper follow approaches in which these above mentioned dependencies are explored in order to obtain information about the EAS development. In this section, the azimuthal dependence of t 1/2 will be explored.
It is observed that the average t 1/2 as a function of the station azimuthal angle (ζ) presents a maximum at the point in which the shower front reaches the ground the earliest, where the relative contribution of the electromagnetic component to the signal is the largest. This point is defined as ζ = 0. The shape of t 1/2 as a function of ζ can be well de-scribed by t 1/2 /r = a + b cos ζ + c cos 2 ζ, where the factor r is included to account for the almost linear dependence of t 1/2 on the distance r from the shower. The degree of asymmetry of t 1/2 can be quantified by the asymmetry factor defined as b/(a + c). For larger asymmetry factors, the peak of t 1/2 at ζ = 0 is more pronounced.
In nearly vertical events, the asymmetry factor is expected to vanish for symmetry reasons. The same behavior is expected for very inclined showers, in which the shower front is dominated by muons that are very weakly attenuated and consequently present very similar time structures for all azimuthal regions. In between these two extremes, the asymmetry factor presents a maximum. It is observed that the zenith angle which corresponds to this maximum (θ max ) carries information about the relative amount of muons in EAS and consequently can be used as an observable sensitive to composition and hadronic interaction properties. The full description of the analysis and results can be found in Ref. [18].
In Fig. 4, we show the energy evolution of sec θ max derived in energy bins. The predictions from simulations are also shown for two hadronic interaction models for proton and iron primaries. The left (right) plot shows the results for the sub-dataset of stations with 500 < r/m < 1000 (1000 < r/m < 2000). The composition interpretation of the sec θ max measurements imply, in general, the presence of heavier primaries on average, when compared to the composition expected from the X max measurements (see Ref. [18]). The difference is of order of ∆ ln A ≈ 1. The only case in which both observables are consistent is for QGSJet II-04 at 1000 < r/m < 2000. The inconsistency with the X max composition and between the two different axis distance ranges, for QGSJet II-04, are indications of a problem with the description of the t 1/2 asymmetry by the models. Since the t 1/2 is sensitive to the muon content of EAS, these problems can be related to the muon deficit.

The Delta method
The second analysis based on risetime to be presented in this paper makes use of the dependence of t 1/2 on the distance r from the shower axis to derive a new event parameter ∆ s . From the t 1/2 measured for each SD station within an event, the parameter ∆ i = (t 1/2 − t bench 1/2 )/σ t 1/2 is computed, where t bench 1/2 is obtained from a benchmark and σ t 1/2 is the uncertainty on t 1/2 . ∆ s is defined as the average of ∆ i over N selected stations, ∆ s = i ∆ i /N.
Two intermediate steps are important for the derivation of ∆ s : the definition of the benchmarks and the estimation of the t 1/2 uncertainties. First, the benchmark is defined as the average behavior of t 1/2 as a function of r observed in data for a given reference primary energy. By construction, the ∆ s vanishes at the reference energy. The benchmark curve was parametrized accounting for its dependence on the zenith angle. Secondly, the t 1/2 uncertainties are determined through a detailed study of the subsets of SD twins stations, which are pairs of stations located effectively at the same position, and pairs, which are stations within the same event at different positions but at approximately the same distance from the shower axis.
A detailed parametrization of σ t 1/2 was then performed by accounting for its dependences on r, θ, and the station signal S . The full description of the method and the results can be found in Ref. [19].
In Fig. 5, we show the ∆ s in energy bins as a function of the primary energy together with the predictions by simulations with two hadronic models. The two plots show the results for the two SD arrays, with SD stations spaced by 750 m on the left and 1500 m on the right. The energy bins chosen for the definition of the benchmarks are 17.7 < lg(E/eV) < 17.8 and 19.1 < lg(E/eV) < 19.2 for the 750 and 1500 m array, respectively. The composition interpretation from the comparison of the measured ∆ s with the predictions from simulations is not consistent with the one from the X max measurements (see Ref. [19]).
Although strong similarities can be seen on the energy evolution of ln A , a systematic difference of ∆ ln A ≈ 1 to 1.5 is observed, where the composition derived from ∆ s is heavier than the one from X max . These differences on ln A are of the same order as the differences observed in the t 1/2 azimuthal asymmetry analysis (see Sec. 5). Thus, the ∆ s method provides us with further evidence that the hadronic models cannot properly describe the t 1/2 parameter, assuming that the uncertainties on the X max predictions are relatively small. Again, this can be related to the muon deficit problem, since the t 1/2 is sensitive to both electromagnetic and muonic components.
It is worthwhile noting that, apart from testing the predictions of the hadronic models, the ∆ s method can also be effectively used for composition inferences through a cross calibration procedure with the X max measured by the FD (see Ref. [19] for the composition results).

Summary
The experimental setup of the Pierre Auger Observatory offers great capability for testing hadronic interactions. While the calorimetric energy and the longitudinal development of the showers (X max ) can be measured by the FD, the SD provides the lateral distribution of the particles on the ground as well as the particle arrival time distributions. By combining measurements of both detectors, it is possible to reconstruct a number of properties of EAS that can be used to test the consistency of the EAS description with the simulations with hadronic interaction models.
We have presented here a selection of five analyses. In Secs. 2 and 3, the muon content of EAS was measured with two different approaches and using datasets with complementary zenith angle ranges. The muon deficit in simulations was observed in both cases. An increase of at least 30% on the number of muons in simulation predictions is required to make the muon content measurements consistent with the composition interpretation from the X max measurements.
In Sec. 4, the measured arrival time distributions of the particles at the ground were used to access the longitudinal development of the muonic component. The maximum of the muon production depth was reconstructed and, by comparing its moments to simulation predictions, it is observed that the hadronic models cannot provide a consistent description of the muon production profile and the X max .
In Secs. 5 and 6, two analyses based on the t 1/2 parameter were presented: the t 1/2 azimuthal asymmetry analysis and the ∆ s method. The results from both analyses show again an inconsistency in the descriptions of these observables when using the X max as a reference. The discrepancy in this case is not as large as in the case of the pure muonic observables (N µ and X µ max ) and it might be also related to the deficit of muons in simulations.
More precise constraints on the hadronic models will be possible in the future with AugerPrime, the upgrade of the Pierre Auger Observatory [20]. By using dedicated muon detectors of two types (SSD and AMIGA), more information about the muonic component will be obtained and valuable experimental input will be provided for the study of hadronic interactions.  Figure 5. Results of the ∆ s method analysis [19].