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 > 1018 eV) from the measurements of the depth of shower maximum (Xmax) 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 Xmax 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 Xmax 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 Xmax 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 1018.5 eV − 1019.0 eV. The test consists in fitting the observed two-dimensional distributions of Xmax 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 Xmax scale shift (∆Xmax) 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).
Figure 1: The energy evolution of the mean of the Xmax distributions measured by the Pierre Auger Observatory using fluorescence (solid circles) and surface (open squares) detectors. The results of this paper for the modified Xmax scales of models are shown with shaded bands with the heights corresponding to the systematic uncertainties. The original non-modified model predictions for different primary species are shown with lines for the entire energy range.
The surprising outcome is that the best fit of the data distributions is achieved for the Xmax scales of all three hadronic models shifted by ~20-50 g/cm2 to deeper values, see figure 1. As a consequence of the shifts, the main differences in model predictions in both average Xmax 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 Xmax 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.
Figure 2: Values of the modification parameters (points) for all possible combinations of experimental systematic uncertainties on the energy (±14%), Xmax (+8−9 g/cm2 ) and S(1000) (±5%). The color of the points shows the difference in log-likelihood expressions in the case of no modifications and in the case of the assumed template modifications, including the differences higher than 50 (note the slightly different scale between models). The value estimated using Wilks’ theorem in the Likelihood-ratio test for the nested model at the level of 5σ corresponds to ∆ ln L ≈ 16.62. The results for no systematic shift of the data are highlighted by stars. Dashed lines outline the contour of the plane from the best fit to the points. The closest approach to the non-modified (Rhad(θ) = 1, ∆Xmax = 0 g/cm2) model predictions using a dense scan of linear combinations of experimental systematic uncertainties is connected with this point by a black line. The animated rotated views are available at [https://doi.org/10.5281/zenodo.10653685].
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 Xmax 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 Xmax 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.
Figure 3: Correlations between ∆Xmax and Rhad (θmax ≈55˚) modifications of the model predictions obtained from the data fits. The contours correspond to 1σ, 3σ, and 5σ statistical uncertainties. The gray rectangles are the projections of the total systematic uncertainties.
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 Xmax 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 Xmax-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]
