**Testing hadronic interaction models using hybrid data of the Pierre Auger Observatory**

In several experimental results, an inconsistency is observed between inferences on the mass composition of ultrahigh-energy cosmic rays (UHECR, energies > 10^{18} eV) from the measurements of the depth of shower maximum (X_{max}) and signals in surface detectors (SD). To describe the SD data, sensitive to the muon air-shower content, heavier primary masses are needed than for fitting the X_{max} data only. This discrepancy is generally attributed to the lack of muons ('muon puzzle') in the air-shower simulations. In such an interpretation, one explicitly assumes that the uncertainties on the predicted X_{max} scale are delimited by the current hadronic interaction models tuned to the LHC data (EPOS-LHC, QGSJet II-04, Sibyll 2.3d). In the approach used in the paper, we disregard this assumption and test predictions of these models considering both the SD signals and the X_{max} scale.

We use 2239 high-quality hybrid events, i.e. the events recorded with both the fluorescence detector and the SD of the Pierre Auger Observatory, with the energies 10^{18.5} eV − 10^{19.0} eV. The test consists in fitting the observed two-dimensional distributions of X_{max} and SD signal at 1000 m from the shower core, S(1000), with simulated templates simultaneously in five zenith-angle (θ) ranges. The free parameters of the fit are the X_{max} scale shift (∆X_{max}) and the rescaling of the hadronic part of S(1000) at two extreme zenith angles, Rhad (θ_{min} ≈ 28˚) and Rhad (θ_{max} ≈ 55˚), of a corresponding hadronic interaction model, and the fractions of four primary nuclei (hydrogen, helium, oxygen and iron).

The surprising outcome is that the best fit of the data distributions is achieved for the X_{max} scales of all three hadronic models shifted by ~20-50 g/cm^{2} to deeper values, see figure 1. As a consequence of the shifts, the main differences in model predictions in both average X_{max} and S(1000) can be reduced leading to the similar inferences on the mass composition within this method. It is remarkable that for all three models, a consistent description of the X_{max} and SD data can be achieved only with the simultaneous modification of two scale parameters at once; modifications of either hadronic signal or Xmax scale alone do not lead to a good data description.

For the first time, we find a 5σ tension between the Auger data and hadronic model predictions, even when accounting for combinations of experimental systematic uncertainties, see figure 2. The shifts of the predicted X_{max} scales lead to heavier primary mass compositions and as a consequence to a smaller, 15% - 25%, deficit in the simulated hadronic signal (dominated by muons), see figure 3, compared to previous much larger estimations of ~30% - 60% that assumed the X_{max} scales predicted by the models. The difference in Rhad (θ) at the two extreme zenith angles for QGSJet II-04 might indicate that softer muon spectra at 1000m from the shower core than in this model would describe the Auger data better.

The specific ways to produce the required changes in the models, which might consist of combinations of modifications of integral (cross-section, multiplicity, elasticity, etc.) or differential (secondary particle energy spectra) characteristics of the hadronic interactions, are beyond the scope of the paper.

In conclusion, our results suggest that the inferences on the mass composition of UHECR based on the X_{max} data might not be bracketed by the current models of hadronic interactions and that the UHECR mass composition can be heavier than in currently existing X_{max}-based estimations.

**Related paper:**

**Testing Hadronic-Model Predictions of Depth of Maximum of Air-Shower Profiles and Ground-Particle Signals using Hybrid Data of the Pierre Auger Observatory**

The Pierre Auger Collaboration, Phys. Rev. D 109, 102001 (2024)

[arxiv.org/abs/2401.10740] [doi: 10.1103/PhysRevD.109.102001]