Updated Air-Shower $X_{\rm max}$ Moment Parametrizations for UHECR Composition with Latest Hadronic Interaction Models

TL;DR

This work delivers an updated, self-consistent interface between air-shower simulations and UHECR composition analyses by parametrizing $X_{\max}$ moments and full distributions across three modern hadronic interaction models. It provides four equivalent representations: (i) $\langle X_{\max}\rangle$ and $\sigma^2(X_{\max})$ as functions of energy and mass with two variance forms (polynomial and exponential); (ii) a mapping from $X_{\max}$ moments to $\ln A$ moments via intrinsic-shower fluctuations; (iii) a generalized three-parameter Gumbel model $\mathcal{G}(X;\mu,\sigma,\lambda)$ for the full $X_{\max}$ distribution with $\epsilon$ and $Y$-dependent parameters, fitted by unbinned likelihood; and (iv) validated internal consistency showing residuals at the few g cm$^{-2}$ level for the mean and a few (g cm$^{-2}$)$^2$ for the variance. The results enable precise forward-folding analyses and robust cross-model comparisons, and the authors publicly release data products to facilitate adoption in composition studies. These updates address hadronic-model uncertainties and reflect the latest HIM developments, with explicit comparisons to older model versions highlighting the impact on $X_{\max}$ observables and composition inferences.

Abstract

The mass composition of ultra-high-energy cosmic rays (UHECRs) is commonly inferred from the first two moments of the depth of shower maximum, $X_{\rm max}$, measured by fluorescence and hybrid detectors. Such analyses require fast and accurate mappings between the moments of $X_{\rm max}$ and those of the logarithmic mass, $\ln A$, based on realistic air-shower simulations. In this work we provide updated parametrizations of the $X_{\rm max}$ moments and distributions for air showers initiated by nuclei from proton to iron, simulated with CONEX for three state-of-the-art hadronic interaction models: Epos LHC-R, Sibyll 2.3e, and QGSJet-III-01. We parametrize the mean depth $\langle X_{\rm max}\rangle$ and the variance $σ^2(X_{\rm max})$ as functions of energy and mass. For the variance we compare a second-order polynomial model with an exponential model. In addition, we model the full $X_{\rm max}$ distributions with a three-parameter generalized Gumbel function. The Gumbel parameters are fitted using an unbinned likelihood and are validated by comparing the implied mean and variance with the raw CONEX samples and with the moment parametrizations. Across the full energy range considered, residuals between the parametrizations (or the Gumbel representation) and the simulations are at the level of a few g cm$^{-2}$ for the mean and a few (g cm$^{-2}$)$^2$ for the variance, making these parametrizations suitable for precision UHECR composition studies and forward-folding analyses of $X_{\rm max}$ distributions.

Figures (10)

• Figure 1: Mean $X_{\mathrm{max}}$ as a function of energy for protons (left) and iron (right), as obtained for the three hadronic interaction models: Sibyll 2.3e, EPOS LHC-R, and QGSJet III-01. The lower panels display the systematic differences of Sibyll 2.3e and QGSJet III-01 with respect to EPOS LHC-R.
• Figure 2: Standard deviation of $X_{\mathrm{max}}$ as a function of energy for protons (left) and iron (right), as obtained for the three hadronic interaction models: Sibyll 2.3e, EPOS LHC-R, and QGSJet III-01. The lower panels display the systematic differences of Sibyll 2.3e and QGSJet III-01 with respect to EPOS LHC-R.
• Figure 3: Comparison between the parametrizations and the CONEX simulations for $\langle X_{\mathrm{max}} \rangle$ for the three hadronic interaction models: EPOS LHC-R (left), Sibyll 2.3e (middle), and QGSJet III-01 (right). The upper panels show the simulated values (dots) together with the analytical parametrizations (solid lines) for representative primary nuclei. The lower panels display the residuals, defined as in Eq. \ref{['eq:bias']}.
• Figure 4: Same as Fig. \ref{['fig:xmax_residual']}, but for $\sigma(X_{\mathrm{max}})$ using the second-order polynomial parametrization of Eq. \ref{['eq:xmax_sigma2']}.
• Figure 5: Same as Fig. \ref{['fig:xmax_residual']}, but for $\sigma(X_{\mathrm{max}})$ using the exponential parametrization of Eq. \ref{['eq:xmax_sigma2_exp']}.