# 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  The Pierre Auger Observatory, Av. San Martín
Norte 306, 5613 Malargüe, Mendoza, Argentina;
http://www.auger.org<EMAIL_ADDRESS>
###### Abstract

We test the predictions of hadronic interaction models regarding the depth of
maximum of air-shower profiles, $X_{\text{max}}$, and ground-particle signals
in water-Cherenkov detectors at 1000 m from the shower core, $S(1000)$, using
the data from the fluorescence and surface detectors of the Pierre Auger
Observatory. The test consists in fitting the measured two-dimensional
($S(1000)$, $X_{\text{max}}$) distributions using templates for simulated air
showers produced with hadronic interaction models Epos-LHC, QGSJet-II-04,
Sibyll 2.3d and leaving the scales of predicted $X_{\text{max}}$ and the
signals from hadronic component at ground as free fit parameters. The method
relies on the assumption that the mass composition remains the same at all
zenith angles, while the longitudinal shower development and attenuation of
ground signal depend on the mass composition in a correlated way.

The analysis was applied to 2239 events detected by both the fluorescence and
surface detectors of the Pierre Auger Observatory with energies between
$10^{18.5}$ to $10^{19.0}\,\text{eV}$ and zenith angles below $60^{\circ}$. We
found, that within the assumptions of the method, the best description of the
data is achieved if the predictions of the hadronic interaction models are
shifted to deeper $X_{\text{max}}$ values and larger hadronic signals at all
zenith angles. Given the magnitude of the shifts and the data sample size, the
statistical significance of the improvement of data description using the
modifications considered in the paper is larger than $5\sigma$ even for any
linear combination of experimental systematic uncertainties.

††preprint: APS/123-QED

##  I Introduction

The dominant contribution to uncertainties in the determination of the mass
composition of ultra-high-energy cosmic rays (UHECR, energy
$E>10^{18.0}\,\text{eV}$) comes from the modelling of extensive air showers.
Modern hadronic interaction models used for this purpose are based on
extrapolation of interaction parameters like cross-sections, multiplicities,
elasticities, etc. measured at accelerators at lower beam energies up to
$\sqrt{s}=13$ TeV for proton-proton collisions at the LHC and pseudorapidities
$|\eta|\lesssim 5$ compared to energies $\sqrt{s}\gtrsim 50$ TeV and
pseudorapidities111The forward calorimeters at LHC can cover $\eta=8.4$ to 15,
but of neutral particles only. $\eta\approx 7$ to 11 driving the energy flow
of the first interactions of UHECR in the atmosphere, where the target is
different (mostly oxygen and nitrogen nuclei). Therefore, improvements in the
description of the LHC data, implemented in the modern models do not
necessarily lead to unambiguous, nearly hadronic-model independent,
predictions for the mass-sensitive air-shower observables. For instance, at
$10^{18.7}\,\text{eV}$ the span in predictions for the mean depth of maximum
of air-shower profiles, $\langle{X_{\text{max}}}\rangle$, between models used
by the UHECR community (Epos-LHC [1], QGSJet-II-04 [2], Sibyll 2.3d [3]) is
${\sim}25$ g/cm2, nearly independently of the primary particle mass and
energy. Such a difference can be considered only a lower limit on the
systematic uncertainty of the predicted $X_{\text{max}}$ scale. This is about
one-quarter of the difference between $\langle X_{\text{max}}\rangle$ values
of the two astrophysical extremes, protons and iron nuclei. As a consequence,
the mass composition of cosmic rays can be referred only with respect to
$X_{\text{max}}$ scale predicted by a particular model. The largest
differences in the predicted $\langle{X_{\text{max}}}\rangle$ between the
models, with a minimal impact on the elongation rate, come from the properties
of the first hadronic interaction and production of nucleons-antinucleons in
pion-air and kaon-air interactions that are not known well in the relevant
kinematic region (for more detailed discussion, see e.g. Ref. [4]). The
general comparison of properties and different treatments of hadronic
interactions by the three models used in this work is summarized in Ref. [5].
Note that although the models used in this work treat differently the
properties of hadronic interactions, the range of air-shower properties
predicted by these models does not need to include all the possibilities.

The difference in the standard deviation of $X_{\text{max}}$ distributions,
$\sigma(X_{\text{max}})$, between the models is within ${\sim}5$ g/cm2,
whereas the difference between the $X_{\text{max}}$ fluctuations of protons
and iron nuclei is ${\sim}40$ g/cm2. Therefore, the difference of
$\sigma(X_{\text{max}})$ in model predictions has a smaller effect on the mass
composition inferences compared to the difference in predictions of
$\langle{X_{\text{max}}}\rangle$. Note, that there is no direct correspondence
of the model scales of $X_{\text{max}}$ and $\sigma(X_{\text{max}})$, i.e.
differences in $\sigma(X_{\text{max}})$ are not a mere consequence of the
differences in $\langle{X_{\text{max}}}\rangle$.

In general, the signal produced by air-shower particles reaching the ground
shows much lower sensitivity to the mass composition than in the case of
$X_{\text{max}}$. The model differences in predictions of the ground-particle
signal, for instance, at 1000 m from the impact point of an air shower at
$10^{18.7}\,\text{eV}$ detected at the Pierre Auger Observatory (Auger) [6]
are at the level of ${\sim}3$ VEM222This unit is the signal produced by a muon
traversing the station on a vertical trajectory., whereas the difference
between protons and iron nuclei is about twice this value. The fluctuations of
the ground signal are dominated by the detector resolution, suppressing any
significant sensitivity to the model or primary mass [7].

Table 1: Indicative summary of the results of tests of models using Auger data. In the case of Sibyll 2.3d, we also show estimations based on the previous version of the model Sibyll when available in the literature. The red shade of the boxes reflects the size of the tension between the models and measurements. test | energy/EeV | $\theta/^{\circ}$ |  Epos-LHC |  QGSJet-II-04 |  Sibyll 2.3d
---|---|---|---|---|---
$X_{\text{max}}$ moments [8, 9, 10, 11] |  ${\sim}3$ to 50 | 0 to 80 | no tension $\blacksquare$ | tension $\blacksquare$ | no tension (2.3c) $\blacksquare$
$X_{\text{max}}$:$S(1000)$ correlation [12, 11] | 3 to 10 | 0 to 60 | no tension $\blacksquare$ | tension $\blacksquare$ | no tension (2.3c) $\blacksquare$
mean muon number [13, 14] | ${\sim}10$ | ${\sim}67$ | tension $\blacksquare$ | tension $\blacksquare$ | tension $\blacksquare$
mean muon number [15] | 0.2 to 2 | 0 to 45 | tension $\blacksquare$ | tension $\blacksquare$ | —
fluctuation of muon number [14] | 4 to 40 | ${\sim}67$ | no tension $\blacksquare$ | no tension $\blacksquare$ | no tension $\blacksquare$
muon production depth [16] | 20 to 70 | ${\sim}60$ | tension $\blacksquare$ | no tension $\blacksquare$ | —
$S(1000)$ [17] | ${\sim}10$ | 0 to 60 | tension $\blacksquare$ | tension $\blacksquare$ | —

###  I.1 Problems in the Description of Data with Hadronic Interaction Models

The correctness of the model predictions can be tested using data from air-
shower experiments. At the Pierre Auger Observatory, for instance, negative
variance of the logarithm of primary masses ($\sigma^{2}(\ln A)$) and poor
description of the measured $X_{\text{max}}$ distributions are obtained when
simulations with QGSJet-II-04 are used for the interpretation of the measured
$X_{\text{max}}$ moments in terms of the moments of $\ln A$ for $E\gtrsim
10^{18.5}\,\text{eV}$ [8, 9, 10, 11]. The problem is that, due to relatively
shallow $\langle X_{\text{max}}\rangle$ predictions of QGSJet-II-04, the best
possible description of the data is achieved with proton-helium mixes, but for
these mixes the modeled $X_{\text{max}}$ distributions are broader than the
observed ones. These findings are further supported by the observation of the
negative correlation between $X_{\text{max}}$ and the signal in surface
detector (SD) [12, 11] in the Auger data near the “ankle” feature ($E\approx
10^{18.7}\,\text{eV}$) of the UHECR spectrum. That can be achieved only for
the mixed primary composition containing particles heavier than helium nuclei,
which is a robust statement based on the general phenomenology of air-shower
development and, therefore, not sensitive to the properties of a particular
hadronic model.

The deficiencies of the models are more evident for the SD observables, where,
in some cases, the data are not even bracketed by the Monte Carlo (MC)
predictions for protons and iron nuclei. The muon deficit in simulations,
known as the “muon puzzle” [5], is the best-known example of this kind. In
particular, this was observed by Auger for inclined showers (zenith angles
$\theta=62^{\circ}$ to $80^{\circ}$), dominated by the muon component, for
Epos-LHC and QGSJet-II-04 at ${\sim}10^{19.0}\,\text{eV}$ [13], and in direct
measurements with underground muon detectors at $E=10^{17.3}\text{ to
}10^{18.3}\,\text{eV}$ and $\theta<45^{\circ}$ [15]. At the same time, the
fluctuations of the muon signal measured by Auger [14] are consistent with the
MC predictions, indicating that the muon deficit might originate from the
accumulation of small deviations from the model predictions during the
development of a shower, rather than be caused by a strong deviation in the
first interaction. The range of predictions for the muon production depth of
protons and iron nuclei is outside of the measured values for Epos-LHC above
$10^{19.3}\,\text{eV}$ [16].

Less directly, the muon deficit in simulations was observed as a deficit of
the total signal at 1000 m from the shower core, $S(1000)$, in the Auger SD
stations [17] for vertical ($\theta<60^{\circ}$) showers with energies around
$10^{19.0}\,\text{eV}$. In this analysis, within the assumption that the
electromagnetic (em) component, $X_{\text{max}}$, and, correspondingly, the
mass composition inferred from the $X_{\text{max}}$ measurements, are
predicted correctly by a particular model, the deficit of $S(1000)$ was
interpreted as the underestimation of the hadronic signal (dominated by muons)
in simulations by $(33\pm 16)\%$ for Epos-LHC and by $(61\pm 21)\%$ for
QGSJet-II-04, with a strong dependence on the energy scale. An indicative
summary of all these tests for the three models used in this work is given in
Table 1.

The use of an incorrect MC $X_{\text{max}}$ scale would lead one to a biased
inference on the mass composition and, through this, to a biased estimate of
the muon deficit, since the muon content in a shower scales with the primary
mass $\propto A^{1-\beta}$, $\beta\approx 0.9$ [18]. Therefore, in a more
comprehensive approach a modification of the MC scales of both
$X_{\text{max}}$ and SD signals, going along with the fitting of the primary
mass composition accounting for these modifications, should be considered.

###  I.2 Progressive Testing of Hadronic Interaction Models

In this work, progressive testing of the model predictions is performed.
First, we allow for a rescaling of the signal on the ground produced by the
hadronic shower component at 1000 m with a factor $R_{\text{had}}$. Then we
add a zenith-angle dependence of $R_{\text{had}}$ defining two parameters
$R_{\text{had}}(\theta_{\text{min}})$ and
$R_{\text{had}}(\theta_{\text{max}})$ at the two extreme zenith-angle bins. We
assume a linear dependence of $R_{\text{had}}$ on the distance of
$X_{\text{max}}$ to the ground in atmospheric depth units to relate
$R_{\text{had}}$ to different zenith angles. Finally, we consider a shift in
the predicted $X_{\text{max}}$ distributions ($\Delta X_{\text{max}}$) that is
assumed to be independent of the primary mass and energy. In this way, we
consider freedom not only in the scale of the simulated hadronic part of the
ground signal but also in the simulated $X_{\text{max}}$ scale. Consequently,
the main differences in $X_{\text{max}}$ and $S(1000)$ predictions of the
models are reduced, and similar mass composition inferences for the Auger data
are obtained. It is remarkable, that a consistent description of the Auger
data with all 3 models can be achieved with the modification of only two scale
parameters.

The analysis is performed for the energy region around the ankle ($E\approx
10^{18.7}\,\text{eV}$) in the energy spectrum where the UHECR mass composition
is mixed [9, 12, 11]. Specifically, we find the values of $\Delta
X_{\text{max}}$, $R_{\text{had}}(\theta)$ and fractions of four primary
particles (protons, and helium, oxygen, iron nuclei) for which the best fit of
the measured two-dimensional distributions of ($S(1000)$, $X_{\text{max}}$) is
achieved. The remaining differences between the predictions of models with a
smaller effect on MC templates like the fluctuations of $X_{\text{max}}$ and
hadronic signal, and the mass composition dependence of
$R_{\text{had}}(\theta)$ and $\Delta X_{\text{max}}$ are not considered in
this work.

In Section II, we give a detailed description of the method. In Section III,
the method is applied to the data of the Pierre Auger Observatory. The results
of the data analysis are discussed in Section IV followed by a summary of our
findings in Conclusions.

##  II Method

The method stems from Ref. [19] and it was first introduced in Ref. [20]
followed by a slight modification in the approach to the ground-signal
rescaling. A preliminary application of the method to the Auger data was
presented in Ref. [21].

We perform a binned maximum-likelihood fit of the measured two-dimensional
distributions ($S$, $X$) with MC templates for four primary species
simultaneously in five zenith-angle bins. The observables are $S(1000)$ and
$X_{\text{max}}$ corrected for the energy evolution using the fluorescence
detector (FD) energy $E_{\text{FD}}$,

$S=S(1000)\,\left(\frac{E^{\text{ref}}}{E_{\text{FD}}}\right)^{1/B},$ (1)

and

$X=X_{\text{max}}+D\,\lg\left(\frac{E^{\text{ref}}}{E_{\text{FD}}}\right),$
(2)

where $B=1.031\pm 0.004$ is the SD energy calibration parameter [22] and the
elongation rate of a single primary $D=58$ g/cm2/decade is taken as the
average value over the four primary particles and three models used in this
work. The value of the elongation rate varies only within $\pm 2$ g/cm2 around
the mean, in accordance with the universality with respect to the primary mass
predicted within the phenomenological model in Ref. [23]. This is a
consequence of energy dependencies of multiplicity, elasticity, and cross-
section assumed in this phenomenological model. We chose the reference energy
$E^{\text{ref}}=10^{18.7}\,\text{eV}$ for the analyzed energy range of
$E_{\text{FD}}=10^{18.5}\text{ to }10^{19.0}\,\text{eV}$.

The signal $S(1000)$ is assumed to be composed of the hadronic
($S_{\text{had}}$) and electromagnetic ($S_{\text{em}}$) components. The
signal $S_{\text{had}}$ is produced by muons, em particles from muon decays
and low-energy neutral pions as in [17] according to the four-component shower
universality model [24, 25]. The signal $S_{\text{em}}$ is produced by em
particles originating from high-energy neutral pions.

###  II.1 Monte Carlo Templates

The MC templates were prepared using simulated showers from a library produced
within the Auger collaboration [26]. Air showers were generated with Corsika
7.7400 [27, 28, 29] with a flat zenith-angle distribution in $\sin^{2}\theta$
(for $\theta\in(0^{\circ},65^{\circ})$) and the energy distribution $\propto
E^{-1}$. Subsequently, we re-weight the events to match the measured energy
spectrum [22]. Four different atmospheric profiles are considered to represent
the typical variations of atmospheric conditions at the location of the Pierre
Auger Observatory. The simulation set includes three models (Epos-LHC, QGSJet-
II-04, Sibyll 2.3d) and four primary particles (p, He, O, Fe). In productions
with Epos-LHC and Sibyll 2.3d, the low-energy ($E_{\text{kin}}<80$ GeV)
interactions were simulated with Urqmd [30], while QGSJet-II-04 was used in
combination with Fluka [31, 32]. No significant dependency of the results on
the choice of the low-energy model was found.

The detector simulations and event reconstruction were performed with the
Auger $\overline{\textrm{Off}}$ $\underline{\textrm{line}}$ Framework [33]. In
the standard event processing chain, not all effects related to the detector
calibration, atmospheric conditions, long-term performances, etc. are included
in the simulations and reconstruction. To account for them, for the FD part an
additional smearing of $X_{\text{max}}$ distributions by $\approx$9 g/cm2 is
applied [9]. In the case of the SD part, the smearing of $S(1000)$ by 9%,
corresponding to the maximal expected contribution from realistic operational
conditions, was tested without finding any statistically significant effect on
the results. The event selection is the same as applied to the data (see
Section III). After the selection, the MC templates contain $\approx 15000$
showers per primary specie and model.

The analysis was carried out by splitting these simulated air-showers in five
zenith-angle ranges containing nearly the same number of events, namely
$(0^{\circ},33^{\circ})$, $(33^{\circ},39^{\circ})$,
$(39^{\circ},45^{\circ})$, $(45^{\circ},51^{\circ})$,
$(51^{\circ},60^{\circ})$. Examples of the ($S$, $X$) distributions for
protons and iron nuclei in the most vertical and most inclined angular ranges
are shown in Fig. 1. Such two-dimensional distributions are then normalized
and fitted with the ansatz function $\Phi$ described in detail in Appendix A.
This function is a convolution of the generalized Gumbel distribution of $X$
and the Gaussian distribution of $S$ with the mean value linearly changing
with $X$, reflecting in this way their correlation indicated in Fig. 1. A set
of these trial functions for each model, primary particle, and zenith-angle
range is used as MC templates in the following fitting procedure.

Figure 1: Examples of two-dimensional distributions of $S$ and $X$ for protons
(left) and iron nuclei (right) generated with Epos-LHC for zenith angles
between $0^{\circ}$ and $33^{\circ}$ (top) and $51^{\circ}$ and $60^{\circ}$
(bottom). The red points indicate the mean values of $S$. $E=10^{18.5}\text{
to }10^{19.0}\,\text{eV}$.

Figure 2: The total ground signal at 1000 m from the shower core (black) and
its hadronic (brown) and em (turquoise) components as a function of the
distance from $X_{\text{max}}$ to the ground in atmospheric depth units for
protons (left) and iron nuclei (right) for different models. The bands contain
the statistical uncertainty. $E=10^{18.5}\text{ to }10^{19.0}\,\text{eV}$,
$\theta<60^{\circ}$.

###  II.2 Fitting Procedure

For each model, we search for the most likely combination of the composition
mix of the four primary species, the zenith-dependent rescaling parameter of
the hadronic signal $R_{\text{had}}(\theta)$, and the constant shift of the
depth of shower maximum $\Delta X_{\text{max}}$ in all the MC templates. The
fitting method is a generalization of the fitting procedure used in Ref. [10]
in the case of the $X_{\text{max}}$ distribution and applied here to the two-
dimensional ($X$, $S$) distributions in five zenith-angle ranges
simultaneously.

The negative log-likelihood-ratio expression that is minimized for a given
model is of the form

$\ln\mathcal{L}=\begin{cases}\sum\limits_{k}\sum\limits_{j}(C_{jk}-n_{jk}+n_{jk}\,\ln\frac{n_{jk}}{C_{jk}}),&n_{jk}>0,\\\
\sum\limits_{k}\sum\limits_{j}C_{jk},&n_{jk}=0,\end{cases}$ (3)

with the sums running over the two-dimensional bins $j$ in ($X$, $S$) for the
five $\theta$-bins $k$. The corresponding number of showers measured in bins
$j$, $k$ is denoted by $n_{jk}$ and the predicted number of MC showers by
$C_{jk}$. The latter number is obtained using the total number of measured
showers $N_{\text{data}}^{k}$ as

$C_{jk}=N_{\text{data}}^{k}\,\sum_{i}f_{i}\,\Phi_{i,k}(X^{\prime}_{\text{j,k}},S^{\prime}_{\text{j,k}}),$
(4)

where $\Phi_{i,k}$ denotes the template function $\Phi$ in $\theta$-bin $k$
for a given model and primary particle $i$ with relative fraction $f_{i}$. The
modified $X$ prediction is of the form

$X^{\prime}_{\text{j,k}}=X_{\text{j,k}}+\Delta X_{\text{max}},$ (5)

and the rescaled predicted ground signal

$S^{\prime}_{\text{j,k}}=S_{\text{j,k}}\,f_{\text{SD}}^{k}(R_{\text{had}}(\theta),\Delta
X_{\text{max}}),$ (6)

where $X_{\text{j,k}}$ and $S_{\text{j,k}}$ are the center bin values of $X$
and $S$, respectively, of the original MC distribution ($X$, $S$).

Figure 3: The average fraction of hadronic signal at 1000 m from the shower
core as a function of the reconstructed zenith angle for different models and
primary masses. $E=10^{18.5}\text{ to }10^{19.0}\,\text{eV}$.

The rescaling parameter $f_{\text{SD}}^{k}$ of all signals $S_{\text{j,k}}$ is
calculated as

$\displaystyle f_{\text{SD}}^{k}=R_{\text{had}}(\theta)\,$
$\displaystyle\frac{(E^{\text{ref}})^{\beta-1/B}}{\langle{E_{\text{FD}}^{\beta-1/B}}\rangle_{k}}\,g_{\text{had},k}\,f_{\text{had},k}~{}+$
(7) $\displaystyle+~{}$
$\displaystyle\frac{(E^{\text{ref}})^{1-1/B}}{\langle{E_{\text{FD}}^{1-1/B}}\rangle_{k}}\,g_{\text{em},k}\,(1-f_{\text{had},k}),$

where $S(1000)\propto E_{\text{FD}}^{B}$ is assumed to be composed only of
$S_{\text{had}}\propto E_{\text{FD}}^{\beta}$ and $S_{\text{em}}\propto
E_{\text{FD}}$, and $f_{\text{had}}=S_{\text{had}}/S(1000)$. The parameter
$\beta=0.92$ is chosen following Ref. [23]. The mean energy factors
$\langle{E_{\text{FD}}^{\beta-1/B}}\rangle_{k}$ and
$\langle{E_{\text{FD}}^{1-1/B}}\rangle_{k}$ are calculated from all measured
showers in the energy range $10^{18.5}\text{ to }10^{19.0}\,\text{eV}$ and
$\theta$-bin $k$.333Note that the choice of
$E^{\text{ref}}\approx\langle{E_{\text{FD}}}\rangle$ minimizes the effect of
the energy factors:
$(E^{\text{ref}})^{\beta-1/B}/\langle{E_{\text{FD}}^{\beta-1/B}}\rangle_{k}\approx
1$ and
$(E^{\text{ref}})^{1-1/B}/\langle{E_{\text{FD}}^{1-1/B}}\rangle_{k}\approx 1$.
The mean hadronic fraction, see also Fig. 3, in the $\theta$-bin $k$,
$f_{\text{had},k}$, is calculated using the average hadronic fractions
$f_{\text{had},k,i}$ for simulated showers induced by a primary $i$ and
weighted over the relative primary fractions $f_{i}$ as

$f_{\text{had},k}=\sum_{i}f_{i}\,f_{\text{had},k,i}.$ (8)

The average effect of the $X_{\text{max}}$ change on the ground signal is
incorporated through the separate effects on the em $g_{\text{em},k}$ and
hadronic $g_{\text{had},k}$ signals, see Appendix C for details. We
parameterized the evolution of the mean ground signal parts with the distance
of $X_{\text{max}}$ to the ground in atmospheric depth units,
$X_{\text{atm}}-X_{\text{max}}$, where
$X_{\text{atm}}=880\,\text{g/cm}^{2}/\cos\theta$, see the examples in Fig. 2.
In this way, the total ground signal is estimated to be modified (via
$g_{\text{em},k}$, $g_{\text{had},k}$) at most by about 7% for a change of
$X_{\text{max}}$ by 50 g/cm2. This is in accordance with the functional
dependencies in Fig. 2 weighted over the relative contributions of hadronic
and em signals (see Fig. 3) and over the primary fractions.

Figure 4: Distributions of $S$ and $X$ for the data of the Pierre Auger
Observatory in the energy range $10^{18.5}\text{ to }10^{19.0}\,\text{eV}$ in
five zenith-angle bins.

As a consequence of the four-component shower universality approach, the em
signal is very similar in all three models, see Fig. 2. The differences in the
total signal stem from the size of the hadronic signal at different zenith
angles, corresponding to different $X_{\text{atm}}-X_{\text{max}}$ values.
Therefore, the freedom in $R_{\text{had}}(\theta)$ and also $\Delta
X_{\text{max}}$ (see Section I) removes the main differences in predictions of
$S$ and $X$ for the three models. We assumed a linear dependence of
$R_{\text{had}}(\theta)$ on $X_{\text{atm}}-X_{\text{max}}$ and defined
rescaling parameters of the hadronic signal at two extreme zenith angles,
$R_{\text{had}}(\theta_{\text{min}})$ and
$R_{\text{had}}(\theta_{\text{max}})$, for ${\sim}28^{\circ}$ and
${\sim}55^{\circ}$, respectively, see Appendix C for the definition.

We have verified using MC-MC tests, see Appendix B, that the method is
performing well and the observed biases were taken as sources of systematic
uncertainty of the results.

Table 2: Evolution of the minimum values of the log-likelihood expression (see Eq. 3) fitting the data with different modifications of the model predictions. In all cases, except $R_{\text{had}}=\text{const.}$ and $\Delta X_{\text{max}}$, the significance of improvement of data description with the $R_{\text{had}}(\theta)$ and $\Delta X_{\text{max}}$ fit is above $5\sigma$ using the Likelihood-ratio test applying the Wilks’ theorem [34] for nested model [35]. In the case of $R_{\text{had}}=\text{const.}$ and $\Delta X_{\text{max}}$, the improvement of data description with $R_{\text{had}}(\theta)$ and $\Delta X_{\text{max}}$ fit is ${\sim}0.1\sigma$, ${\sim}4.4\sigma$ and ${\sim}2.0\sigma$ for Epos-LHC, QGSJet-II-04 and Sibyll 2.3d, respectively. $\ln\mathcal{L}_{\text{min}}$ | Epos-LHC |  QGSJet-II-04 |  Sibyll 2.3d
---|---|---|---
none | 2022.9 | 4508.0 | 2496.5
$\Delta X_{\text{max}}$ | 738.6 | 1674.8 | 1015.7
$R_{\text{had}}=\text{const.}$ | 489.2 | 684.4 | 521.6
$R_{\text{had}}(\theta)$ | 489.2 | 673.9 | 517.6
$R_{\text{had}}=\text{const.}$ and $\Delta X_{\text{max}}$ | 452.2 | 486.7 | 454.2
$R_{\text{had}}(\theta)$ and $\Delta X_{\text{max}}$ | 451.9 | 476.3 | 451.6

##  III Data Analysis

###  III.1 Data Selection

The analysis is applied to hybrid events, i.e. events detected with both SD
and FD of the Pierre Auger Observatory between 1 January 2004 and 31 December
2018. To select high-quality events, for the FD part we apply the selection
criteria used for the $X_{\text{max}}$ analysis [9, 11], for the SD part the
selection is the same as in the measurements of the SD energy spectrum [22],
additionally removing events with saturated SD stations. In this way, we
ensure an accurate estimation of the observables $X_{\text{max}}$,
$E_{\text{FD}}$, $\theta_{\text{FD}}$, and $S(1000)$ used as the inputs to the
method. The selection efficiency is similar for all four primary masses and
all zenith-angle bins, therefore it does not introduce mass-composition
biases. In total, 2239 hybrid events were selected in the FD energy range
$10^{18.5}\text{ to }10^{19.0}\,\text{eV}$
($\langle{E_{\text{FD}}}\rangle\approx 10^{18.7}\,\text{eV}$) and zenith
angles between $0^{\circ}$ and $60^{\circ}$. The data were divided into five
zenith-angle bins containing nearly the same ($N=425-500$) and sufficiently
large number of events, see Fig. 4.

To take into account long-term performance of the FD and SD [36], we applied a
time-dependent correction to $E_{\text{FD}}$ with a negligible systematic
effect on the final results. The signal $S(1000)$ was corrected for seasonal
atmospheric effects [37].

###  III.2 Results

Figure 5: Distributions of $X$ (top-left) and $S$ (bottom) for the Auger data
(points) and for the best fits of ($S$, $X$) distributions with
$R_{\text{had}}(\theta)$ and $\Delta X_{\text{max}}$ modifications in
predictions of three models (denoted by † as templates were modified from
original predictions). In the top-right panel, the dependence of the Gideon-
Hollister coefficient of correlation [38] between $S$ and $X$ on the zenith
angle is shown. The $\chi^{2}$ probabilities characterize the compatibility
between measurements and MC predictions in the individual plots.

The minimized values of the log-likelihood expression (see Eq. 3) are
summarized in Table 2, progressively applying different modifications to the
models. At all stages, the fits of the mass composition are performed with (p,
He, O, Fe) fractions as the free-fit parameters. First, the model predictions
without any modifications are used for fitting the data. Then we perform fits
adding freedom in only one of the modifications $\Delta X_{\text{max}}$,
$R_{\text{had}}$ (independent on zenith angle), $R_{\text{had}}(\theta)$
(zenith-angle dependent). Finally, the fits are performed using combinations
of ($R_{\text{had}}$, $\Delta X_{\text{max}}$) and ($R_{\text{had}}(\theta)$,
$\Delta X_{\text{max}}$). From the minimum values of the log-likelihood
expression for these scenarios, and from the Likelihood-ratio test for the
nested model, one can see that the major improvements in the data description
are achieved due $R_{\text{had}}$ and then $\Delta X_{\text{max}}$
modifications. The improvement in the data description due to the introduction
of the zenith-angle dependence in $R_{\text{had}}$ is statistically
significant for QGSJet-II-04, and less significant for Sibyll 2.3d , with a
negligible effect in the case of Epos-LHC.

Figure 6: _Left:_ Correlations between $\Delta X_{\text{max}}$ and
$R_{\text{had}}(\theta_{\text{max}}{\approx}55^{\circ})$ modifications of the
model predictions obtained from the data fits. The contours correspond to
1$\sigma$, 3$\sigma$, and 5$\sigma$ statistical uncertainties. The gray
rectangles are the projections of the total systematic uncertainties. _Right:_
The most likely primary fractions of the four components from the data fits
using $\Delta X_{\text{max}}$ and $R_{\text{had}}(\theta)$. The height of the
gray bands shows the size of projected total systematic uncertainties.

Figure 7: _Left:_ Correlations between $\Delta X_{\text{max}}$ and $R_{\text{had}}(\theta_{\text{min}}{\approx}28^{\circ})$ modifications of the model predictions obtained from the data fits. _Right:_ Correlation between $R_{\text{had}}(\theta_{\text{max}}{\approx}55^{\circ})$ and $R_{\text{had}}(\theta_{\text{min}}{\approx}28^{\circ})$. The contours correspond to 1$\sigma$, 3$\sigma$, and 5$\sigma$ statistical uncertainties. The gray rectangles are the projections of the total systematic uncertainties. Table 3: Modifications of the model predictions and primary fractions in the energy range $10^{18.5}\text{ to }10^{19.0}\,\text{eV}$ with statistical and systematic uncertainties for the best data fits and the $p$-values obtained using MC-MC tests. | $R_{\text{had}}(\theta_{\text{min}})$ | $R_{\text{had}}(\theta_{\text{max}})$ |  $\Delta X_{\text{max}}$/(g/cm2) |  $f_{\text{p}}$ (%) |  $f_{\text{He}}$ (%) |  $f_{\text{O}}$ (%) |  $f_{\text{Fe}}$ (%) |  $p$-value (%)
---|---|---|---|---|---|---|---|---
Epos-LHC | $1.15\pm 0.01~{}^{+0.20}_{-0.16}$ | $1.16\pm 0.01~{}^{+0.14}_{-0.10}$ | $22\pm 3~{}^{+11}_{-14}$ | $21\pm 3~{}^{+14}_{-11}$ | $20\pm 4~{}^{+15}_{-6}$ | $44\pm 5~{}^{+15}_{-6}$ | $15\pm 4~{}^{+7}_{-15}$ | 10.6
QGSJet-II-04 | $1.24\pm 0.01~{}^{+0.22}_{-0.19}$ | $1.18\pm 0.01~{}^{+0.15}_{-0.12}$ | $47^{+2}_{-1}~{}^{+9}_{-11}$ | $16\pm 2~{}^{+8}_{-10}$ | $11\pm 4~{}^{+20}_{-7}$ | $36\pm 5~{}^{+21}_{-5}$ | $37\pm 5~{}^{+6}_{-22}$ | 19.8
Sibyll 2.3d | $1.18\pm 0.01~{}^{+0.21}_{-0.17}$ | $1.15\pm 0.01~{}^{+0.15}_{-0.11}$ | $29\pm 2~{}^{+10}_{-13}$ | $13\pm 2~{}^{+18}_{-5}$ | $15\pm 4~{}^{+15}_{-12}$ | $40\pm 5~{}^{+22}_{-5}$ | $32\pm 5~{}^{+3}_{-25}$ | 32.6

The measured two-dimensional ($S$, $X$) distributions are described acceptably
by MC templates modified by $R_{\text{had}}(\theta)$ and $\Delta
X_{\text{max}}$ in all five zenith-angle bins with $p$-values estimated to be
higher than 10% for all three models. These probabilities are obtained from
MC-MC tests (see Appendix B) using distributions of values of log-likelihood
expression from the fitting of 500 random samples consisting of 2239 simulated
showers. Each MC sample had a mass composition and artificially modified
$X_{\text{max}}$ and $S_{\text{had}}(\theta)$ following the values obtained
from the best fits to data. We plot in Fig. 5 projected distributions of $X$
and $S$ (together with plots in the logarithmic scale to stress the
consistency also for the tails of distributions) and their zenith-angle
dependent correlation using the Gideon-Hollister correlation coefficient [38]
for the data and MC templates with the modifications of
($R_{\text{had}}(\theta)$, $\Delta X_{\text{max}}$). In this way, we
demonstrate that with the modified MC templates, we achieve a consistent
simultaneous description of $X_{\text{max}}$, $S(1000)$ and their correlation
for the data shown in Fig. 4. The models without modifications or only with
the zenith-angle independent modification of $R_{\text{had}}$ do not describe
the data equivalently well (see Appendix D).

The resulting parameters of the data fits with ($R_{\text{had}}(\theta)$,
$\Delta X_{\text{max}}$) modifications are presented in Table 3. For all three
models, a deeper $X_{\text{max}}$ scale is favored (see also left panels in
Figs. 6 and 7) that would result in heavier primary mass composition derived
from the $X_{\text{max}}$ data compared to the inferences with non-modified
predictions of the models [10, 39]. The modifications $\Delta X_{\text{max}}$
are such that they reduce the difference in $X_{\text{max}}$ scales between
the models (see Fig. 10 in Section IV) and, as a consequence, similar
estimations of the fractions of the primary nuclei can be inferred from the
data, see the right panel of Fig. 6.

(a)  (b)

Figure 8: _Left:_ The total rescaling of $S$ (dashed lines) broken down into
the contributions from the rescaling of the hadronic signal
$R_{\text{had}}(\theta)$ (solid lines) and the change of the predicted
$X_{\text{max}}$ scale ($\Delta X_{\text{max}}$) (dash-dotted lines) for
different models. _Right:_ The dependence of the hadronic signal at 1000 m on
the distance of $X_{\text{max}}$ to the ground in atmospheric depth units as
predicted for proton showers generated using non-modified models (dashed
lines) and accounting for the $R_{\text{had}}(\theta)$, $\Delta
X_{\text{max}}$ modifications (solid lines). The height of the gray bands
shows the size of projected total systematic uncertainties.

Due to the shift of the $X_{\text{max}}$ predictions deeper in the atmosphere,
more shower particles reach the ground producing a few percent larger SD
signals compared to the non-modified models. Therefore, the total increase of
$S$ for the modified MC predictions consists of contributions from both
$\Delta X_{\text{max}}$ (${\sim}0$ to 7%) and $R_{\text{had}}(\theta)$
(${\sim}10$ to 16%) as shown in the left panel of Fig. 8. The two effects from
the modification of the $X_{\text{max}}$ scale — heavier primary mass
composition and a larger number of shower particles reaching the ground — lead
to the values of $R_{\text{had}}(\theta)$ (Table 3) that are smaller than the
values previously found for Epos-LHC and QGSJet-II-04 in Ref. [17]. The
modified and original $S_{\text{had}}$ scales are shown for the three models
in the right panel of Fig. 8.

###  III.3 Systematic Uncertainties

There are four dominant sources of systematic uncertainties in the fitted
parameters:

  * •

the uncertainty in the FD energy scale $\pm 14\%$ [22],

  * •

the uncertainty in the $X_{\text{max}}$ measurements ${}^{+8}_{-9}$ g/cm2 [9],

  * •

the uncertainty in the $S(1000)$ measurements $\pm 5\%$ [7],

  * •

the biases of the method estimated from the MC-MC tests (see Appendix C for
the results of these tests).

Since the $X_{\text{max}}$ systematic uncertainty is strongly correlated with
the modification $\Delta X_{\text{max}}$, its effect on the nuclei fractions
is nearly cancelled out by the corresponding change of $\Delta
X_{\text{max}}$. In general, the nuclei fractions, and therefore the
inferences on the mass composition, are weakly sensitive to all experimental
systematic errors due to the simultaneous fitting of $\Delta X_{\text{max}}$
and $R_{\text{had}}$ in the method. To explore the effect of those
systematics, and as a simplifying ansatz, the total systematic uncertainties
on the fit parameters are obtained by summing all four contributions in
quadrature (see Appendix E for the size of individual contributions).

Figure 9: Values of the modification parameters (points) for all possible
combinations of experimental systematic uncertainties on the energy ($\pm
14\%$), $X_{\text{max}}$ (${}^{+8}_{-9}$ g/cm2) and $S(1000)$ ($\pm 5\%$). The
color of the points shows the difference in log-likelihood expressions
($\Delta\ln\mathcal{L}$) 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 results (see Table 3)
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 ($R_{\text{had}}(\theta)=1$, $\Delta
X_{\text{max}}=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].

One can also note, that within systematic errors no significant dependence of
$R_{\text{had}}$ on the zenith angle was found. The difference between
$R_{\text{had}}(\theta_{\text{min}})$ and
$R_{\text{had}}(\theta_{\text{max}})$ shows a tight correlation with the
uncertainty on the energy scale. However, given all the experimental
uncertainties in the case of QGSJet-II-04, the measured data prefers within
the method rather flatter attenuation of the hadronic signal at 1000 m than
predicted by the model, indicating too hard spectra of muons predicted by this
model.

The systematic uncertainties on the parameters $B$, $D$, $\beta$ used in Eqs.
1 and 2 for the energy correction of $S$, $X$, as well as corrections of the
long-term performance and other effects related to the operation of the SD and
FD, have a negligible contribution to the systematic uncertainties. We could
not identify any significant dependencies of the results on the zenith angle
or energy in the studied ranges.

###  III.4 Significance of Improvement in Data Description

In Fig. 9, the results of our method for $\Delta X_{\text{max}}$ and
$R_{\text{had}}(\theta)$ applying also all possible combinations of the
systematic uncertainties on $E_{\text{FD}}$, $X_{\text{max}}$ and $S(1000)$
are shown with the full points. These points are located approximately in a
plane, contour outlined with a dashed line, due to a correlation between the
modification parameters through the mass composition describing the data (see
the left panels of Figs. 6 and 7, e.g. increase of $\Delta X_{\text{max}}$
leads to a heavier fitted mass composition and consequently to a decrease of
$R_{\text{had}}$). The plane is tilted with respect to the
[$R_{\text{had}}(\theta_{\text{min}})$, $R_{\text{had}}(\theta_{\text{max}})$]
plane. It is a consequence of the effect of $\Delta X_{\text{max}}$ on the
ground signal $S$ at different zenith angles, see Eqs. in Appendix C, and
consequently on the fitted $R_{\text{had}}(\theta_{\text{min}})$,
$R_{\text{had}}(\theta_{\text{max}})$. The color of the points corresponds to
the difference in fitted log-likelihood expressions ($\Delta\ln\mathcal{L}$)
in case of no modifications and in case of the assumed template modifications.
The closest approach to the point of no modification is estimated through a
dense scan of linear combinations of experimental systematic uncertainties for
the lowest values of $\Delta\ln\mathcal{L}$ that are, in some cases, even
beyond the range of systematic uncertainties quoted by Auger, see Appendix F
for more details. For all three models, the closest approach of
$\Delta\ln\mathcal{L}$ (indicated by a line in Fig. 9 connecting point [1,1,0
g/cm2] with the plane) is $>19$ which is still higher than the value estimated
using the Wilks’ theorem in the Likelihood-ratio test for nested model at the
level of 5$\sigma$ ($\Delta\ln\mathcal{L}\approx 16.62$). This confirms that
the modifications of $X_{\text{max}}$ and $R_{\text{had}}$ scales are not an
artifact of the systematic uncertainties of our measurements but a needed
change in the model descriptions. A correction of the results for the biases
seen in the MC-MC tests leads to even larger significance values.

##  IV Discussion

(a)  (b)

Figure 10: The energy evolution of the mean (left) and the standard deviation
(right) of the $X_{\text{max}}$ distributions measured by the Pierre Auger
Observatory using FD [11] (solid circles) and SD [40] (open squares). The
results of this paper for the modified $X_{\text{max}}$ scales (left panel)
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.

###  IV.1 Implications for Inferences on Mass Composition

One straightforward consequence of the $X_{\text{max}}$ shift deeper in the
atmosphere is the solution of the problem with the negative variance of the
logarithmic mass $\sigma^{2}(\ln A)$ derived with QGSJet-II-04 from the
measured $\langle{X_{\text{max}}}\rangle$ and $\sigma(X_{\text{max}})$ (as
discussed in Section I). After application of the corresponding $\Delta
X_{\text{max}}$ shifts, one finds $\sigma^{2}(\ln A)\approx 0.5$ to 2.5 in the
energy range $10^{18.5}\text{ to }10^{19.0}\,\text{eV}$ for all models used in
this work. These values are consistent with the degree of mixing of the
primary composition found in the analysis of the correlation between
$X_{\text{max}}$ and $S(1000)$ with non-modified models, since the correlation
analysis relies on the general phenomenology of air showers and this way is
weakly sensitive to the uncertainties in the description of hadronic
interactions [12, 11].

Another outcome of the method can be foreseen by taking into account the
quasi-universal behavior of the $X_{\text{max}}$ elongation rate for all pure
beams and models with values staying within 54 to 61 g/cm2/decade range.
Changing the elongation rate within these limits introduces an energy-
dependent uncertainty on the MC $X_{\text{max}}$ scale of about $4$
g/cm2/decade at most. Under the assumption that the difference $\Delta
X_{\text{max}}$ between the models and data remains nearly independent of the
primary energy, i.e. that there is no new physics that can significantly
change the predictions for the $X_{\text{max}}$ elongation rate of single
primary species, one could speculate that at the highest energies ($E\gtrsim
10^{19.5}$ eV) the Auger $X_{\text{max}}$ measurements, see Fig. 10, can be
described with a heavy mass composition having a low degree of mixing due to
$\langle{X_{\text{max}}}\rangle$ and $\sigma(X_{\text{max}})$ staying between
the extrapolated predictions of the modified models for oxygen and iron
nuclei.

###  IV.2 Primary Species in the Cosmic-ray Beam

We checked if the shape of the data distribution of ground signal and
$X_{\text{max}}$ and its zenith-angle evolution can be fitted better with an
artificially reduced range of the masses in the MC templates and thus with
different values of $R_{\text{had}}(\theta)$ and $\Delta X_{\text{max}}$. Due
to the presence of deep events and $\sigma(X_{\text{max}})$ values close to
the predictions for protons in $10^{18.5}\text{ to }10^{19.0}\,\text{eV}$
range, we kept protons in the MC templates and used (p, He, O) and (p, He)
mixes for the data fits. In both cases, the quality of data fits was found to
be inferior compared to the fits with (p, He, O, Fe) nuclei (see Table 4).
This was confirmed by the MC-MC tests and was observed using the fits to the
measured data for all three models that the obtained $X_{\text{max}}$ scale
decreases by about 5 to 7, 10 to 17 and 30 to 40 g/cm2 and the hadronic signal
scale by about 2 to 5%, 4 to 9% and 15 to 20% when the heaviest primary Fe is
replaced by Si, O and He in the fit, respectively.

For the five-component (p, He, O, Si, Fe) fits, the values of
$R_{\text{had}}(\theta)$ and $\Delta X_{\text{max}}$ remain well within
statistical errors from the values obtained from the (p, He, O, Fe) fits.
Though the fractions of silicon and iron nuclei are strongly anticorrelated in
such fits, the silicon fraction remains low, ${<}5$%.

Table 4: Minimum values of the log-likelihood expression (see Eq. 3) for different number of primaries assumed to be present in the model predictions. In all cases, the significance of improvement of data description with the $R_{\text{had}}(\theta)$, $\Delta X_{\text{max}}$ and p, He, O, Fe fit using the Likelihood-ratio test applying the Wilks’ theorem for nested model is above $5\sigma$. $\ln\mathcal{L}_{\text{min}}$ | Epos-LHC |  QGSJet-II-04 |  Sibyll 2.3d
---|---|---|---
p He | 518.3 | 633.5 | 563.5
p He O | 467.5 | 523.3 | 486.6
p He O Fe | 451.9 | 476.3 | 451.6

###  IV.3 Modifications of Hadronic Interactions

Assuming the same modifications of hadronic interaction features as in Ref.
[41], the observed increase of the $X_{\text{max}}$ scale in MC predictions
could be explained by an increase of the elasticity or a decrease of the
multiplicity and cross-section. A detailed study of the combinations of such
modifications at the energy equivalent to this study is ongoing [42].

The elasticity is a very good candidate for a potential source of
$X_{\text{max}}$ modification because there are no precise data to constrain
this in models at high energy. It could be precisely measured only in the
region where no detector at the LHC experiments exists. Consequently, the
different models have relatively different predictions with large
uncertainties.

A reduction of the multiplicity can also increase $X_{\text{max}}$ since the
available energy is shared between fewer particles leading to higher
$\uppi^{0}$ energy at the first interaction. But at the same time, the same
effect will reduce the energy available for muon production effectively
increasing the tension in $R_{\text{had}}$.

The superposition model [43] makes an ad-hoc modification of p–p cross-section
that would explain a single change of $X_{\text{max}}$ scale for all primaries
rather difficult. The change in the case of iron nuclei would correspond to a
change of p–air cross-section at the energy 56-times smaller,
${\sim}9{\times}10^{16}$ eV, which starts to be in tension with the
corresponding LHC measurements of p–p collisions [44]. Remarkably, the most
recent cross-section measurements by the ALFA experiment are lower and more
precise [45] than the one used to tune the models, possibly indicating an
overestimation of the current cross-section of the air-shower simulations.

In the case of the multiplicity and elasticity, we can consider the model of
shower development from Ref. [23] also assuming the energy increase of the
total multiplicity as $N\propto N_{0}\,E^{\alpha}$ and a decrease of the
elasticity with energy $\kappa\propto\kappa_{0}\,E^{-\omega}$. We see that an
ad-hoc change of the normalization of multiplicity ($N_{0}$) or elasticity
($\kappa_{0}$) would modify the $X_{\text{max}}$ independently of the primary
and energy,

$\displaystyle X_{\text{max}}^{A}$
$\displaystyle=X_{1}^{A}+X_{0}\ln\frac{\kappa E}{A\,2N\xi^{\uppi}_{\text{c}}}$
(9)
$\displaystyle=X_{1}^{A}+X_{0}\left[(1-\alpha-\omega)\ln\frac{E}{A\,\xi^{\uppi}_{\text{c}}}+\ln\frac{\kappa_{0}}{N_{0}}\right],$

where $\xi^{\uppi}_{\text{c}}$ is the critical energy of pions at which their
decay and interaction lengths are equal.

The hadronic signal can be increased by increasing the multiplicity or
decreasing the ratio of neutral pions according to Ref. [41]. But, as
previously described, an increase in the multiplicity would decrease
$X_{\text{max}}$ increasing the tension there. A change in the ratios of
different hadrons, in particular with more strange particles, is a more likely
possible explanation according to recent LHC data [46, 47].

The energy spectrum of muons naturally influences the dependence of
$S_{\text{had}}$ on $X_{\text{atm}}-X_{\text{max}}$. Though the systematic
uncertainty of the energy scale limits a significant conclusion deduced from
the method (see
$R_{\text{had}}(\theta_{\text{min}})-R_{\text{had}}(\theta_{\text{max}})$ in
Appendix E), the QGSJet-II-04 model seems to predict too hard spectra of
muons. For instance, a larger fraction of strange particles like kaons in the
shower would lead to harder muon spectra because of the larger critical energy
(smaller lifetime) of strange particles.

There are of course other possible modifications of hadronic interactions that
could influence the observed differences between the predictions from the
models and the measured data like cross-sections of low-mass mesons, energy
spectra of pions etc.; see e.g. Ref. [5].

###  IV.4 Limitations of the Method

There are remaining differences between the models and, correspondingly, with
the data due to various limitations of the method. Though our approach leads
to a reduction of the differences between models in
$\langle{X_{\text{max}}}\rangle$, the obtained modifications $\Delta
X_{\text{max}}$ do not cancel out these differences completely. In particular,
the modified $X_{\text{max}}$ scale for Epos-LHC is shallower compared to
QGSJet-II-04 and Sibyll 2.3d (see the left panel of Fig. 10). We found, that a
large part of this difference (${\sim}10$ g/cm2) can be removed if an
additional smearing of $X_{\text{max}}$ is applied to Epos-LHC showers to
compensate for the smaller $X_{\text{max}}$ fluctuations predicted by this
model in comparison to the other two (possibly partly due to strong
defragmentation of nuclei in Epos-LHC [48], as recently confirmed by newer
model EPOS-LHCR [49]). The remaining difference of about 5 g/cm2 between the
models might be due to the statistical errors on $\Delta X_{\text{max}}$ and
differences between the models in the separation of the primary species in
$\langle{X_{\text{max}}}\rangle$.

Possible dependencies of $\Delta X_{\text{max}}$, $R_{\text{had}}(\theta)$,
and fluctuations of $X_{\text{max}}$ and $S(1000)$ on the primary mass are out
of the scope of this paper. We checked that adding a linear mass-dependence of
$\Delta X_{\text{max}}$ into the method did not improve the fit significantly.
The assumption of a mass-independent $\Delta X_{\text{max}}$ used in this
paper was mainly motivated by the similar differences in predicted
$\langle{X_{\text{max}}}\rangle$ for different primaries [49]. Such a quasi-
universal difference is a consequence of the very similar energy dependencies
of interaction features like multiplicity, elasticity, and cross-section
assumed in the models [50]. It means that an offset in these features would
lead to an approximately mass-independent difference in the $X_{\text{max}}$
predictions. An ad-hoc modification of the energy dependence of these
features, like in [41, 42], would lead to mass-dependent $X_{\text{max}}$
shift.

##  V Conclusions

In this paper, we tested the predictions of the models QGSJet-II-04, Sibyll
2.3d, Epos-LHC regarding the depths of maximum of air-shower profiles
$X_{\text{max}}$ and the signal produced by air-shower particles in water-
Cherenkov stations at 1000 m from the shower core, $S(1000)$, composed of
electromagnetic and hadronic ($S_{\text{had}}$) parts. The test consisted of
fitting with MC templates of two-dimensional distributions of ($S(1000)$,
$X_{\text{max}}$) measured at the Pierre Auger Observatory and obtaining the
scales of $X_{\text{max}}$ and $S_{\text{had}}(\theta)$ predicted by the
models, as well as the fractions of four primary nuclei (p, He, O, Fe).

We found that for the best description of the data distributions in the energy
range $10^{18.5}\text{ to }10^{19.0}\,\text{eV}$ for $\theta<60^{\circ}$ the
MC predictions of $X_{\text{max}}$ should be deeper in the atmosphere by about
20 to 50 g/cm2, and the hadronic signal should be increased by about 15 to 25%
in all three models. These modifications reduce the differences between the
models in $X_{\text{max}}$ and $S(1000)$, and as a consequence, lead to
smaller uncertainties on the estimated fractions of the primary nuclei. Due to
the deeper MC $X_{\text{max}}$ scale and, correspondingly, a heavier mass
composition inferred from the data compared with non-modified models, the
scaling factors for the hadronic signal are found to be smaller than in
previous estimations not considering any modifications to the MC
$X_{\text{max}}$ scales. The statistical significance of the improvement in
the data description using the assumed modifications to the MC templates is
above $5\sigma$ for all three models even accounting for all possible linear
combinations of experimental systematic uncertainties. The difference in
$R_{\text{had}}(\theta)$ at the two extreme zenith angles implies an
indication that softer spectra of muons generated by QGSJet-II-04 in 1000 m
from the shower core would describe the Auger data better.

The specific ways to produce the required changes in the models might consist
in combinations of modifications of integral (cross-section, multiplicity,
elasticity etc.) or differential (secondary particles energy spectra)
characteristics of the hadronic interactions as discussed e.g. in Refs. [41,
5, 42]. Our method addresses only the first-order differences in the mean
values of $X_{\text{max}}$ and hadronic signal without taking into account
their possible dependencies on the primary mass or energy. These dependencies,
as well as the investigation of a modification of fluctuations of air-shower
observables, need to be studied further to corroborate the scale adjustments
found in this paper, supplemented also by future data with increased mass-
sensitivity from AugerPrime [51], the upgrade of the Pierre Auger Observatory.

###### Acknowledgements.

The successful installation, commissioning, and operation of the Pierre Auger
Observatory would not have been possible without the strong commitment and
effort from the technical and administrative staff in Malargüe. We are very
grateful to the following agencies and organizations for financial support:
Argentina – Comisión Nacional de Energía Atómica; Agencia Nacional de
Promoción Científica y Tecnológica (ANPCyT); Consejo Nacional de
Investigaciones Científicas y Técnicas (CONICET); Gobierno de la Provincia de
Mendoza; Municipalidad de Malargüe; NDM Holdings and Valle Las Leñas; in
gratitude for their continuing cooperation over land access; Australia – the
Australian Research Council; Belgium – Fonds de la Recherche Scientifique
(FNRS); Research Foundation Flanders (FWO), Marie Curie Action of the European
Union Grant No. 101107047; Brazil – Conselho Nacional de Desenvolvimento
Científico e Tecnológico (CNPq); Financiadora de Estudos e Projetos (FINEP);
Fundação de Amparo à Pesquisa do Estado de Rio de Janeiro (FAPERJ); São Paulo
Research Foundation (FAPESP) Grants No. 2019/10151-2, No. 2010/07359-6 and No.
1999/05404-3; Ministério da Ciência, Tecnologia, Inovações e Comunicações
(MCTIC); Czech Republic – GACR 24-13049S, CAS LQ100102401, MEYS LM2023032,
CZ.02.1.01/0.0/0.0/16_013/0001402, CZ.02.1.01/0.0/0.0/18_046/0016010 and
CZ.02.1.01/0.0/0.0/17_049/0008422 and CZ.02.01.01/00/22_008/0004632; France –
Centre de Calcul IN2P3/CNRS; Centre National de la Recherche Scientifique
(CNRS); Conseil Régional Ile-de-France; Département Physique Nucléaire et
Corpusculaire (PNC-IN2P3/CNRS); Département Sciences de l’Univers (SDU-
INSU/CNRS); Institut Lagrange de Paris (ILP) Grant No. LABEX ANR-10-LABX-63
within the Investissements d’Avenir Programme Grant No. ANR-11-IDEX-0004-02;
Germany – Bundesministerium für Bildung und Forschung (BMBF); Deutsche
Forschungsgemeinschaft (DFG); Finanzministerium Baden-Württemberg; Helmholtz
Alliance for Astroparticle Physics (HAP); Helmholtz-Gemeinschaft Deutscher
Forschungszentren (HGF); Ministerium für Kultur und Wissenschaft des Landes
Nordrhein-Westfalen; Ministerium für Wissenschaft, Forschung und Kunst des
Landes Baden-Württemberg; Italy – Istituto Nazionale di Fisica Nucleare
(INFN); Istituto Nazionale di Astrofisica (INAF); Ministero dell’Università e
della Ricerca (MUR); CETEMPS Center of Excellence; Ministero degli Affari
Esteri (MAE), ICSC Centro Nazionale di Ricerca in High Performance Computing,
Big Data and Quantum Computing, funded by European Union NextGenerationEU,
reference code CN_00000013; México – Consejo Nacional de Ciencia y Tecnología
(CONACYT) No. 167733; Universidad Nacional Autónoma de México (UNAM); PAPIIT
DGAPA-UNAM; The Netherlands – Ministry of Education, Culture and Science;
Netherlands Organisation for Scientific Research (NWO); Dutch national
e-infrastructure with the support of SURF Cooperative; Poland – Ministry of
Education and Science, grants No. DIR/WK/2018/11 and 2022/WK/12; National
Science Centre, grants No. 2016/22/M/ST9/00198, 2016/23/B/ST9/01635,
2020/39/B/ST9/01398, and 2022/45/B/ST9/02163; Portugal – Portuguese national
funds and FEDER funds within Programa Operacional Factores de Competitividade
through Fundação para a Ciência e a Tecnologia (COMPETE); Romania – Ministry
of Research, Innovation and Digitization, CNCS-UEFISCDI, contract no. 30N/2023
under Romanian National Core Program LAPLAS VII, grant no. PN 23 21 01 02 and
project number PN-III-P1-1.1-TE-2021-0924/TE57/2022, within PNCDI III;
Slovenia – Slovenian Research Agency, grants P1-0031, P1-0385, I0-0033,
N1-0111; Spain – Ministerio de Economía, Industria y Competitividad
(FPA2017-85114-P and PID2019-104676GB-C32), Xunta de Galicia (ED431C 2017/07),
Junta de Andalucía (SOMM17/6104/UGR, P18-FR-4314) Feder Funds, RENATA Red
Nacional Temática de Astropartículas (FPA2015-68783-REDT) and María de Maeztu
Unit of Excellence (MDM-2016-0692); USA – Department of Energy, Contracts No.
DE-AC02-07CH11359, No. DE-FR02-04ER41300, No. DE-FG02-99ER41107 and No. DE-
SC0011689; National Science Foundation, Grant No. 0450696; The Grainger
Foundation; Marie Curie-IRSES/EPLANET; European Particle Physics Latin
American Network; and UNESCO.

## References

  * Pierog  _et al._ [2015] T. Pierog, I. Karpenko, J. M. Katzy, E. Yatsenko, and K. Werner, EPOS LHC: Test of collective hadronization with data measured at the cern large hadron collider, Phys. Rev. C 92, 034906 (2015).
  * Ostapchenko [2011] S. Ostapchenko, Monte Carlo treatment of hadronic interactions in enhanced Pomeron scheme: QGSJET-II model, Phys. Rev. D 83, 014018 (2011).
  * Riehn  _et al._ [2020] F. Riehn, R. Engel, A. Fedynitch, T. K. Gaisser, and T. Stanev, Hadronic interaction model Sibyll 2.3d and extensive air showers, Phys. Rev. D 102, 063002 (2020).
  * Ostapchenko [2016] S. Ostapchenko, LHC results and hadronic interaction models, ArXiv e-prints 1612.09461 [astro-ph.HE] (2016), arXiv:1612.09461 [astro-ph.HE] .
  * Albrecht  _et al._ [2022] J. Albrecht _et al._ , The Muon Puzzle in cosmic-ray induced air showers and its connection to the Large Hadron Collider, Astrophysics and Space Science 367, 27 (2022).
  * Aab A. et al. [2015] Aab A. et al. (Pierre Auger Collaboration), The Pierre Auger Cosmic Ray Observatory, Nucl. Instrum. Methods Phys. Res. A 798, 172 (2015).
  * Aab  _et al._ [2020a] A. Aab _et al._ (Pierre Auger Collaboration), Reconstruction of events recorded with the surface detector of the Pierre Auger Observatory, Journal of Instrumentation 15 (10), P10021.
  * Abreu  _et al._ [2013] P. Abreu _et al._ (Pierre Auger Collaboration), Interpretation of the depths of maximum of extensive air showers measured by the Pierre Auger Observatory, Journal of Cosmology and Astroparticle Physics 2013 (02), 026.
  * Aab  _et al._ [2014] A. Aab _et al._ (Pierre Auger Collaboration), Depth of maximum of air-shower profiles at the Pierre Auger Observatory. I. Measurements at energies above $10^{17.8}$ eV, Phys. Rev. D 90, 122005 (2014).
  * Aab  _et al._ [2014a] A. Aab _et al._ (Pierre Auger Collaboration), Depth of maximum of air-shower profiles at the Pierre Auger Observatory. II. Composition implications, Phys. Rev. D 90, 122006 (2014a).
  * Yushkov [2019] A. Yushkov (Pierre Auger Collaboration), Mass composition of cosmic rays with energies above $10^{17.2}$ eV from the hybrid data of the Pierre Auger Observatory,  _Proceedings, 36th International Cosmic Ray Conference (ICRC 2019): Madison, WI, U.S.A., July 24th - August 1st 2019_ , PoS ICRC2019, 482 (2019).
  * Aab  _et al._ [2016a] A. Aab _et al._ (Pierre Auger Collaboration), Evidence for a mixed mass composition at the ‘ankle’ in the cosmic-ray spectrum, Physics Letters B 762, 288 (2016a).
  * Aab  _et al._ [2015] A. Aab _et al._ (Pierre Auger Collaboration), Muons in air showers at the Pierre Auger Observatory: Mean number in highly inclined events, Phys. Rev. D 91, 032003 (2015).
  * Aab  _et al._ [2021] A. Aab _et al._ (Pierre Auger Collaboration), Measurement of the fluctuations in the number of muons in extensive air showers with the pierre auger observatory, Phys. Rev. Lett. 126, 152002 (2021).
  * Aab  _et al._ [2020b] A. Aab _et al._ (Pierre Auger Collaboration), Direct measurement of the muonic content of extensive air showers between $2\cdot 10^{17}$ and $2\cdot 10^{18}$ eV at the Pierre Auger Observatory, Eur. Phys. J. C 210, 751 (2020b).
  * Aab  _et al._ [2014b] A. Aab _et al._ (Pierre Auger Collaboration), Muons in air showers at the pierre auger observatory: Measurement of atmospheric production depth, Phys. Rev. D 90, 012012 (2014b).
  * Aab  _et al._ [2016b] A. Aab _et al._ (Pierre Auger Collaboration), Testing Hadronic Interactions at Ultrahigh Energies with Air Showers Measured by the Pierre Auger Observatory, Phys. Rev. Lett. 117, 192001 (2016b).
  * Matthews [2005] J. Matthews, A Heitler model of extensive air showers, Astroparticle Physics 22, 387 (2005).
  * Vícha [2016] J. Vícha,  _Analysis of Air Showers with respect to Primary Composition of Cosmic Rays_ (PhD thesis, Czech Technical University, Prague, 2016) pp. 83 – 114.
  * Vícha  _et al._ [2019] J. Vícha _et al._ , Testing Hadronic Interactions Using Hybrid Observables,  _Proceedings, 36th International Cosmic Ray Conference (ICRC 2019): Madison, WI, U.S.A., July 24th - August 1st 2019_ , PoS ICRC2019, 452 (2019).
  * Vícha [2021] J. Vícha (Pierre Auger Collaboration), Adjustments to Model Predictions of Depth of Shower Maximum and Signals at Ground Level using Hybrid Events of the Pierre Auger Observatory,  _Proceedings, 37th International Cosmic Ray Conference (ICRC 2021): Berlin, Germany, July 12th - 23rd 2021_ , PoS ICRC2021, 310 (2021).
  * Aab  _et al._ [2020c] A. Aab _et al._ (Pierre Auger Collaboration), Measurement of the cosmic-ray energy spectrum above $2.5\times{}{10}^{18}\text{ }\text{ }\mathrm{eV}$ using the pierre auger observatory, Phys. Rev. D 102, 062005 (2020c).
  * Kampert and Unger [2012] K.-H. Kampert and M. Unger, Measurements of the cosmic ray composition with air shower experiments, Astroparticle Physics 35, 660 (2012).
  * Ave  _et al._ [2017a] M. Ave, R. Engel, M. Roth, and A. Schulz, A generalized description of the signal size in extensive air shower detectors and its applications, Astroparticle Physics 87, 23 (2017a).
  * Ave  _et al._ [2017b] M. Ave, M. Roth, and A. Schulz, A generalized description of the time dependent signals in extensive air shower detectors and its applications, Astropart. Phys. 88, 46 (2017b).
  * dos Santos [2021] E. dos Santos (Pierre Auger Collaboration), Monte Carlo simulations for the Pierre Auger Observatory using the VO auger grid resources,  _Proceedings, 37th International Cosmic Ray Conference (ICRC 2021): Berlin, Germany, July 12th - 23rd 2021_ , PoS ICRC2021, 232 (2021).
  * Heck  _et al._ [1998] D. Heck _et al._ , CORSIKA: A Monte Carlo code to simulate extensive air showers, Report FZKA Forschungszentrum Karlsruhe 6019 (1998).
  * Heck and Knapp [1998] D. Heck and J. Knapp, Upgrade of the Monte Carlo code CORSIKA to simulate extensive air showers with energies $>10^{20}$ eV, Report FZKA Forschungszentrum Karlsruhe 6097 (1998).
  * Heck [2006] D. Heck, Air shower simulation with CORSIKA at arbitrary direction of incidence, Report FZKA Forschungszentrum Karlsruhe 7254 (2006).
  * Bleicher  _et al._ [1999] M. Bleicher, E. Zabrodin, C. Spieles, S. A. Bass, C. Ernst, S. Soff, L. Bravina, M. Belkacem, H. Weber, H. Stöcker, and W. Greiner, Relativistic hadron-hadron collisions in the ultra-relativistic quantum molecular dynamics model, Journal of Physics G: Nuclear and Particle Physics 25, 1859 (1999).
  * Battistoni  _et al._ [2015] G. Battistoni _et al._ , Overview of the FLUKA code, Annals Nucl. Energy 82, 10 (2015).
  * Böhlen  _et al._ [2014] T. Böhlen, F. Cerutti, M. Chin, A. Fassò, A. Ferrari, P. Ortega, A. Mairani, P. Sala, G. Smirnov, and V. Vlachoudis, The fluka code: Developments and challenges for high energy and medical applications, Nuclear Data Sheets 120, 211 (2014).
  * Argirò  _et al._ [2007] S. Argirò, S. Barroso, J. Gonzalez, L. Nellen, T. Paul, T. Porter, L. Prado Jr., M. Roth, R. Ulrich, and D. Veberič, The offline software framework of the pierre auger observatory, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 580, 1485 (2007).
  * Wilks [1938] S. S. Wilks, The Large-Sample Distribution of the Likelihood Ratio for Testing Composite Hypotheses, Annals Math. Statist. 9, 60 (1938).
  * Fan  _et al._ [2001] J. Fan, C. Zhang, and J. Zhang, Generalized Likelihood Ratio Statistics and Wilks Phenomenon, The Annals of Statistics 29, 153 (2001).
  * Choi [2021] K. Choi (Pierre Auger Collaboration), Long Term Performance of the Pierre Auger Observatory,  _Proceedings, 36th International Cosmic Ray Conference (ICRC 2019): Madison, WI, U.S.A., July 24th - August 1st 2019_ , PoS ICRC2019, 222 (2021).
  * Aab  _et al._ [2017] A. Aab _et al._ (Pierre Auger Collaboration), Impact of Atmospheric Effects on the Energy Reconstruction of Air Showers Observed by the Surface Detectors of the Pierre Auger Observatory, JINST 12 (02), P02006, arXiv:1702.02835 [astro-ph.IM] .
  * Gideon and Hollister [1987] R. A. Gideon and R. A. Hollister, A rank correlation coefficient resistant to outliers, Journal of the American Statistical Association 82, 656 (1987).
  * J. Bellido for the Pierre Auger Collaboration [2018] J. Bellido for the Pierre Auger Collaboration, Depth of maximum of air-shower profiles at the Pierre Auger Observatory: Measurements above 1017.2 eV and Composition Implications, in  _Proceedings of the 35th International Cosmic Ray Conference ICRC 2017, PoS(ICRC2017)506_ (2018).
  * Peixoto [2019] C. J. T. Peixoto (Pierre Auger Collaboration), Estimating the Depth of Shower Maximum using the Surface Detectors of the Pierre Auger Observatory,  _Proceedings, 36th International Cosmic Ray Conference (ICRC 2019): Madison, WI, U.S.A., July 24th - August 1st 2019_ , PoS ICRC2019, 440 (2019).
  * Ulrich  _et al._ [2011] R. Ulrich, R. Engel, and M. Unger, Hadronic multiparticle production at ultrahigh energies and extensive air showers, Phys. Rev. D 83, 054026 (2011).
  * Ebr  _et al._ [2023] J. Ebr _et al._ , Impact of modified characteristics of hadronic interactions on cosmic-ray observables for proton and nuclear primaries,  _Proceedings, 38th International Cosmic Ray Conference (ICRC 2023): Nagoya, Japan, July 26th - August 3rd 2023_ , PoS ICRC2023, 245 (2023).
  * Engel  _et al._ [1992] J. Engel, T. K. Gaisser, P. Lipari, and T. Stanev, Nucleus-nucleus collisions and interpretation of cosmic-ray cascades, Phys. Rev. D 46, 5013 (1992).
  * Antchev  _et al._ [2019] G. Antchev _et al._ (TOTEM Collaboration), First measurement of elastic, inelastic and total cross-section at $\sqrt{s}=13$ tev by totem and overview of cross-section data at lhc energies, European Physical Journal C 79 (2019).
  * Grafstrom [2015] P. Grafstrom (ATLAS Collaboration), Measurements of the elastic, inelastic and total cross sections in pp collisions with ATLAS sub-detectors, PoS EPS-HEP2015, 486 (2015).
  * Baur  _et al._ [2023] S. Baur, H. Dembinski, M. Perlin, T. Pierog, R. Ulrich, and K. Werner, Core-corona effect in hadron collisions and muon production in air showers, Phys. Rev. D 107, 094031 (2023).
  * Anchordoqui  _et al._ [2022] L. A. Anchordoqui, C. G. Canal, F. Kling, S. J. Sciutto, and J. F. Soriano, An explanation of the muon puzzle of ultrahigh-energy cosmic rays and the role of the Forward Physics Facility for model improvement, JHEAp 34, 19 (2022), arXiv:2202.03095 [hep-ph] .
  * Ostapchenko [2018] S. Ostapchenko,  _High Energy Cosmic Ray Interactions and UHECR Composition Problem, talk at the UHECR 2018, Paris, October 08-12, 2018 (Paris)_ (2018 (accessed Oct, 2018)).
  * Werner [2023] T. P. . K. Werner, EPOS LHC-R : up-to-date hadronic model for EAS simulations,  _Proceedings, 38th International Cosmic Ray Conference (ICRC 2023): Nagoya, Japan, July 26th - August 3rd 2023_ , PoS ICRC2023, 230 (2023).
  * Alvarez-Muñiz  _et al._ [2002] J. Alvarez-Muñiz, R. Engel, T. K. Gaisser, J. A. Ortiz, and T. Stanev, Hybrid simulations of extensive air showers, Phys. Rev. D 66, 033011 (2002).
  * Castellina [2019] A. Castellina (Pierre Auger Collaboration), AugerPrime: the Pierre Auger Observatory Upgrade, EPJ Web Conf. 210, 06002 (2019), arXiv:1905.04472 [astro-ph.HE] .
  * Domenico  _et al._ [2013] M. D. Domenico, M. Settimo, S. Riggi, and E. Bertin, Reinterpreting the development of extensive air showers initiated by nuclei and photons, Journal of Cosmology and Astroparticle Physics 2013 (07), 050.
  * Justel  _et al._ [1997] A. Justel, D. Pena, and R. Zamar, A multivariate Kolmogorov-Smirnov test of goodness of fit, Statistics & Probability Letters 35, 251 (1997).
  * Gaisser and Hillas [1977] T. K. Gaisser and A. M. Hillas, Reliability of the Method of Constant Intensity Cuts for Reconstructing the Average Development of Vertical Showers, in  _International Cosmic Ray Conference_ , International Cosmic Ray Conference, Vol. 8 (1977) p. 353.

$\bullet$

A. Abdul Halim13, P. Abreu73, M. Aglietta55,53, I. Allekotte1, K. Almeida
Cheminant71, A. Almela7,12, R. Aloisio46,47, J. Alvarez-Muñiz79, J. Ammerman
Yebra79, G.A. Anastasi59,48, L. Anchordoqui86, B. Andrada7, S. Andringa73, L.
Apollonio60,50, C. Aramo51, P.R. Araújo Ferreira43, E. Arnone64,53, J.C.
Arteaga Velázquez68, P. Assis73, G. Avila11, E. Avocone58,47, A. Bakalova33,
F. Barbato46,47, A. Bartz Mocellin85, J.A. Bellido13,70, C. Berat37, M.E.
Bertaina64,53, G. Bhatta71, M. Bianciotto64,53, P.L. Biermanni, V. Binet5, K.
Bismark40,7, T. Bister80,81, J. Biteau38,b, J. Blazek33, C. Bleve37, J.
Blümer42, M. Boháčová33, D. Boncioli58,47, C. Bonifazi8,27, L. Bonneau
Arbeletche22, N. Borodai71, J. Brackk, P.G. Brichetto Orchera7, F.L.
Briechle43, A. Bueno78, S. Buitink15, M. Buscemi48,62, M. Büsken40,7, A.
Bwembya80,81, K.S. Caballero-Mora67, S. Cabana-Freire79, L. Caccianiga60,50,
F. Campuzano6, R. Caruso59,48, A. Castellina55,53, F. Catalani19, G.
Cataldi49, L. Cazon79, M. Cerda10, A. Cermenati46,47, J.A. Chinellato22, J.
Chudoba33, L. Chytka34, R.W. Clay13, A.C. Cobos Cerutti6, R. Colalillo61,51,
M.R. Coluccia49, R. Conceição73, A. Condorelli38, G. Consolati50,56, M.
Conte57,49, F. Convenga58,47, D. Correia dos Santos29, P.J. Costa73, C.E.
Covault84, M. Cristinziani45, C.S. Cruz Sanchez3, S. Dasso4,2, K. Daumiller42,
B.R. Dawson13, R.M. de Almeida29, J. de Jesús7,42, S.J. de Jong80,81, J.R.T.
de Mello Neto27,28, I. De Mitri46,47, J. de Oliveira18, D. de Oliveira
Franco49, F. de Palma57,49, V. de Souza20, B.P. de Souza de Errico27, E. De
Vito57,49, A. Del Popolo59,48, O. Deligny35, N. Denner33, L. Deval42,7, A. di
Matteo53, M. Dobre74, C. Dobrigkeit22, J.C. D’Olivo69, L.M. Domingues
Mendes73,16, Q. Dorosti45, J.C. dos Anjos16, R.C. dos Anjos26, J. Ebr33, F.
Ellwanger42, M. Emam80,81, R. Engel40,42, I. Epicoco57,49, M. Erdmann43, A.
Etchegoyen7,12, C. Evoli46,47, H. Falcke80,82,81, G. Farrar88, A.C. Fauth22,
N. Fazzinig, F. Feldbusch41, F. Fenu42,f, A. Fernandes73, B. Fick87, J.M.
Figueira7, A. Filipčič77,76, T. Fitoussi42, B. Flaggs90, T. Fodran80, T.
Fujii89,h, A. Fuster7,12, C. Galea80, C. Galelli60,50, B. García6, C. Gaudu39,
H. Gemmeke41, F. Gesualdi7,42, A. Gherghel-Lascu74, P.L. Ghia35, U.
Giaccari49, J. Glombitza43,d, F. Gobbi10, F. Gollan7, G. Golup1, M. Gómez
Berisso1, P.F. Gómez Vitale11, J.P. Gongora11, J.M. González1, N. González7,
D. Góra71, A. Gorgi55,53, M. Gottowik79, T.D. Grubb13, F. Guarino61,51, G.P.
Guedes23, E. Guido45, L. Gülzow42, S. Hahn40, P. Hamal33, M.R. Hampel7, P.
Hansen3, D. Harari1, V.M. Harvey13, A. Haungs42, T. Hebbeker43, C. Hojvatg,
J.R. Hörandel80,81, P. Horvath34, M. Hrabovský34, T. Huege42,15, A.
Insolia59,48, P.G. Isar75, P. Janecek33, V. Jilek33, J.A. Johnsen85, J.
Jurysek33, K.-H. Kampert39, B. Keilhauer42, A. Khakurdikar80, V.V. Kizakke
Covilakam7,42, H.O. Klages42, M. Kleifges41, F. Knapp40, J. Köhler42, N.
Kunka41, B.L. Lago17, N. Langner43, M.A. Leigui de Oliveira25, Y. Lema-
Capeans79, A. Letessier-Selvon36, I. Lhenry-Yvon35, L. Lopes73, L. Lu91, Q.
Luce40, J.P. Lundquist76, A. Machado Payeras22, M. Majercakova33, D. Mandat33,
B.C. Manning13, P. Mantschg, F.M. Mariani60,50, A.G. Mariazzi3, I.C. Mariş14,
G. Marsella62,48, D. Martello57,49, S. Martinelli42,7, O. Martínez Bravo65,
M.A. Martins79, H.-J. Mathes42, J. Matthewsa, G. Matthiae63,52, E.
Mayotte85,39, S. Mayotte85, P.O. Mazurg, G. Medina-Tanco69, J. Meinert39, D.
Melo7, A. Menshikov41, C. Merx42, S. Michal33, M.I. Micheletti5, L.
Miramonti60,50, S. Mollerach1, F. Montanet37, L. Morejon39, C. Morello55,53,
K. Mulrey80,81, R. Mussa53, W.M. Namasaka39, S. Negi33, L. Nellen69, K.
Nguyen87, G. Nicora9, M. Niechciol45, D. Nitz87, D. Nosek32, V. Novotny32, L.
Nožka34, A. Nucita57,49, L.A. Núñez31, C. Oliveira20, M. Palatka33, J.
Pallotta9, S. Panja33, G. Parente79, T. Paulsen39, J. Pawlowsky39, M. Pech33,
J. Pȩkala71, R. Pelayo66, L.A.S. Pereira24, E.E. Pereira Martins40,7, J. Perez
Armand21, C. Pérez Bertolli7,42, L. Perrone57,49, S. Petrera46,47, C.
Petrucci58,47, T. Pierog42, M. Pimenta73, M. Platino7, B. Pont80, M.
Pothast81,80, M. Pourmohammad Shahvar62,48, P. Privitera89, M. Prouza33, S.
Querchfeld39, J. Rautenberg39, D. Ravignani7, J.V. Reginatto Akim22, M.
Reininghaus40, J. Ridky33, F. Riehn79, M. Risse45, V. Rizi58,47, W. Rodrigues
de Carvalho80, E. Rodriguez7,42, J. Rodriguez Rojo11, M.J. Roncoroni7, S.
Rossoni44, M. Roth42, E. Roulet1, A.C. Rovero4, P. Ruehl45, A. Saftoiu74, M.
Saharan80, F. Salamida58,47, H. Salazar65, G. Salina52, J.D. Sanabria Gomez31,
F. Sánchez7, E.M. Santos21, E. Santos33, F. Sarazin85, R. Sarmento73, R.
Sato11, P. Savina91, C.M. Schäfer40, V. Scherini57,49, H. Schieler42, M.
Schimassek35, M. Schimp39, D. Schmidt42, O. Scholten15,j, H.
Schoorlemmer80,81, P. Schovánek33, F.G. Schröder90,42, J. Schulte43, T.
Schulz42, S.J. Sciutto3, M. Scornavacche7,42, A. Sedoski7, A. Segreto54,48, S.
Sehgal39, S.U. Shivashankara76, G. Sigl44, G. Silli7, O. Sima74,c, K.
Simkova15, F. Simon41, R. Smau74, R. Šmída89, P. Sommersl, J.F. Soriano86, R.
Squartini10, M. Stadelmaier50,60,42, S. Stanič76, J. Stasielak71, P. Stassi37,
S. Strähnz40, M. Straub43, T. Suomijärvi38, A.D. Supanitsky7, Z.
Svozilikova33, Z. Szadkowski72, F. Tairli13, A. Tapia30, C. Taricco64,53, C.
Timmermans81,80, O. Tkachenko42, P. Tobiska33, C.J. Todero Peixoto19, B.
Tomé73, Z. Torrès37, A. Travaini10, P. Travnicek33, C. Trimarelli58,47, M.
Tueros3, M. Unger42, L. Vaclavek34, M. Vacula34, J.F. Valdés Galicia69, L.
Valore61,51, E. Varela65, A. Vásquez-Ramírez31, D. Veberič42, C. Ventura28,
I.D. Vergara Quispe3, V. Verzi52, J. Vicha33, J. Vink83, S. Vorobiov76, C.
Watanabe27, A.A. Watsone, A. Weindl42, L. Wiencke85, H. Wilczyński71, D.
Wittkowski39, B. Wundheiler7, B. Yue39, A. Yushkov33, O. Zapparrata14, E.
Zas79, D. Zavrtanik76,77, M. Zavrtanik77,76

The Pierre Auger Collaboration

1

Centro Atómico Bariloche and Instituto Balseiro (CNEA-UNCuyo-CONICET), San
Carlos de Bariloche, Argentina

2

Departamento de Física and Departamento de Ciencias de la Atmósfera y los
Océanos, FCEyN, Universidad de Buenos Aires and CONICET, Buenos Aires,
Argentina

3

IFLP, Universidad Nacional de La Plata and CONICET, La Plata, Argentina

4

Instituto de Astronomía y Física del Espacio (IAFE, CONICET-UBA), Buenos
Aires, Argentina

5

Instituto de Física de Rosario (IFIR) – CONICET/U.N.R. and Facultad de
Ciencias Bioquímicas y Farmacéuticas U.N.R., Rosario, Argentina

6

Instituto de Tecnologías en Detección y Astropartículas (CNEA, CONICET,
UNSAM), and Universidad Tecnológica Nacional – Facultad Regional Mendoza
(CONICET/CNEA), Mendoza, Argentina

7

Instituto de Tecnologías en Detección y Astropartículas (CNEA, CONICET,
UNSAM), Buenos Aires, Argentina

8

International Center of Advanced Studies and Instituto de Ciencias Físicas,
ECyT-UNSAM and CONICET, Campus Miguelete – San Martín, Buenos Aires, Argentina

9

Laboratorio Atmósfera – Departamento de Investigaciones en Láseres y sus
Aplicaciones – UNIDEF (CITEDEF-CONICET), Argentina

10

Observatorio Pierre Auger, Malargüe, Argentina

11

Observatorio Pierre Auger and Comisión Nacional de Energía Atómica, Malargüe,
Argentina

12

Universidad Tecnológica Nacional – Facultad Regional Buenos Aires, Buenos
Aires, Argentina

13

University of Adelaide, Adelaide, S.A., Australia

14

Université Libre de Bruxelles (ULB), Brussels, Belgium

15

Vrije Universiteit Brussels, Brussels, Belgium

16

Centro Brasileiro de Pesquisas Fisicas, Rio de Janeiro, RJ, Brazil

17

Centro Federal de Educação Tecnológica Celso Suckow da Fonseca, Petropolis,
Brazil

18

Instituto Federal de Educação, Ciência e Tecnologia do Rio de Janeiro (IFRJ),
Brazil

19

Universidade de São Paulo, Escola de Engenharia de Lorena, Lorena, SP, Brazil

20

Universidade de São Paulo, Instituto de Física de São Carlos, São Carlos, SP,
Brazil

21

Universidade de São Paulo, Instituto de Física, São Paulo, SP, Brazil

22

Universidade Estadual de Campinas (UNICAMP), IFGW, Campinas, SP, Brazil

23

Universidade Estadual de Feira de Santana, Feira de Santana, Brazil

24

Universidade Federal de Campina Grande, Centro de Ciencias e Tecnologia,
Campina Grande, Brazil

25

Universidade Federal do ABC, Santo André, SP, Brazil

26

Universidade Federal do Paraná, Setor Palotina, Palotina, Brazil

27

Universidade Federal do Rio de Janeiro, Instituto de Física, Rio de Janeiro,
RJ, Brazil

28

Universidade Federal do Rio de Janeiro (UFRJ), Observatório do Valongo, Rio de
Janeiro, RJ, Brazil

29

Universidade Federal Fluminense, EEIMVR, Volta Redonda, RJ, Brazil

30

Universidad de Medellín, Medellín, Colombia

31

Universidad Industrial de Santander, Bucaramanga, Colombia

32

Charles University, Faculty of Mathematics and Physics, Institute of Particle
and Nuclear Physics, Prague, Czech Republic

33

Institute of Physics of the Czech Academy of Sciences, Prague, Czech Republic

34

Palacky University, Olomouc, Czech Republic

35

CNRS/IN2P3, IJCLab, Université Paris-Saclay, Orsay, France

36

Laboratoire de Physique Nucléaire et de Hautes Energies (LPNHE), Sorbonne
Université, Université de Paris, CNRS-IN2P3, Paris, France

37

Univ. Grenoble Alpes, CNRS, Grenoble Institute of Engineering Univ. Grenoble
Alpes, LPSC-IN2P3, 38000 Grenoble, France

38

Université Paris-Saclay, CNRS/IN2P3, IJCLab, Orsay, France

39

Bergische Universität Wuppertal, Department of Physics, Wuppertal, Germany

40

Karlsruhe Institute of Technology (KIT), Institute for Experimental Particle
Physics, Karlsruhe, Germany

41

Karlsruhe Institute of Technology (KIT), Institut für Prozessdatenverarbeitung
und Elektronik, Karlsruhe, Germany

42

Karlsruhe Institute of Technology (KIT), Institute for Astroparticle Physics,
Karlsruhe, Germany

43

RWTH Aachen University, III. Physikalisches Institut A, Aachen, Germany

44

Universität Hamburg, II. Institut für Theoretische Physik, Hamburg, Germany

45

Universität Siegen, Department Physik – Experimentelle Teilchenphysik, Siegen,
Germany

46

Gran Sasso Science Institute, L’Aquila, Italy

47

INFN Laboratori Nazionali del Gran Sasso, Assergi (L’Aquila), Italy

48

INFN, Sezione di Catania, Catania, Italy

49

INFN, Sezione di Lecce, Lecce, Italy

50

INFN, Sezione di Milano, Milano, Italy

51

INFN, Sezione di Napoli, Napoli, Italy

52

INFN, Sezione di Roma “Tor Vergata”, Roma, Italy

53

INFN, Sezione di Torino, Torino, Italy

54

Istituto di Astrofisica Spaziale e Fisica Cosmica di Palermo (INAF), Palermo,
Italy

55

Osservatorio Astrofisico di Torino (INAF), Torino, Italy

56

Politecnico di Milano, Dipartimento di Scienze e Tecnologie Aerospaziali ,
Milano, Italy

57

Università del Salento, Dipartimento di Matematica e Fisica “E. De Giorgi”,
Lecce, Italy

58

Università dell’Aquila, Dipartimento di Scienze Fisiche e Chimiche, L’Aquila,
Italy

59

Università di Catania, Dipartimento di Fisica e Astronomia “Ettore Majorana“,
Catania, Italy

60

Università di Milano, Dipartimento di Fisica, Milano, Italy

61

Università di Napoli “Federico II”, Dipartimento di Fisica “Ettore Pancini”,
Napoli, Italy

62

Università di Palermo, Dipartimento di Fisica e Chimica ”E. Segrè”, Palermo,
Italy

63

Università di Roma “Tor Vergata”, Dipartimento di Fisica, Roma, Italy

64

Università Torino, Dipartimento di Fisica, Torino, Italy

65

Benemérita Universidad Autónoma de Puebla, Puebla, México

66

Unidad Profesional Interdisciplinaria en Ingeniería y Tecnologías Avanzadas
del Instituto Politécnico Nacional (UPIITA-IPN), México, D.F., México

67

Universidad Autónoma de Chiapas, Tuxtla Gutiérrez, Chiapas, México

68

Universidad Michoacana de San Nicolás de Hidalgo, Morelia, Michoacán, México

69

Universidad Nacional Autónoma de México, México, D.F., México

70

Universidad Nacional de San Agustin de Arequipa, Facultad de Ciencias
Naturales y Formales, Arequipa, Peru

71

Institute of Nuclear Physics PAN, Krakow, Poland

72

University of Łódź, Faculty of High-Energy Astrophysics,Łódź, Poland

73

Laboratório de Instrumentação e Física Experimental de Partículas – LIP and
Instituto Superior Técnico – IST, Universidade de Lisboa – UL, Lisboa,
Portugal

74

“Horia Hulubei” National Institute for Physics and Nuclear Engineering,
Bucharest-Magurele, Romania

75

Institute of Space Science, Bucharest-Magurele, Romania

76

Center for Astrophysics and Cosmology (CAC), University of Nova Gorica, Nova
Gorica, Slovenia

77

Experimental Particle Physics Department, J. Stefan Institute, Ljubljana,
Slovenia

78

Universidad de Granada and C.A.F.P.E., Granada, Spain

79

Instituto Galego de Física de Altas Enerxías (IGFAE), Universidade de Santiago
de Compostela, Santiago de Compostela, Spain

80

IMAPP, Radboud University Nijmegen, Nijmegen, The Netherlands

81

Nationaal Instituut voor Kernfysica en Hoge Energie Fysica (NIKHEF), Science
Park, Amsterdam, The Netherlands

82

Stichting Astronomisch Onderzoek in Nederland (ASTRON), Dwingeloo, The
Netherlands

83

Universiteit van Amsterdam, Faculty of Science, Amsterdam, The Netherlands

84

Case Western Reserve University, Cleveland, OH, USA

85

Colorado School of Mines, Golden, CO, USA

86

Department of Physics and Astronomy, Lehman College, City University of New
York, Bronx, NY, USA

87

Michigan Technological University, Houghton, MI, USA

88

New York University, New York, NY, USA

89

University of Chicago, Enrico Fermi Institute, Chicago, IL, USA

90

University of Delaware, Department of Physics and Astronomy, Bartol Research
Institute, Newark, DE, USA

91

University of Wisconsin-Madison, Department of Physics and WIPAC, Madison, WI,
USA

—–

a

Louisiana State University, Baton Rouge, LA, USA

b

Institut universitaire de France (IUF), France

c

also at University of Bucharest, Physics Department, Bucharest, Romania

d

now at ECAP, Erlangen, Germany

e

School of Physics and Astronomy, University of Leeds, Leeds, United Kingdom

f

now at Agenzia Spaziale Italiana (ASI). Via del Politecnico 00133, Roma, Italy

g

Fermi National Accelerator Laboratory, Fermilab, Batavia, IL, USA

h

now at Graduate School of Science, Osaka Metropolitan University, Osaka, Japan

i

Max-Planck-Institut für Radioastronomie, Bonn, Germany

j

also at Kapteyn Institute, University of Groningen, Groningen, The Netherlands

k

Colorado State University, Fort Collins, CO, USA

l

Pennsylvania State University, University Park, PA, USA

##  Appendix A Parameterization of MC templates

The MC templates of the two-dimensional distributions of ($S$, $X$) normalized
by the total number of simulated showers $N_{\text{MC}}$ and weighted to
correspond to the measured energy spectrum are fitted using the following
function:

$\displaystyle\Phi$
$\displaystyle=\frac{\mathrm{d}N}{N_{\text{MC}}\,\mathrm{d}X\,\mathrm{d}S}$
$\displaystyle=A_{\text{Gauss}}\,A_{\text{Gumbel}}\,f_{\text{Gumbel}}(X)\,f_{\text{Gauss}}(X,S).$
(10)

The $X_{\text{max}}$ part is described by the generalized Gumbel distribution
[52]

$f_{\text{Gumbel}}(X)=\exp[-\lambda(x-\exp x))],$ (11)

where $x=(X-m)/s$ and the ground-signal part is assumed to follow the Gaussian
distribution with the mean value linearly dependent on $X$

$f_{\text{Gauss}}(X,S)=\exp\left[-\frac{y^{2}}{2r^{2}}\right],$ (12)

where $y=S-p\,X-q$. The normalization terms are given by

$A_{\text{Gauss}}=\frac{1}{\sqrt{2\pi}r}\quad\text{and}\quad
A_{\text{Gumbel}}=\frac{\lambda^{\lambda}}{s\,\Gamma(\lambda)},$ (13)

where $\Gamma$ is the gamma function. The six free parameters in each of the
MC template fits are $m$, $s$, $\lambda$ of the generalized Gumbel
distribution, and $p$, $q$, $r$ of the Gaussian part. These fitted parameters
for the three models used in this work are listed in Tables 5, 6 and 7.

In Fig. 11, we show the mean and the standard deviation of the pull
distribution for the description of each bin of the two-dimensional
distribution ($S$, $X$), see examples in Fig. 1, by the fits with a function
Eq. 10. The goodness of the description of the MC templates with these fits is
also tested with the two-dimensional Kolmogorov-Smirnov test [53]
demonstrating very good consistency between the MC templates and MC
distributions of ($S$, $X$).

Figure 11: The mean and the standard deviation of the MC distribution of $\text{PULL}=(N_{\text{fit}}-N_{\text{mc}})/\sqrt{N_{\text{mc}}}$ where $N_{\text{fit}}$ is the value of the fitted function Eq. 10 and $N_{\text{mc}}$ is the value for each bin of the two-dimensional MC distribution ($S$, $X$). The probability of consistency between the parameterized function and MC distribution (red) is tested with the two-dimensional Kolmogorov-Smirnov (K-S) test. $E=10^{18.5}\text{ to }10^{19.0}\,\text{eV}$. Table 5: Parameters of MC templates for air showers generated with model Epos-LHC and initiated by a primary particle $i$, see Eq. 10. $~{}~{}~{}i~{}~{}~{}$ | $~{}~{}~{}~{}~{}~{}~{}\theta~{}~{}~{}~{}~{}~{}~{}$ | $~{}~{}~{}~{}f_{\text{had}}~{}~{}~{}~{}$ |  $~{}~{}m$/(g/cm2) |  $~{}~{}s$/(g/cm2) | $~{}~{}~{}~{}~{}\lambda~{}~{}~{}~{}~{}$ |  $~{}~{}p$/(VEM/(g/cm2)) |  $~{}~{}q$/VEM |  $~{}~{}r$/VEM
---|---|---|---|---|---|---|---|---
p | $(0^{\circ},33^{\circ})$ | $0.64$ | $757$ | $49$ | $1.2$ | $-1.6\times 10^{-2}$ | $3.10$ | $34.3$
p | $(33^{\circ},39^{\circ})$ | $0.67$ | $755$ | $44$ | $1.1$ | $-3\times 10^{-3}$ | $2.85$ | $22.9$
p | $(39^{\circ},45^{\circ})$ | $0.71$ | $755$ | $43$ | $8.9\times 10^{-1}$ | $6.7\times 10^{-3}$ | $2.88$ | $13.8$
p | $(45^{\circ},51^{\circ})$ | $0.77$ | $754$ | $47$ | $9.8\times 10^{-1}$ | $1.2\times 10^{-2}$ | $2.83$ | $6.3$
p | $(51^{\circ},60^{\circ})$ | $0.88$ | $757$ | $53$ | $1.3$ | $1.1\times 10^{-2}$ | $2.68$ | $3.0$
He | $(0^{\circ},33^{\circ})$ | $0.67$ | $741$ | $51$ | $2.1$ | $-2.7\times 10^{-3}$ | $2.78$ | $25.3$
He | $(33^{\circ},39^{\circ})$ | $0.70$ | $739$ | $49$ | $2.0$ | $7.7\times 10^{-3}$ | $2.60$ | $16.1$
He | $(39^{\circ},45^{\circ})$ | $0.75$ | $742$ | $61$ | $2.6$ | $1.7\times 10^{-2}$ | $2.47$ | $7.2$
He | $(45^{\circ},51^{\circ})$ | $0.81$ | $741$ | $57$ | $2.5$ | $1.7\times 10^{-2}$ | $2.55$ | $3.5$
He | $(51^{\circ},60^{\circ})$ | $0.91$ | $739$ | $64$ | $2.7$ | $1.3\times 10^{-2}$ | $2.50$ | $2.6$
O | $(0^{\circ},33^{\circ})$ | $0.71$ | $714$ | $63$ | $5.7$ | $-1.6\times 10^{-3}$ | $2.61$ | $26.2$
O | $(33^{\circ},39^{\circ})$ | $0.74$ | $714$ | $65$ | $5.4$ | $8.8\times 10^{-3}$ | $2.45$ | $16.9$
O | $(39^{\circ},45^{\circ})$ | $0.78$ | $714$ | $66$ | $5.5$ | $1.9\times 10^{-2}$ | $2.41$ | $6.8$
O | $(45^{\circ},51^{\circ})$ | $0.84$ | $711$ | $48$ | $2.9$ | $2.2\times 10^{-2}$ | $2.37$ | $1.8$
O | $(51^{\circ},60^{\circ})$ | $0.93$ | $714$ | $82$ | $7.5$ | $1.7\times 10^{-2}$ | $2.33$ | $1.0$
Fe | $(0^{\circ},33^{\circ})$ | $0.74$ | $688$ | $76$ | $13.2$ | $1.4\times 10^{-2}$ | $2.65$ | $17.4$
Fe | $(33^{\circ},39^{\circ})$ | $0.77$ | $688$ | $100$ | $19.4$ | $1.3\times 10^{-2}$ | $2.62$ | $15.7$
Fe | $(39^{\circ},45^{\circ})$ | $0.81$ | $687$ | $81$ | $13.4$ | $1.4\times 10^{-2}$ | $2.51$ | $11.9$
Fe | $(45^{\circ},51^{\circ})$ | $0.87$ | $688$ | $134$ | $32.7$ | $1.8\times 10^{-2}$ | $2.41$ | $6.2$
Fe | $(51^{\circ},60^{\circ})$ | $0.94$ | $687$ | $176$ | $48.1$ | $1.3\times 10^{-2}$ | $2.34$ | $4.9$
Table 6: Same as in Table 5, but for model QGSJet-II-04. $~{}~{}~{}i~{}~{}~{}$ | $~{}~{}~{}~{}~{}~{}~{}\theta~{}~{}~{}~{}~{}~{}~{}$ | $~{}~{}~{}~{}f_{\text{had}}~{}~{}~{}~{}$ |  $~{}~{}m$/(g/cm2) |  $~{}~{}s$/(g/cm2) | $~{}~{}~{}~{}~{}\lambda~{}~{}~{}~{}~{}$ |  $~{}~{}p$/(VEM/(g/cm2)) |  $~{}~{}q$/VEM |  $~{}~{}r$/VEM
---|---|---|---|---|---|---|---|---
p | $(0^{\circ},33^{\circ})$ | $0.62$ | $741$ | $49$ | $1.1$ | $-1\times 10^{-2}$ | $2.67$ | $27.6$
p | $(33^{\circ},39^{\circ})$ | $0.65$ | $748$ | $70$ | $1.9$ | $1.1\times 10^{-3}$ | $2.53$ | $18.0$
p | $(39^{\circ},45^{\circ})$ | $0.69$ | $743$ | $56$ | $1.3$ | $1.2\times 10^{-2}$ | $2.54$ | $8.0$
p | $(45^{\circ},51^{\circ})$ | $0.76$ | $739$ | $47$ | $9.5\times 10^{-1}$ | $1.6\times 10^{-2}$ | $2.51$ | $2.2$
p | $(51^{\circ},60^{\circ})$ | $0.88$ | $737$ | $50$ | $1.0$ | $1\times 10^{-2}$ | $2.35$ | $2.7$
He | $(0^{\circ},33^{\circ})$ | $0.65$ | $727$ | $68$ | $2.7$ | $3.3\times 10^{-3}$ | $2.49$ | $19.0$
He | $(33^{\circ},39^{\circ})$ | $0.69$ | $729$ | $88$ | $4.7$ | $1.2\times 10^{-2}$ | $2.35$ | $11.0$
He | $(39^{\circ},45^{\circ})$ | $0.74$ | $727$ | $74$ | $3.1$ | $1.8\times 10^{-2}$ | $2.38$ | $4.5$
He | $(45^{\circ},51^{\circ})$ | $0.80$ | $725$ | $60$ | $2.3$ | $2\times 10^{-2}$ | $2.20$ | $4.6\times 10^{-1}$
He | $(51^{\circ},60^{\circ})$ | $0.91$ | $724$ | $74$ | $3.0$ | $1.1\times 10^{-2}$ | $2.21$ | $3.2$
O | $(0^{\circ},33^{\circ})$ | $0.69$ | $700$ | $99$ | $9.5$ | $1.1\times 10^{-2}$ | $2.40$ | $15.1$
O | $(33^{\circ},39^{\circ})$ | $0.73$ | $702$ | $87$ | $6.6$ | $1.5\times 10^{-2}$ | $2.38$ | $10.9$
O | $(39^{\circ},45^{\circ})$ | $0.77$ | $702$ | $91$ | $7.5$ | $2.1\times 10^{-2}$ | $2.27$ | $4.0$
O | $(45^{\circ},51^{\circ})$ | $0.84$ | $704$ | $102$ | $8.0$ | $2.1\times 10^{-2}$ | $2.22$ | $8.4\times 10^{-1}$
O | $(51^{\circ},60^{\circ})$ | $0.93$ | $701$ | $84$ | $5.6$ | $1.5\times 10^{-2}$ | $2.15$ | $1.8$
Fe | $(0^{\circ},33^{\circ})$ | $0.72$ | $675$ | $68$ | $6.7$ | $1.6\times 10^{-2}$ | $2.41$ | $13.9$
Fe | $(33^{\circ},39^{\circ})$ | $0.76$ | $676$ | $77$ | $8.2$ | $2.1\times 10^{-2}$ | $2.37$ | $8.2$
Fe | $(39^{\circ},45^{\circ})$ | $0.81$ | $676$ | $91$ | $11.9$ | $2.2\times 10^{-2}$ | $2.37$ | $4.7$
Fe | $(45^{\circ},51^{\circ})$ | $0.87$ | $675$ | $87$ | $10.7$ | $1.6\times 10^{-2}$ | $2.24$ | $6.4$
Fe | $(51^{\circ},60^{\circ})$ | $0.94$ | $675$ | $162$ | $31.6$ | $1.1\times 10^{-2}$ | $2.19$ | $5.6$
Table 7: Same as in Table 5, but for model Sibyll 2.3d. $~{}~{}~{}i~{}~{}~{}$ | $~{}~{}~{}~{}~{}~{}~{}\theta~{}~{}~{}~{}~{}~{}~{}$ | $~{}~{}~{}~{}f_{\text{had}}~{}~{}~{}~{}$ |  $~{}~{}m$/(g/cm2) |  $~{}~{}s$/(g/cm2) | $~{}~{}~{}~{}~{}\lambda~{}~{}~{}~{}~{}$ |  $~{}~{}p$/(VEM/(g/cm2)) |  $~{}~{}q$/VEM |  $~{}~{}r$/VEM
---|---|---|---|---|---|---|---|---
p | $0^{\circ}-33^{\circ}$ | $0.63$ | $769$ | $70$ | $2.1$ | $-9\times 10^{-3}$ | $2.97$ | $27.7$
p | $33^{\circ}-39^{\circ}$ | $0.65$ | $771$ | $68$ | $1.8$ | $7.8\times 10^{-4}$ | $2.96$ | $19.0$
p | $39^{\circ}-45^{\circ}$ | $0.69$ | $772$ | $71$ | $2.0$ | $1.1\times 10^{-2}$ | $2.84$ | $9.0$
p | $45^{\circ}-51^{\circ}$ | $0.76$ | $766$ | $62$ | $1.5$ | $1.4\times 10^{-2}$ | $2.73$ | $4.4$
p | $51^{\circ}-60^{\circ}$ | $0.87$ | $766$ | $52$ | $1.1$ | $9.2\times 10^{-3}$ | $2.77$ | $4.0$
He | $0^{\circ}-33^{\circ}$ | $0.66$ | $747$ | $56$ | $2.2$ | $-1.3\times 10^{-3}$ | $2.68$ | $23.3$
He | $33^{\circ}-39^{\circ}$ | $0.69$ | $745$ | $50$ | $1.8$ | $1\times 10^{-2}$ | $2.55$ | $13.1$
He | $39^{\circ}-45^{\circ}$ | $0.73$ | $742$ | $54$ | $1.8$ | $1.6\times 10^{-2}$ | $2.51$ | $6.8$
He | $45^{\circ}-51^{\circ}$ | $0.80$ | $744$ | $58$ | $2.1$ | $2\times 10^{-2}$ | $2.41$ | $7.7\times 10^{-1}$
He | $51^{\circ}-60^{\circ}$ | $0.90$ | $743$ | $58$ | $2.1$ | $1.5\times 10^{-2}$ | $2.41$ | $6\times 10^{-1}$
O | $0^{\circ}-33^{\circ}$ | $0.70$ | $715$ | $48$ | $2.3$ | $8.2\times 10^{-3}$ | $2.44$ | $18.1$
O | $33^{\circ}-39^{\circ}$ | $0.73$ | $716$ | $65$ | $4.3$ | $1.6\times 10^{-2}$ | $2.41$ | $10.8$
O | $39^{\circ}-45^{\circ}$ | $0.77$ | $716$ | $62$ | $3.6$ | $2.1\times 10^{-2}$ | $2.36$ | $4.9$
O | $45^{\circ}-51^{\circ}$ | $0.84$ | $717$ | $57$ | $3.0$ | $1.8\times 10^{-2}$ | $2.26$ | $3.8$
O | $51^{\circ}-60^{\circ}$ | $0.93$ | $717$ | $67$ | $3.9$ | $1.5\times 10^{-2}$ | $2.25$ | $1.8$
Fe | $0^{\circ}-33^{\circ}$ | $0.73$ | $688$ | $60$ | $5.6$ | $1.4\times 10^{-2}$ | $2.61$ | $16.2$
Fe | $33^{\circ}-39^{\circ}$ | $0.77$ | $690$ | $62$ | $5.8$ | $2.2\times 10^{-2}$ | $2.40$ | $8.6$
Fe | $39^{\circ}-45^{\circ}$ | $0.81$ | $690$ | $60$ | $5.5$ | $2.3\times 10^{-2}$ | $2.40$ | $5.5$
Fe | $45^{\circ}-51^{\circ}$ | $0.87$ | $690$ | $97$ | $12.1$ | $2\times 10^{-2}$ | $2.35$ | $3.8$
Fe | $51^{\circ}-60^{\circ}$ | $0.94$ | $690$ | $96$ | $11.8$ | $1.9\times 10^{-2}$ | $2.24$ | $8.7\times 10^{-1}$

##  Appendix B MC-MC Tests of the Method

To evaluate the precision and reliability of the method, we performed tests on
the MC simulations with artificially modified predictions on the event-by-
event basis using the shift $\Delta X_{\text{max}}$ and the factors
$R_{\text{had}}(\theta_{\text{min}})$, $R_{\text{had}}(\theta_{\text{max}})$
obtained from the fits to the measured data, see Table 3.

For each model, all possible mixes of (p, He, O, Fe) with relative fractions
changing in 0.1 steps are used to estimate the biases on modification
parameters as a function of the mean logarithmic mass of the primary beam
$\langle{\ln A}\rangle$. In the tests, we use only the composition mixes
containing primary iron nuclei, as there can be additional biases stemming
from the assumption on the presence of different species in the primary beam,
see discussion in Section IV.2. For each composition mix, five different
samples of simulated showers are randomly selected with the same number of
events as in the data ($N=2239$) and following the shape of the measured
energy spectrum. The biases of the method for each composition mix are
calculated as an average over the biases for these five random sets. The
results of the MC-MC tests are summarized in Fig. 12. We find that for
$2<\langle{\ln A}\rangle<3$, the range of the masses inferred from the data
analysis (Section III.2), the maximum overestimation of the fitted
$X_{\text{max}}$ scale with the method is 5 g/cm2, and biases on the fitted
$\langle{S_{\text{had}}}\rangle$ scale are within 3%. These resulting values
of these MC-MC tests are considered as systematic uncertainties of the method
on the modification parameters, see Fig. 16. The systematic uncertainties of
the method on the primary fractions are within $\pm 20\%$, which corresponds
to the maximum bias of the method on the fractions inferred from ${\sim}90$%
of the studied MC samples, see Fig. 13.

The performance of the method [20] is compatible with these estimations of
systematic uncertainties also for a test when one of the models is used for
the analysis of the MC samples prepared with other models. This is in
agreement with the expectations that the main differences in MC templates
between the models are due to the differences in $X_{\text{max}}$ and
$S_{\text{had}}$ scales. Finally, we checked that the number of zenith-angle
bins does not affect the results unless it is too small (2 bins are not enough
to disentangle the hadronic and em parts of the ground signal) or too large
(${>}10$) when the event statistics per zenith-angle bin becomes too low.

Figure 12: The values of modification parameters found by the method (points)
in the MC samples with artificially modified $\Delta X_{\text{max}}$ and
$R_{\text{had}}(\theta)$. The values of artificial modifications, shown with
solid horizontal lines, are applied to each simulated shower individually. The
MC samples contain all possible combinations of primary nuclei with 0.1 steps
in relative fractions, iron nuclei are present in all samples (relative
fraction ${\geq}0.1$), and the mean logarithmic masses of the samples are
marked on horizontal axis. Hadronic interaction models used in the tests:
Epos-LHC (top row, red), QGSJet-II-04 (middle row, blue), Sibyll 2.3d (bottom
row, green). The systematic uncertainties are taken as the maximum biases
within the $\langle{\ln A}\rangle=2$ to 3 range indicated with dashed lines.

Figure 13: Biases of the method on the (p, He, O, Fe) fractions for the MC
samples with artificially modified $\Delta X_{\text{max}}$ and
$R_{\text{had}}(\theta)$ as described in Fig. 12. Additionally, we select the
fitted relative fractions with $0.0\leq f_{\text{p}}^{\text{FIT}}\leq 0.4$,
$0.0\leq f_{\text{He}}^{\text{FIT}}\leq 0.3$, $0.2\leq
f_{\text{O}}^{\text{FIT}}\leq 0.6$ and $0.0\leq f_{\text{Fe}}^{\text{FIT}}\leq
0.5$ which correspond to the ranges obtained by the data fits. The range
considered as a contribution to the systematic uncertainty is indicated by
vertical lines. Hadronic interaction models used in the tests: Epos-LHC (top
left, red), QGSJet-II-04 (top-right, blue), Sibyll 2.3d (bottom, green).

##  Appendix C Parameterization of Attenuation of Ground Signals

The electromagnetic and hadronic signals at 1000 m were corrected for the
energy evolution as
$S_{\text{em}}^{\text{ref}}=S_{\text{em}}\,(E^{\text{ref}}/E_{\text{FD}})$ and
$S_{\text{had}}^{\text{ref}}=S_{\text{had}}\,(E^{\text{ref}}/E_{\text{FD}})^{1/\beta}$.
The dependence of these average signals on the distance of $X_{\text{max}}$ to
the ground in atmospheric depth units,
$t=X_{\text{atm}}(\theta)-X_{\text{max}}$, see Fig. 2, was parameterized with
the Gaisser-Hillas function [54] allowing its vertical offset,

$\langle{S_{\alpha}^{\text{ref}}}\rangle(t)=S_{\alpha}^{0}\,\left(\frac{t-t_{\alpha}}{u_{\alpha}-t_{\alpha}}\right)^{Z(t_{\alpha})}\,\exp(Z(t))+w_{\alpha},$
(14)

where $\alpha=\text{had}$ or em, and normalization scaling
$Z(x)=(u_{\alpha}-x)/v_{\alpha}$. $t_{\alpha}$ is the value of
$X_{\text{atm}}(\theta)-X_{\text{max}}$ where the function reaches its
maximum, $u_{\alpha}$ and $v_{\alpha}$ are parameters without a
straightforward physics interpretation, and $S^{0}_{\alpha}$, $w_{\alpha}$ are
the rescale and offset parameters, respectively. The fitted parameters are
listed in Table 8 and Table 9 for em and muon signal, respectively.

The factor $g_{\text{em},k}$ reflecting the average change of em signal due to
the change of $X_{\text{max}}$ scale for a mix of primary species $i$ with
relative fractions $f_{i}$ in zenith-angle bin $k$ is calculated as

$g_{\text{em},k}=\sum_{i}f_{i}\,\frac{S_{\text{em}}(\langle{t}\rangle_{k}-\Delta
X_{\text{max}})}{S_{\text{em}}(\langle{t}\rangle_{k})},$ (15)

where $\langle{t}\rangle_{k}$ is the average measured
$X_{\text{atm}}(\theta)-X_{\text{max}}$ in a zenith-angle bin $k$: $\sim
241.2$ g/cm2, $\sim 335.5$ g/cm2, $\sim 429.4$ g/cm2, $\sim 556.1$ g/cm2,
$\sim 777.6$ g/cm2, for respective increasing average values of the zenith
angle. The factor $g_{\text{had}}$ for the hadronic signal is obtained as

$g_{\text{had},k}=\sum_{i}f_{i}\,\frac{S_{\text{had}}(\langle{t}\rangle_{k}-\Delta
X_{\text{max}})}{S_{\text{had}}(\langle{t}\rangle_{k})}\,\frac{R_{\text{had}}(\langle{t}\rangle_{k}-\Delta
X_{\text{max}})}{R_{\text{had}}(\langle{t}\rangle_{k})},$ (16)

where

$\displaystyle R_{\text{had}}(\theta)$ $\displaystyle\equiv
R_{\text{had}}(\langle{t}\rangle_{k})=R_{\text{had}}(\theta_{\text{min}})+{}$
(17)
$\displaystyle+(R_{\text{had}}(\theta_{\text{max}})-R_{\text{had}}(\theta_{\text{min}}))\,\frac{\langle{t}\rangle_{k}-\langle{t}\rangle_{\text{min}}}{\langle{t}\rangle_{\text{max}}-\langle{t}\rangle_{\text{min}}},$
(18)

for $\langle{t}\rangle_{\text{min}}$ and $\langle{t}\rangle_{\text{max}}$
corresponding to minimum and maximum zenith-angle bins, respectively.

Table 8: Fitted parameters of functional dependence of the average em signal part of $S(1000)$ vs. distance of $X_{\text{max}}$ to the ground in atmospheric depth units for air showers initiated by a primary particle $i$, see Eq. 14. | $~{}~{}~{}i~{}~{}~{}$ |  $~{}~{}S_{\text{em}}^{\text{0}}$/VEM |  $~{}~{}t_{\text{em}}$/(g/cm2) |  $~{}~{}u_{\text{em}}$/(g/cm2) |  $~{}~{}v_{\text{em}}$/(g/cm2) |  $~{}~{}w_{\text{em}}$/VEM
---|---|---|---|---|---|---
Epos-LHC | p | $7.58$ | $-1{,}638$ | $172$ | $42$ | $1\times 10^{-1}$
Epos-LHC | He | $7.76$ | $-580$ | $179$ | $81$ | $6\times 10^{-2}$
Epos-LHC | O | $7.83$ | $-1{,}412$ | $161$ | $48$ | $9\times 10^{-2}$
Epos-LHC | Fe | $7.64$ | $-407$ | $204$ | $86$ | $7\times 10^{-2}$
QGSJet-II-04 | p | $7.56$ | $-2{,}006$ | $175$ | $35$ | $1\times 10^{-1}$
QGSJet-II-04 | He | $7.71$ | $-495$ | $183$ | $87$ | $3\times 10^{-2}$
QGSJet-II-04 | O | $7.77$ | $-763$ | $167$ | $72$ | $7\times 10^{-2}$
QGSJet-II-04 | Fe | $7.52$ | $-278$ | $206$ | $103$ | $2\times 10^{-3}$
Sibyll 2.3d | p | $7.57$ | $-1{,}152$ | $176$ | $53$ | $1\times 10^{-1}$
Sibyll 2.3d | He | $7.62$ | $-490$ | $185$ | $88$ | $3\times 10^{-2}$
Sibyll 2.3d | O | $7.54$ | $-665$ | $187$ | $73$ | $6\times 10^{-2}$
Sibyll 2.3d | Fe | $7.74$ | $-1{,}892$ | $152$ | $39$ | $9\times 10^{-2}$
Table 9: Same as in Table 8, but for the average hadronic-signal part of $S(1000)$. | $~{}~{}~{}i~{}~{}~{}$ |  $~{}~{}S_{\text{had}}^{\text{0}}$/VEM |  $~{}~{}~{}~{}t_{\text{had}}$/(g/cm2) |  $~{}~{}u_{\text{had}}$/(g/cm2) |  $~{}~{}v_{\text{had}}$/(g/cm2) |  $~{}~{}w_{\text{had}}$/VEM
---|---|---|---|---|---|---
Epos-LHC | p | $5.6$ | $-9{,}999$ | $312$ | $9$ | $8.3$
Epos-LHC | He | $7.2$ | $-9{,}419$ | $232$ | $16$ | $8.4$
Epos-LHC | O | $9.0$ | $-1{,}298$ | $228$ | $98$ | $8.7$
Epos-LHC | Fe | $10.7$ | $-8{,}304$ | $171$ | $24$ | $9.6$
QGSJet-II-04 | p | $5.1$ | $-1{,}880$ | $314$ | $47$ | $7.2$
QGSJet-II-04 | He | $7.2$ | $-300$ | $255$ | $222$ | $6.8$
QGSJet-II-04 | O | $7.2$ | $-4{,}934$ | $205$ | $34$ | $8.7$
QGSJet-II-04 | Fe | $9.6$ | $-1{,}614$ | $166$ | $115$ | $8.7$
Sibyll 2.3d | p | $4.8$ | $-8{,}519$ | $283$ | $11$ | $8.1$
Sibyll 2.3d | He | $6.7$ | $-3{,}557$ | $225$ | $44$ | $8.0$
Sibyll 2.3d | O | $8.4$ | $-1{,}521$ | $207$ | $100$ | $8.4$
Sibyll 2.3d | Fe | $8.9$ | $-5{,}078$ | $194$ | $35$ | $10.1$

##  Appendix D Data description using fits with less freedom in MC templates

The description of Auger data using fits without any modification to MC
templates and with zenith-independent $R_{\text{had}}$ are shown in Figs. 14
and 15, respectively.

Figure 14: Same as in Fig. 5, but for the data fits performed without any
modification of the predictions of models.

Figure 15: Same as in Fig. 5, but for the data fits performed with only
zenith-angle independent modification $R_{\text{had}}$ of the predictions of
models (denoted by ∗ as templates were modified from the original
predictions).

##  Appendix E Systematic Uncertainties

The individual contributions to the total systematic uncertainties are plotted
in Fig. 16.

Figure 16: Contributions of individual experimental systematic uncertainties
on energy, $X_{\text{max}}$, $S(1000)$ and of the method biases inferred from
the MC-MC tests (see Appendix C) to the total systematic uncertainties on
fractions of primary nuclei (top left), $\Delta X_{\text{max}}$ (top right),
$R_{\text{had}}(\theta)$ (bottom) and zenith-angle difference of
$R_{\text{had}}$ (top right). Coloured bands are the total systematic
uncertainties obtained by summing individual contributions in quadrature with
the best-fit results indicated by dashed horizontal lines.

##  Appendix F Scan for the linear combinations of experimental systematic
uncertainties most favorable for the models

Figure 17: Number of scans in linear combinations of experimental systematic
uncertainties for fits using modifications $\Delta X_{\text{max}}$ and
$R_{\text{had}}(\theta)$ to estimate the closest approach to the point
($R_{\text{had}}(\theta)=1$, $\Delta X_{\text{max}}=0$ g/cm2) from fits with
non-modified MC templates in the difference of log-likelihood
$\Delta\ln\mathcal{L}$. The value estimated using Wilks’ theorem in the
Likelihood-ratio test for the nested model at the level of 5$\sigma$
($\Delta\ln\mathcal{L}\approx 16.62$) is indicated by the arrow.

In Fig. 17, the histogram of difference in log-likelihood expressions
($\Delta\ln\mathcal{L}$) for fits using non-modified and $\Delta
X_{\text{max}}$, $R_{\text{had}}(\theta)$ modifications is shown for dense
scans in linear combinations of experimental systematic uncertainties. The
ranges of these uncertainties were selected in a way to estimate the closest
approach to the non-modified value ($R_{\text{had}}(\theta)=1$, $\Delta
X_{\text{max}}=0$ g/cm2) even for cases when the uncertainties were out of the
range of uncertainties quoted by Auger. It was not possible to find such a
linear combination of experimental systematic uncertainties that would
decrease the significance of improvement in data description below $5\sigma$.