# New physics above 50 TeV: probing its phenomenology through UHECR air-shower
simulations

Stylianos Romanopoulos<EMAIL_ADDRESS>Vasiliki Pavliou
<EMAIL_ADDRESS>Department of Physics & Institute of Theoretical and
Computational Physics, University of Crete, GR-70013, Heraklion, Greece
Institute of Astrophysics, Foundation for Research and Technology-Hellas,
Vasilika Vouton, GR-70013 Heraklion, Greece  Theodore Tomaras  deceased
Department of Physics & Institute of Theoretical and Computational Physics,
University of Crete, GR-70013, Heraklion, Greece

###### Abstract

Ground based observations appear to indicate that Ultra High Energy Cosmic
Rays (UHECR) of the highest energies ($>10^{18.7}\rm eV$) consist of heavy
particles – shower depth and muon production data both pointing towards this
conclusion. On the other hand, cosmic-ray arrival directions at energies
$>10^{18.9}{\rm\,eV}$ exhibit a dipole anisotropy, which disfavors heavy
composition, since higher-Z nuclei are strongly deflected by the Galactic
magnetic field, suppressing anisotropy. This is the composition problem of
UHECR. One solution could be the existence of yet-unknown effects in proton
interactions at center-of-mass (CM) energies $\gtrsim 50$ TeV, which would
alter the interaction cross section and the multiplicity of interaction
products, mimicking heavy primaries. We study the impact of such changes on
cosmic-ray observables using simulations of Extensive Air-Shower (EAS), in
order to place constrains on the phenomenology of any new effects for high
energy proton interactions that could be probed by $\sqrt{s}>50$ TeV
collisions. We simulate showers of primaries with energies in the range
$10^{17}-10^{20}{\rm eV}$ using the CORSIKA code, modified to implement a
possible increase in cross-section and multiplicity in hadronic collisions
exceeding a CM energy threshold of $50$ TeV. We study the composition-
sensitive shower observables (shower depth, muons) as a function of cross-
section, multiplicity, and primary energy. We find that in order to match the
Auger shower depth measurements by means of new hadronic collision effects
alone (if extragalactic UHECR are all protons even at the highest energies),
the cross-section of proton-air interactions has to be $\sim$ 800 mb at $140$
TeV CM energy, accompanied by an increase of a factor of 2-3 in secondary
particles. We also study the muon production of the showers in the same
scenario. Although the muon production does increase, this increase is not
enough to resolve the muon problem of UHECRs if the distribution of
secondaries among different species remains unchanged with respect to the
Standard Model prediction.

††preprint: APS/123-QED

##  I Introduction

Cosmic Rays (CR) are the most energetic particles in the universe. Over 100
years have passed since they were first discovered by Hess [1], yet still
today their composition, origin and acceleration mechanism are subjects of
debate. This lingering uncertainty stems not only from observational
limitations but also from particle-physics uncertainties: the first collisions
of ultra-high-energy cosmic rays (UHECRs $E\gtrsim 10^{18}$ eV) with the
Earth’s atmosphere occur at center-of-mass (CM) energies exceeding $40$ TeV
and reaching $300$ TeV; for comparison, current lab tests of hadronic physics
(in the Large Hadron Collider, LHC) only reach CM energies of 14 TeV.

Observationally, the greatest challenges in studying UHECRs are their low flux
and their inability to penetrate the atmosphere: detection has to be indirect,
tracing either the development of the extensive particle air shower (EAS)
caused by a CR’s collision with the atmosphere, or the EAS products reaching
the ground, or both; and dedicated UHECR observatories need collective areas
of thousands of km2.

Despite significant progress in recent decades in improving statistics of
UHECR detections thanks to very lagre facilities such as the Pierre Auger
Observatory in Argentina [2] and the Telescope Array in the United States [3],
important questions regarding UHECR astrophysics remain open. In particular,
the distribution of UHECR arrival directions on the sky, their origin, and
their composition constitute three persistent, interconnected puzzles.

Because cosmic rays are charged, it is not possible to directly associate
their arrival directions with the sources that accelerate them. However, there
are astrophysical arguments about possible classes and cosmic locations of
cosmic-ray sources, with important implications for the resulting cosmic-ray
properties.

At lower energies ($<10^{12}$ eV), it is fairly certain that cosmic rays
originate in Galactic sources. This is evidenced by differences in the cosmic
ray fluxes estimated for different galaxies through observations of gamma rays
originating in the decay of neutral pions produced through collisions of
cosmic rays with interstellar gas (e.g., [4, 5, 6]). At very high energies
($>8\times 10^{18}$ eV) we are equally confident that cosmic rays are
extragalactic, since the anisotropic distribution of arrival directions that
starts emerging at these energies is not correlated with the Galactic plane or
the Galactic center [7].

The exact energy at which the transition from Galactic to extragalactic cosmic
rays occurs is still under debate (see, e.g., [8] for a recent review),
however most recent works assume that it happens somewhere between $10^{17}$
and $10^{18.5}$ eV. Hints for this are seen both in the spectrum and in the
composition of cosmic rays at this energy range. The steepening of the
spectrum seen at the ”second knee” and the transition, at the same energies,
to a heavier composition [9] point towards a population originating in
magnetically-confining accelerators reaching its maximum possible energy (as
per the Hillas criterion, [10, 11]).

At energies between $\sim 10^{17}$ eV and $\sim 10^{18.5}$ eV, composition-
sensitive observables from KASCADE-Grande, Auger, and Telescope Array indicate
a transition back to lighter composition [12, 13, 14]. At $\sim
10^{18.7}{\rm\,eV}$ the spectrum also becomes harder, a feature known as the
ankle [15, 14]. There are two interpretations of the ankle. In the first, the
ankle marks the transition to UHECR of extragalactic origin (e.g., [16, 17,
18]). In the second, this transition is assumed to have already happened at
somewhat lower energies. As a result, the composition is already light at the
ankle, and the ankle spectral feature is actually a ”dip” caused by electron-
positron losses (e.g., [19, 20, 21]).

The debate of ankle-versus-dip is not, however, the only controversy at these
energies. A second one is the so-called composition problem of UHECR. This is
summarized as follows: At energies of $~{}10^{18.7}$ eV, composition-sensitive
variables, taken at face-value, indicate a transition back to heavier
composition [22]. However, there are certain astrophysical indicators against
a heavy composition: (a) The spectrum at this energy is transitioning to a
shallower slope, not a steeper one - i.e, there is no coincident spectral
indication that the UHECR accelerators are reaching their maximum energy [23]
(b) Anisotropies start to emerge at these energies [7, 24, 25, 26, 27]. This
might not be so severe a problem if the Galactic magnetic field has the
overall low strengths indicated, e.g., by [28]. However, recent studies of the
Galactic magnetic field have shown that it is approximately an order of
magnitude stronger than previously thought [29] in a small region near the
reported hotspot from TA [26, 27]. If indeed the average Galactic magnetic
field is proven to be just a few times stronger than the existing models,
combined with the dipole anisotropy at high energies, we can conclude that
UHECR are light nuclei. The reason is that heavy nuclei are strongly deflected
from Galactic magnetic fields and would spread over all the sky, eliminating
all evidence of anisotropy. (c) Heavier nuclei photodissociate fast during
propagation (e.g., [30, 31, 18]) – with the exception of iron – so the
composition becomes lighter during propagation, unless it starts out as pure
iron. However, iron is far from a best-fit to Auger composition-sensitive
observables. Instead, observations can be better fit by a mix of intermediate-
mass nuclei, requiring an astrophysically contrived composition of the
accelerated particles at the source (e.g., [32, 33, 34, 35, 36]). In contrast,
models that are more natural astrophysically are not in as good agreement with
composition-sensitive observables [37, 38].

In addition to these astrophysical considerations, there are particle-physics
considerations that add to the composition problem. Hybrid detectors such as
Auger and Telescope Array measure composition indirectly, in two ways: from
the atmospheric slant depth at which the shower reaches a maximum, $X_{\rm
max}$ (measured through fluorescent detectors); and from the number of muons
reaching the ground (measured by surface array tanks). At a fixed primary CR
energy, heavier nuclei will typically give a lower and less variable $X_{\rm
max}$; and they will produce more muons. The observations of $X_{\rm max}$ and
muon numbers are then compared against the predictions from air-shower
simulations. However, the best-fit compositions from muons and $X_{\rm max}$
do not match [39]: too many muons are produced on the ground compared to what
would be expected from the best-fit composition obtained from $X_{\rm max}$
alone.

The latter two problems indicate that the air-shower simulations (or rather,
the hadronic collision simulation models on which these are based) may not be
capturing correctly the development of showers. This is not altogether
unexpected, since the first collision of a $10^{17}{\rm\,eV}$ cosmic ray with
a stationary atmospheric proton is already super-LHC: we are simulating
collisions of these primaries with the atmosphere based on theoretical
extrapolations of hadronic behavior to higher energies. This has led several
authors to hypothesize that the problem may lie in the hadronic collision
models themselves. The solution that has been proposed in this context is
that, above a threshold energy $E_{\text{th}}$, the proton-Air interaction
changes due to new new physics beyond the SM. This scenario is widely
recognized both by the Auger Collaboration [40, 41, 42] and other authors [43,
44, 45, 46]. In this scenario, the composition of the primary remains light.
In [46] (hereafter PT19) we calculated analytically phenomenological
constraints on any new effect that would alter hadronic interactions in such a
way as to mimic a heavy composition at the highest UHECR energies. We showed
that if the multiplicity of first-collision products increased over the SM
predictions at a certain rate, the growth of the average $X_{\rm max}$ with
energy can be fully explained while keeping the composition light even at the
highest energies. We also showed that a simultaneous increase in the proton-
air cross-section over the SM prediction would improve agreement of
$\sigma_{X_{\rm max}}$ (the shower-to-shower variation in $X_{\rm max}$) with
the data, although we did not calculate the optimal behavior of the cross-
section to best match Auger observations.

In this paper, we extend the analytic formulation of PT19 using EAS
simulations. For that purpose, we use a widely known program:
CORSIKA111CORSIKA version 7.7402 [47]. CORSIKA uses extrapolations of SM at
post-LHC energies to model hadronic interactions. In this work we used two
such models: EPOS LHC [48] and QGSJETII-04 [49]. For low energy interactions
we used FLUKA 222FLUKA version 2020.0.3 [50] which is a Monte Carlo code used
extensively at CERN. Our main goals are: (a) To test whether an increase in
the multiplicity of first-collision products can indeed yield the changes in
the $X_{\rm max}$ distribution predicted in PT09 – as contrasted with, e.g.
the implementation of such a change by [51], who found that the variance of
$X_{\rm max}$ remains practically unchanged under a change in product
multiplicity. (b) To calculate the optimal change in cross-section that best
matches Auger $X_{\rm max}$ data. In contrast, in PT19 we did not fully
explore the parameter space, but we only argued that an increase in cross
section changes the $X_{\rm max}$ variance in the direction of better
agreement with Auger data. (c) To evaluate the impact of this scenario on the
muon problem, which was not addressed in PT19.

This paper is organized as follows. In Section II we present the mathematical
formulation of our implementation of new physics effects in the first
collision of $E>10^{18}{\rm\,eV}$ CR with the atmosphere. We discuss the
results of our simulations in Section III, and we summarize and discuss our
conclusions in Section IV.

##  II Mathematical Formulation

The ”new physics above 50 TeV” scenario that we explore here using simulations
of EAS implements two phenomenological changes in the first collision of the
incoming primary cosmic ray with the atmosphere (as in PT2019): (1) an
increase in the multiplicity of the first collision products; and (2) an
increase in the cross section of the first interaction with the atmosphere.
Our approach is phenomenological and not tied to any specific new physics
model. However, several candidate particles and new physics mechanisms exist
that might lead to such behavior (see, for example, [52, 53, 54, 43]).

In §II.1 we discuss how we describe quantitatively these phenomenological
changes. We discuss how these changes impact the slant depth of the showers in
§II.2, and the muons reaching the ground in §II.3. In §II.4 we describe how we
implement the changes of §II.1 in air shower simulations using CORSIKA.

In this paper, we assume that all extragalactic cosmic rays reaching the Earth
are protons. However, not all $E>10^{18}{\rm eV}$ cosmic rays detected on
Earth are extragalactic, and, more importantly, the high-energy-end of the
Galactic cosmic ray spectrum has a heavy composition, additionally affecting
the slant depth and muon content. In order therefore to compare our results
with observations, we additionally need a model for the way the Galactic
cosmic ray flux cuts off with energy. We describe this model in §II.5.

###  II.1 Parametrization of changes in cross section and multiplicity

The cross-section of protons with nitrogen has a logarithmic behaviour at high
energies [55]. For that reason, the cross section is usually parameterized as

$\sigma_{\text{p-Air}}=\sigma_{0}+\beta\,\log\varepsilon$ (1)

where $\sigma_{0}$ and $\beta$ are constants. Here we normalize the energy
scale as $\varepsilon=E/E_{\text{th}}$, where $E_{\text{th}}$ is the threshold
energy above which new physics sets in. Based on the arguments in PT2019, we
will take $E_{\text{th}}=10^{18}$ eV, corresponding, for a collision of an
primary proton with a stationary atmospheric proton, to a CM collision energy
of $\sim 50$ TeV. We note that this threshold is neither coincident with the
location of the break seen in Auger $X_{\rm max}$ data, nor fine-tuned (other
choices that satisfy both that the threshold is ultra-LHC and that it lies
below the Auger data break give similarly good results). We can calculate
$\sigma_{0}$ and $\beta$ from hadronic interaction models (here we use EPOS
LHC and QGSJETII-04) and the standard-model extrapolations employed therein.
The results are given in the first 2 lines of Table (1).

| EPOS LHC | QGSJETII-04
---|---|---
$\sigma_{0}$ (mb) | 527.61 | 499.35
$\beta$ (mb) | 49.95 | 37.49
$X_{0}$ (gr$/$cm2) | 707.38 | 692.03
$\alpha$ (gr$/$cm2) | 63.74 | 60.37
Table 1: Parameters for cross section and shower maximum. These parameters
were calculated by performing linear fit in the SM data from EAS simulation
with proton as a projectile.

If new phenomena take place at $\varepsilon>1$, the cross section of the first
interaction may change. Here we assume that such a change will only affect the
value of the coefficient $\beta$ of the energy-dependent term, and that the
cross section will be continuous at $\varepsilon=1$:

$\sigma_{\text{p-Air,new}}=\sigma_{0}+\beta^{\prime}\,\log\varepsilon$ (2)

We will parameterize this change in terms of the fractional change $\delta$ in
the coefficient $\beta$ relative to its standard-model value. Defining then

$\delta=(\beta^{\prime}-\beta)/\beta\,,$ (3)

we obtain

$\sigma_{\text{p-Air,new}}=\sigma_{\text{p-Air}}+\delta\,\beta\,\log\varepsilon\,.$
(4)

In other words, for $\varepsilon>1$, the cross section will deviate
logarithmically from its standard-model–predicted value at the particular
energy, with a coefficient $\delta\beta$. Clearly, $\delta=0$ corresponds to
no change to the cross-section at any energy over the standard-model
prediction. We note however that the uncertainty in the standard-model
predictions for the proton-air cross-section to super-LHC energies is high
(see, e.g., Fig. 2 of [51]), and values of $\delta$ as high as 3.5 could still
be consistent with the standard model within uncertainties.

We also postulate an increase of the number of secondary particles produced
after the first collision of the primary with the atmosphere. We limit the
effect to the first collision since, for energies of interest, secondary
particles will, with very high probability, have energies such that their
collisions with air will occur at CM energies below the threshold for new
physics. We parameterize this increase in first-collision product multiplicity
by

$n(\varepsilon)=N(\varepsilon)/N_{\text{SM}}(\varepsilon)\,,$ (5)

where $N$ is the number of secondaries produced under new physics, and
$N_{\text{SM}}$ is the standard-model prediction for the number of
secondaries. It is possible that new physics may also change the charged-
particle ratio of the products; we do not however implement such a change
here. We further discuss this issue in §III.

###  II.2 Shower maximum

The slant depth of the shower maximum is the air column density traversed by
the EAS front – measured from the top of the atmosphere – for which the shower
front achieves its maximum atmospheric ionization rate:

$X_{\text{max}}=\int\limits_{\infty}^{x_{\text{max}}}\rho(l)\,dl\,.$ (6)

In Eq. (6), $\rho$ is the density of the atmosphere, $l$ is a length measured
along the path the shower particles traverse in it, and $x_{\text{max}}$ is
the height of the shower maximum. $X_{\text{max}}$ quantifies the total
atmospheric column the shower has already encountered at its maximum,
independently of the inclination of the incoming CR. $X_{\text{max}}$ can be
written as the sum of two terms. The first one is the column density after
which the first interaction of the cosmic ray primary takes place,
$X_{\text{int}}$; and the second one is the column density between first
interaction and shower maximum, $X_{\text{long}}$, corresponding to the
”longitudinal” development of the shower:

$X_{\text{max}}=X_{\text{int}}+X_{\text{long}}\,.$ (7)

####  II.2.1 The first interaction: $X_{\rm int}$

The probability that a CR has not interacted with the atmosphere in the
vicinity of height $x$ is

$\exp\Bigg{(}-\frac{\sigma_{\text{CR-
Air}}}{m}\int\limits_{\infty}^{x}\rho(l)\;dl\;\Bigg{)},$

where $m$ is the average mass of the particles in the atmosphere (mainly
nitrogen) and $\sigma_{\text{CR-Air}}$ is the cosmic-ray–air cross section for
the energy of the primary. The average value of the depth of the first
interaction for a given primary energy thus is

$\langle X_{\text{int}}\rangle=\frac{m}{\sigma_{\text{CR-
Air}}(\varepsilon)}\,.$ (8)

Since $X_{\text{int}}$ follows Poisson statistics, its variance will be

$\text{Var}(X_{\text{int}})=\langle
X_{\text{int}}\rangle^{2}=\frac{m^{2}}{\sigma^{2}_{\text{CR-
Air}}(\varepsilon)}\,.$ (9)

If new physics sets in for $\varepsilon>1$, the distribution of
$X_{\text{int}}$ for primaries of a given energy will be affected through the
change in $\sigma_{\text{CR-Air}}$, which (assuming the composition remains
light) will take the value of $\sigma_{\text{p,Air,new}}$ given by Eq. (4):

$\langle
X_{\text{int,new}}\rangle=\frac{m}{\sigma_{\text{p-Air}}(\varepsilon)+\delta\,\beta\,\log\varepsilon}\,,$
(10)

and

$\text{Var}(X_{\text{int,new}})=\frac{m^{2}}{[\sigma_{\text{p-Air}}(\varepsilon)+\delta\,\beta\,\log\varepsilon]^{2}}\,.$
(11)

####  II.2.2 The longitudinal development: $X_{\rm long}$

Following the CR interaction with the atmosphere, secondary and subsequent
generation of particles are produced. As the shower of particles evolve though
the atmosphere, the energy per particle decreases, partly because of
ionization losses and partly because of new particle production. This process
continues until the energy of the shower particles drops below the energy
threshold of new particle production, at which point the shower continues to
evolve by ionization losses alone. This evolution of the shower is observable
through the fluorescence of ionized atoms of the atmosphere, which can be
detected with ground telescopes. The intensity of this fluorescent light
encodes the energy loss rate of the shower front. $X_{\rm long}$ measures the
column density traversed by the shower front after the first interaction until
the energy loss rate reaches its maximum value.

The dependence of the longitudinal column, $X_{\text{long}}$, with energy can
be derived from the simple model of Heitler [56, 57], where the initial
particle produces two daughter particles which split the primary’s energy, and
the process continues until the energy of each particle reaches a critical
energy below which the process cannot continue.

This results to a logarithmic increase of $X_{\text{long}}$ with energy. A
more realistic calculation of $X_{\text{long}}$ is more complex because
additional phenomena take place (e.g. bremsstrahlung, pair production, pion
production, hadronization), and there are significant shower-to-shower
fluctuations. For that purpose numerical simulations are used to follow the
development of EAS (e.g. CORSIKA). The final behavior of $\langle
X_{long}\rangle$ with energy, however, is still logarithmic:

$\langle X_{\text{long}}\rangle=X_{0}+\alpha\,\log\varepsilon$ (12)

where $X_{0}$ and $\alpha$ are constants, while $\sigma_{X_{\text{long}}}$
remains approximately constant with energy.

In the last 2 lines of Table (1) we show the best-fit parameters $X_{0}$ and
$\alpha$ for protons, derived from CORSIKA EAS simulations using two different
models (EPOS LHC and QGSJETII-4) of hadronic interactions. Both models employ
standard-model extrapolations for collisions at super-LHC energies. The
parameters do depend slightly on the hadronic interactions model, however both
models produce qualitatively similar results.

If now new physics sets in for $\varepsilon>1$, the distribution of
$X_{\text{long}}$ will change. To quantify this change, we model empirically
the shower as $n(\epsilon)$ ”component showers”, of energy
$\epsilon/n(\epsilon)$ on average, developing independently. Note that this
approach is conceptually and qualitatively different from that of [51] who
used a multiplicative factor to increase the number of products in each
collision - the most prominent difference being that the presence of
independently developing ”component showers” decreases the shower-to-shower
fluctuations of $X_{\text{long}}$ (i.e. $\sigma_{X_{\text{long}}}$), while a
multiplicative increase of products of identical distribution as the original
shower leaves $\sigma_{X_{\text{long}}}$ unchanged.

Under this change, $\langle X_{\rm long,new}\rangle$ becomes [from Eq. (12)]

$\langle
X_{\text{long,new}}\rangle=X_{0}+\alpha\,\log\frac{\varepsilon}{n(\epsilon)}\,.$
(13)

To produce an analytical estimate the variance of $X_{\text{long}}$ we take,
as in PT19, the average of individual ”component-shower” longitudinal depth,
$\cfrac{1}{n}\sum\limits_{i}X_{\rm long,i}$ to be a reasonable estimation of
$X_{\rm long}$. Then $X_{\rm long}$ is the “sample mean” of $n$ “draws” from
the underlying distribution of $X_{\rm long,i}$ and the distribution of these
“sample means” has a variance that is given by the “error in the mean”
formula,

$\text{Var}(X_{\text{long,new}})=\frac{\text{Var}(X_{\text{long,i}})}{n(\varepsilon)}\,.$
(14)

Here $\text{Var}(X_{\text{long,i}})$ is the variance of $X_{\text{long,i}}$,
and it can be assumed to follow the SM predictions, since each subshower of
the sample will have $\epsilon<1$.

####  II.2.3 New physics mimics a transition to heavy composition

The average value and the variance of $X_{\text{max}}$ are sensitive to the
first interaction - both its cross-section and the multiplicity of its
products. A higher first-interaction cross-section will result to lower
$\langle X_{\text{int}}\rangle$ and $\text{Var}(X_{\text{int}})$, and
therefore lower $\langle X_{\text{max}}\rangle$. A large number of first
collision products will distribute the energy of the primary more widely,
resulting to lower $\langle X_{\text{long}}\rangle$ and therefore to lower
$\langle X_{\text{max}}\rangle$. A larger number of first-collision products
will also result to reduced shower-to-shower fluctuations in $X_{\text{long}}$
and in turn to a lower $\text{Var}(X_{\text{long}})$ and a lower
$\text{Var}(X_{\text{max}})$.

A heavy primary composition will drive booth $\langle X_{\text{max}}\rangle$
and $\text{Var}(X_{\text{max}})$ to lower values for a given energy, through
both an increased cross-section of the primary-air collision, and an increased
first-collision product multiplicity. The new-physics scenario we discuss here
will also move the $X_{\text{max}}$ distribution in the same direction through
similar changes in cross-section and multiplicity, and can thus mimic a
transition to a heavier composition.

####  II.2.4 Constraining $n(\varepsilon)$ and $\delta$

The empirical model we have presented here features a free parameter $\delta$
and a free function $n(\varepsilon)$. However, if we assume that the change of
slope in $\langle X_{\text{max}}(\epsilon)\rangle$ observed by Auger at the
highest energies is, in its entirety, due to new physics, so that all
primaries at these energies are protons, we can completely determine
$n(\epsilon)$ as a function of $\delta$, leaving only a single free parameter
in our model (which can also be optimized by comparison to other moments of
the $X_{\text{max}}$ distribution, as we will see in §III).

The Auger Collaboration reports, for energies $E\geq 10^{18.3}$ eV

$\langle X_{\text{max, Auger}}\rangle=X_{0,\text{
Auger}}+\alpha_{\text{Auger}}\,\log\varepsilon$ (15)

with $X_{0,\text{Auger}}=749.74\;\text{gr}/\text{cm}^{2}$ and
$\alpha_{\text{Auger}}=23.98\;\text{gr}/\text{cm}^{2}$. Simulated data of pure
protons air showers indicate a difference from Auger observations above this
energy. By equating the behaviour of Auger data to the results of the
postulated new physics effects at $\epsilon>1$ [using Eqs. (7), (10), and
(13)] we obtain

$\displaystyle\langle
X_{\text{max,Auger}}\rangle(\varepsilon)=\cfrac{m}{\sigma_{\text{p-Air}}(\varepsilon)+\delta\,\beta\,\log\varepsilon}$
$\displaystyle+X_{0}+\alpha\,\log\cfrac{\varepsilon}{n(\varepsilon)}$ (16)

and solving for the multiplicity we obtain

$\displaystyle\log n(\varepsilon)=\frac{X_{0}-X_{0,\text{Auger}}}{\alpha}$
$\displaystyle+\cfrac{\alpha-\alpha_{\text{Auger}}}{\alpha}\log\varepsilon$
$\displaystyle+\cfrac{1}{\alpha}\cfrac{m}{\sigma_{\text{p-Air}}(\varepsilon)+\delta\,\beta\,\log\varepsilon}$
(17)

with $\sigma_{\text{p-Air}}(\varepsilon)$ given by Eq.(1).

The variance of $X_{max}$ under new physics effects can be estimated through

$\text{Var}(X_{max,new})=\text{Var}(X_{\text{int,
new}})+\text{Var}(X_{\text{long,new}})$ (18)

using Eqs. (11) and (14), yielding

Var $\displaystyle(X_{max,new})=\sigma_{X_{\text{max}}}^{2}$
$\displaystyle=\frac{m^{2}}{[\sigma_{\text{p-Air}}(\varepsilon)+\delta\,\beta\,\log\varepsilon]^{2}}+\frac{\text{Var}(X_{long})}{n(\varepsilon)}\,,$
(19)

which again has onle one free parameter: $\delta$. Clearly, $\delta$ can be
optimized by comparing Eq. (II.2.4) with Auger data on
$\sigma_{X_{\text{max}}}$.

###  II.3 Muons

The muonic part of EAS is also a problem in UHECR physics. Muons are produced
mainly when charged pions or kaons decay, indicating hadronic interactions.
They have a large mean free path and consequently their journey in the
atmosphere is mostly undisturbed. Upon reaching the ground, muons can be
detected via the Cherenkov radiation they produce inside the water tanks
comprising the surface array of hybrid experiments such as Auger and TA. Thus
the energy and spacial distribution of muons can be measured.

A useful parameter used to compare experiment with simulations is the ratio

$R_{\mu}=\cfrac{N_{\mu}}{N_{\mu,19}}$ (20)

where $N_{\mu}$ is the number of muons detected on the ground. The reference
parameter $N_{\mu,19}=2.148\cdot 10^{7}$ is inferred from simulations assuming
proton as the primary particle with energy of $10^{19}$ eV, taking into
consideration the detector’s response for muons above $0.3$ GeV 333This is the
Cherenkov threshold for Auger’s water tanks. that reach Auger’s site at an
altitude of $1425$ m with inclination of $\theta=60^{\circ}$. Aab  _et al._
[58] report that $N_{\mu,19}$ does not depend strongly on the selected high
energy model. They quote a $\sim 11\%$ systematic error introduced in this
way. The reason for selecting so inclined showers is because for inclination
$\theta>60$, EAS are dominated mostly by muons, since the electromagnetic part
is absorbed by the atmosphere. This is in fact the inclination we used in our
simulations.

Auger reports that the average of that ratio depends on energy as

$\langle R_{\mu}\rangle=a\Bigg{(}\cfrac{E}{10^{19}\text{ eV}}\Bigg{)}^{b}$
(21)

where $a=1.841$ and $b=1.029$. They also report a relative standard deviation
of

$\sigma[R_{\mu}]/\langle R_{\mu}\rangle=0.136\,.$ (22)

Auger observes $25-40\%$ more muons at $10^{19}$ eV than high energy models
predict, assuming proton as a primary. However, muons are observed to be in
overabundance even if a heavier composition, consistent with the one that
would produce the observed $\langle X_{\text{max}}\rangle(\epsilon)$ is
assumed for the primaries.

In the scenario where proton interactions change above $10^{18}$ eV
introducing new physics, the muon production will change due to the increase
of multiplicity for the secondary particles.

The total number of muons on the ground will be the sum of muons that each
component-shower produces. Note that each component-shower does not produce
the same amount of muons.

We then expect the average number of muons on the ground to be

$\langle N_{\mu,\;new}(\varepsilon)\rangle=n(\varepsilon)\;\langle
N_{\mu}(\varepsilon/n)\rangle\,.$ (23)

###  II.4 CORSIKA simulations

We now turn to EAS simulations. We wish to: (a) test whether the
implementation, in simulations, of new physics as discussed in §II.1 will
produce the same behavior as our analytic approximations; and (b) determine
the optimal phenomenological parametrization (cross section and the
multiplicity, quantified by $\delta$ and $n(\varepsilon)$ in our description)
that any new proton-air interaction must exhibit in order for EAS to produce
the observed data in Auger if the composition of primaries is to remain light
up to the highest energies.

We simulated showers induced by primaries with energies in the range
$10^{17}-10^{20}$ eV with step in $\log E$ of $0.1$. At each energy bin, we
performed $1000$ EAS simulations. We simulated EAS with the first collision
treated either with SM extrapolations, or with new physics as per our
phenomenological implementation. We performed SM EAS simulations for proton
primaries with $E<E_{\rm th}=10^{18}{\rm\,eV}$, and for the heavier Galactic
primaries (see next section), since their per-nucleon kinetic energy never
exceeds $10^{18}{\rm\,eV}$. We implemented new physics for all proton
primaries with $E>10^{18}{\rm\,eV}$.

Figure 1: CR energy spectrum (flattened by $E^{3}$). Blue: Auger Collaboration
data. Green: Galactic CR flux model, assuming power law at low energies with
an exponential cutoff at $10^{17.9}$ eV. Red: extragalactic CR flux, obtained
as the difference between observed data and Galactic CR. Inset: fraction of
Galactic to total CR flux

CORSIKA EAS simulations using either EPOS-LHC or QGSJETII-04 yielded our SM
results. For low-energy interactions we used FLUKA. We also used the CONEX
hybrid scheme [59] which decreases the simulation time dramatically. Each
simulation generates three output files. The first records the energy
deposited as a function of depth in the atmosphere. We fitted a Geisser-Hillas
function

$\displaystyle\cfrac{dE}{dX}(X)=N_{max}\Bigg{(}\cfrac{X-W_{0}}{X_{long}-W_{0}}\Bigg{)}^{\frac{X_{long}-W_{0}}{\lambda}}$
$\displaystyle\times\exp\Bigg{(}\cfrac{X_{long}-X}{\lambda}\Bigg{)}$ (24)

to the simulated data to estimate the shower longitudinal-development maximum,
$X_{\text{long}}$, for each shower. For each energy bin, we then calculated
the average value of $X_{\text{long}}$ and its variance.

The second output file records information on the cross section of the primary
with the atmosphere and from it we calculated $\langle X_{\text{int}}\rangle$.
The second output file also contains the number of muons detected on the
ground. We fitted a convolution of a Gaussian with an exponential to calculate
the average value and the relevant variance of the number of muons on the
ground. The third output file (stack file) contains information about the
secondary particles produced after the first interaction. The results from
this run are the Standard Model predictions from extrapolated models.

For the new-physics simulations, we used the following approach. For a given
value of $\delta$, we first calculated the multiplicity according to
Eq.(II.2.4) at each energy bin. We then combined stack files from the same
energy bin (produced by SM simulations) according to the calculated
multiplicity: we rounded $n(\varepsilon)$ to the nearest integer, and we
combined as many stack files, accounting for energy and momentum conservation.
To this end, we divided the energy and momentum of each particle by the number
of stacked files. We then used the combined stack files as input to CORSIKA
and continued the simulation of the EAS, obtaining data files for the energy
deposition as a function of depth in the atmosphere and the muon number on the
ground. The cross section was calculated from Eq.(2).

###  II.5 Cosmic Ray Flux

At energies above $10^{18.3}$ eV, both Auger and TA observe a dipole
distribution of CR uncorrelated with the galactic plane [60, 61]. This
indicates that above this energy, CR are of extra-galactic origin. The energy
at which this transition takes place is an important input for our
calculations. The reason is that Galactic CR are heavy particles and thus the
the energy per nucleon will be below the threshold of new physics. As a result
we need not apply any new physics corrections to their EAS simulations.

We model this transition using a simple, phenomenological approach, based on
three assumptions: (1) That above $10^{17}$ eV cosmic rays consist of a
single, fixed-composition Galactic component, and a single, fixed-composition
extragalactic component. The energy-dependent fraction of Galactic CR is
$f(\varepsilon)$. The fraction of extragalactic CR is then $1-f(\varepsilon)$.
(2) That above $10^{17}$ eV we can model the Galactic CR spectrum
(differential particle flux $J(\varepsilon)$) as a power law of slope
$-\gamma_{G}$, cutting off exponentially at a a characteristic energy
$\varepsilon_{G}$, corresponding to the maximum energy of Galactic CR
accelerators:

$J_{G}(\varepsilon)=J_{G,0}\Bigg{(}\cfrac{\varepsilon}{\varepsilon_{17.5}}\Bigg{)}^{-\gamma_{G}}\exp{\Big{(}-\varepsilon/\varepsilon_{G}\Big{)}}$
(25)

where $\varepsilon_{17.5}=10^{17.5}\text{eV}/E_{\text{th}}=10^{-0.5}$. (3)
That CR at energies lower than those where losses (either e+e- or pion
photoproduction) become important, the extragalactic CR flux
$J_{EG}(\varepsilon)$ is a single power law.

At low energies ($E<10^{17.5}$ eV), Auger data constrain $J_{G,0}=4.1\times
10^{-15}\,(\text{km}^{2}\;\text{eV}\;\text{yr}\;\text{sr})^{-1}$, and
$\gamma_{G}=2.9$. By virtue of our third assumption above, we can also
constrain $\varepsilon_{G}$ by demanding that, for
$10^{17.5}\text{eV}/E_{\text{th}}<\varepsilon<10^{18.2}\text{eV}/E_{\text{th}}$
the extragalactic spectrum $J_{EG}(\varepsilon)=J_{\text{total,
Auger}}-J_{G}(\varepsilon)$ is consistent with a single power law. We thus
find $\varepsilon_{G}=10^{17.9}\text{eV}/E_{\text{th}}$, which results in an
extragalactic spectrum consistent with
$J_{EG}(\varepsilon)\propto\varepsilon^{-2.0}$ between $10^{17.5}$ and
$10^{18.2}$ eV (see Fig. 1).

The resulting Galactic CR fraction,
$f(\varepsilon)=J_{G}(\varepsilon)/J_{\text{total,Auger}}(\varepsilon)$ is
shown in the inset of Fig. 1. We note that under the three assumptions adopted
here, extragalactic CR are found to dominate already at $10^{18}$ eV, the
composition at $10^{18.5}$ is light, and the ankle must be an e+e- ”dip”. In
this simple scenario, the probability density function of $X_{\text{max}}$
will be

$p(X_{\text{max}})=f\;p_{G}(X_{\text{max}})+(1-f)p_{EG}(X_{\text{max}})$ (26)

leading to an average shower maximum

$\langle X_{\text{max}}\rangle=f\;\langle
X_{\text{max}}\rangle_{G}+(1-f)\langle X_{\text{max}}\rangle_{EG}$ (27)

and its variance

$\displaystyle Var(X_{\text{max}})$ $\displaystyle=f\;Var(X_{\text{max,G}})$
$\displaystyle+(1-f)\;Var(X_{\text{max,EG}})$ $\displaystyle+f(1-f)\;(\langle
X_{\text{max}}\rangle_{G}-\langle X_{\text{max}}\rangle_{EG})^{2}$ (28)

with subscripts G and EG referring to the Galactic and extragalactic
populations respectively.

At energies around $10^{17}$ eV, CR are mainly of Galactic origin. In this
paper we assume Galactic CR to be one type of nucleus for simplicity. We
assume the Galactic component to be helium for simplicity – lighter than PT09,
who had assumed carbon. Ultimately, the Galactic component should ideally be
simulated with appropriate mixed composition, with each species cutting off at
different energies according to its charged, as per the Hillas criterion [10].

##  III Results

[QGSJETII-04]  [EPOS LHC]

Figure 2: Average shower maximum as a function of energy. Standard model
extrapolations (grey) through QGSJETII-04 (upper) and EPOS LHC (lower) do not
agree with observations from Auger Observatory (blue) above
$E_{\text{th}}=10^{18.5}$ GeV. When we alter the way the cross section and the
first-interaction product multiplicity scale with energy as in Eqs. (4) and
(II.2.4), EAS simulations with protons as a primaries (green and red)
reproduce the observed data well at the highest energies.The simplifying
assumption of a single component Galactic CR is the reason our simulations
deviate from observational data at lower energies.

[QGSJETII-04]  [EPOS LHC]

Figure 3: Standard deviation of shower maxima as a function of energy.
Although the mean of $X_{\text{max}}$ does not depend on the parameter
$\delta$ (see Fig. 2), its standard deviation does. We obtain the best
agreement between simulated $\sigma_{X_{\text{max}}}$ and Auger observations
for $\delta$ near $6$ for QGSJETII-04 (upper) and for $\delta$ near $4$ for
EPOS-LHC (lower). Further increment of the $\delta$ parameter does not result
in significant changes at high energies, as $\sigma_{X_{\text{max}}}$ quickly
reaches an asymptotic behavior.  Figure 4: Agreement between new-physics EAS
simulations and Auger data for $\sigma_{X_{\text{max}}}$ and
$E>10^{18.5}{\rm\,eV}$, quantified through the reduced $\chi^{2}$ statistic,
as a function of the value of the $\delta$ parameter. Overall, EPOS-LHC
simulations produce results more consistent with observational data. To find
the position of the minimum $\chi^{2}$, we perform a parabolic fit to the
datapoints. The locations of the minima are at $\delta=4.8$ for EPOS-LHC and
at $\delta=6.7$ for QGSJETII-04.

[QGSJETII-04]  [EPOS LHC]

Figure 5: Fractional difference of muons on the ground predicted by
simulations relative to the observed data at the Auger Observatory. When new
physics sets in above $E_{\text{th}}$, the number of muons change (green and
red), due to the change in product multiplicity (hence, independently of
$\delta$). This change does not fully reconcile the simulated number of muons
with observations. However, unlike the SM prediction (grey), the discrepancy
does not increase with energy above $10^{18.5}$ eV, but rather stabilizes
around $30\%$ ($40\%$) for EPOS-LHC (QGSJETII-04).

[QGSJETII-04]  [EPOS LHC]

Figure 6: Standard deviation over average of the muon ratio as a function of
energy. When new physics sets in above $E_{\text{th}}$, this quantity becomes
more consistent with observed data than SM predictions (grey), and reproduces
the trend with energy observed by Auger.

We performed CORSIKA simulations as described in §II.4 for
$\delta=0,2.9,3.5,4,6,8$ and $10$, treating the Galactic-to-extragalactic
transition as described in §II.5. The simulation results for $\langle
X_{\text{max}}\rangle$ are (by construction) in excellent agreement with Auger
data444Throughout this paper, we use error bars for systematic uncertainties,
and shaded areas to indicate 1, 2, and 3$\sigma$ statistical uncertainties.
for all values of $\delta$ (Fig. 2). Our simulations are also a much improved
fit to the $\sigma_{x_{\text{max}}}$ data even for $\delta=0$ (no change in
cross section, only multiplicity increases). For higher $\delta$ the agreement
improves further and is optimal for $\delta$ between 4 and 8 (Fig. 3). In
Figs. 2 and 3 we did not plot results for all simulated values of $\delta$ in
order to keep the figures legible.

Figs. 2 and 3 show a deviation of the simulation results from observations at
low energies. We expect that this discrepancy is due to the assumption of a
Galactic component consisting purely of helium. A more reasonable assumption
is a mixture of helium, carbon and oxygen, with ratios that depend on the
energy. However, since the details of Galactic CR composition have little to
no impact on the (dis)agreement between theory and observations of
$X_{\text{max}}$ at the highest energies, we chose, for simplicity, not to
focus on the modeling of the Galactic component in this work. Instead, we
implement the very simple recipe described above aiming only to show the
direction in which $\langle X_{\rm max}\rangle$ and $\sigma_{X_{\rm max}}$
will change due to the Galactic-to-extragalactic transition. We plan to return
to this problem and relax the assumption of a single-species Galactic
component in a future publication.

We next explore the agreement between new-physics $\sigma(X_{\text{max}})$ and
observed Auger data, quantifying it by means of the reduced $\chi^{2}$
statistic. Since we have treated the Galactic component in a very approximate
way, we only use energies at which extragalactic protons have fully dominated
the UHECR flux ($E>10^{18.5}$ eV) for this comparison. In Fig.(4) we plot the
reduced $\chi^{2}$ as a function of $\delta$ for our EAS simulation results
(datapoints), as well as for our analytic approximation for
$\sigma_{X_{\text{max}}}$ given by Eq. (II.2.4) (solid lines). We calculate
the position of minimum $\chi^{2}$ for our simulations by fitting parabolas
(dashed lines) to our datapoints. EPOS-LHC has a minimum of $\chi^{2}=1.5$ at
$\delta=4.8$ and QGSJETII-04 has a minimum of $\chi^{2}=5.2$ at $\delta=6.7$.
EPOS-LHC has overall better performance. At their minimum $\chi^{2}$ the cross
section at an energy of $10^{19}$ GV rises to $830$ mb ($788$ mb) for EPOS-LHC
(QGSJETII-04). Furthermore at the same energy the multiplicity has increased
by a factor of $3$ ($2$).

In Figs. 5 and 6 we show how the changes in cross-section and multiplicity in
our scenario affect the number of muons measured on the ground. Fig. 5 shows
the fractional difference of muons on the ground from simulations relative to
the data observed by Auger. The new-physics scenario produces more muons than
the SM prediction. Although this is an improvement, still a deviation between
$30-37\%$ from observational data persists. The reason is that the muon
production does not depend strongly on the overall multiplicity of the
secondary particles but rather on the ratio of pions and kaons produced after
the first interaction [62]. The fraction of such particles in the products of
the first collision does not change in the implementation of new physics we
have considered here. However, such a change in the charged-particle ratio may
be an important feature in any specific new-physics model that attempts to
fully explain the UHECR composition problem, including the muon problem. We do
however point out that the residual discrepancy between simulated and observed
number of muons in our current implementation is constant with energy, unlike
the SM predictions, where the discrepancy increases with energy. Our
implementation additionally produces a variance for the muon number that is
consistent with observed data (Fig.6), including the correct trend with
energy. We plan to revisit this issue in a future publication, adding the
possibility of a change in the charged-particle ratio.

##  IV Summary and Conclusions

We performed simulations of EAS with CORSIKA, appropriately modified for
interactions above $E_{\text{th}}=10^{18}\text{ eV}$ to feature an increased
proton-air cross-section and first-collision product multiplicity. Ww have
parameterized the increase in cross section through a parameter $\delta$,
defined as the fractional increase of the coefficient of logarithmic growth
with energy of the proton-air cross-section with respect to its standard-model
value (see Eq. 4). We have parameterized the increase in product multiplicity
through a function $n(\varepsilon)$, defined as the ratio of first-collision
products over their SM-predicted number. By demanding that Auger observations
of $\langle X_{\text{max}}\rangle$ as a function of energy above $10^{18.7}$
eV are reproduced under the assumption that all UHECR at these energies are
protons, we can determine $n(\varepsilon)$ for any given value of $\delta$
(see Eq. II.2.4). This then leaves $\delta$ as the only parameter in our
description.

We have shown that these modifications to hadronic interactions at energies
above $10^{18}$ eV are sufficient to reproduce observations of
$X_{\text{max}}$. The growth of $\langle X_{\text{max}}\rangle$ is reproduced
for any value of $\delta$ (by construction), and $\sigma_{X_{\text{max}}}$ is
best reproduced for $\delta$ between 4 and 7. If QGSJETII-04 is used for SM
modeling of high-energy hadronic interactions, the optimal value for $\delta$
is 6.7. For EPOS-LHC, the optimal value for $\delta$ is 4.8. Epos-LHC with
$\delta=4.8$ provides the best overall fit to Auger $X_{\text{max}}$ data at
$E>10^{18.5}$ eV. In each case, product multiplicity increases by a factor
$n(\varepsilon)$ given by Eq. (II.2.4).

These results provide phenomenological constraints on the properties (cross
section, multiplicity) of any new effect beyond the SM that may set it for
collisions at CM energies exceeding $\sim 50$ TeV, if such an effect is to be
held responsible for the change of behavior of the $X_{\text{max}}$
distribution observed by Auger for collisions of primaries more energetic than
$10^{18.7}$ eV with the atmosphere.

As far as the muon problem is concerned, we found that although the change in
multiplicity we have investigated here does improve the agreement between
Auger data and EAS simulations, it is not by itself sufficient to fully
resolve the discrepancy. For this reason, we speculate that any new effect
setting in at 50 TeV may also induce a change in the charged-particle ratio,
which would further increase the number of muons produced and detected on the
ground. We plan to investigate this possibility in a future publication.

Current and planned advances in the tomographic mapping of the Galactic
magnetic field through local measurements [63, 64, 29, 65, 66, 67] are
expected to make feasible an electromagnetic determination of the charge of
UHECR of the highest energies in the near future [68]. If such studies provide
unequivocal evidence that UHECR at the highest energies are indeed protons,
then this will be a strong argument in favor of new physics setting in for
hadronic interactions at CM energies above 50 TeV, with the phenomenology of
any such new effect exhibiting the behavior calculated in this work.

###### Acknowledgements.

S.R. and V.P. would like to dedicate this work to the memory of our friend and
collaborator Theodore Tomaras whom we have lost way too soon. We will always
fondly remember all the exciting and fruitful scientific discussions and
debates that we have had and we will be sorely missing the ones that will now
never take place. We thank Alan Watson, Nicusor Arsene, Konstantina Dolapsaki,
Andreas Tersenov, and Christos Litos for helpful comments and discussions that
improved this manuscript, and Dieter Heck for valuable feedback on the use of
CORSIKA. This work was supported by the Hellenic Foundation for Research and
Innovation (H.F.R.I.) under the “First Call for H.F.R.I. Research Projects to
support Faculty members and Researchers and the procurement of high-cost
research equipment grant” (Project 1552 CIRCE). V.P. acknowledges support from
the Foundation of Research and Technology - Hellas Synergy Grants Program
through project MagMASim, jointly implemented by the Institute of Astrophysics
and the Institute of Applied and Computational Mathematics.

## References

  * Hess [1912] V. F. Hess, Über Beobachtungen der durchdringenden Strahlung bei sieben Freiballonfahrten, Phys. Z. 13, 1084 (1912).
  * Abraham and et al [2004] J. Abraham and et al, Properties and performance of the prototype instrument for the Pierre Auger Observatory, Nuclear Instruments and Methods in Physics Research A 523, 50 (2004).
  * Abu-Zayyad and et al [2012] T. Abu-Zayyad and et al, The surface detector array of the Telescope Array experiment, Nuclear Instruments and Methods in Physics Research A 689, 87 (2012), arXiv:1201.4964 [astro-ph.IM] .
  * Sreekumar and et al [1993] P. Sreekumar and et al, Constraints on the cosmic rays in the Small Magellanic Cloud, Phys. Rev. Lett. 70, 127 (1993).
  * Ackermann and et al [2016] M. Ackermann and et al, Deep view of the Large Magellanic Cloud with six years of Fermi-LAT observations, aap 586, A71 (2016), arXiv:1509.06903 [astro-ph.HE] .
  * Ackermann and et al [2017] M. Ackermann and et al, Observations of M31 and M33 with the Fermi Large Area Telescope: A Galactic Center Excess in Andromeda?, apj 836, 208 (2017), arXiv:1702.08602 [astro-ph.HE] .
  * Pierre Auger Collaboration and et al [2017] A. Pierre Auger Collaboration, Aab and et al, Observation of a large-scale anisotropy in the arrival directions of cosmic rays above 8 × 1018 eV, Science 357, 1266 (2017), arXiv:1709.07321 [astro-ph.HE] .
  * Gabici  _et al._ [2019] S. Gabici, C. Evoli, D. Gaggero, P. Lipari, P. Mertsch, E. Orlando, A. Strong, and A. Vittino, The origin of Galactic cosmic rays: Challenges to the standard paradigm, International Journal of Modern Physics D 28, 1930022-339 (2019), arXiv:1903.11584 [astro-ph.HE] .
  * Chiavassa and et al [2019] A. Chiavassa and et al, Summary of the main results of the KASCADE and KASCADE-Grande experiments, in  _European Physical Journal Web of Conferences_, European Physical Journal Web of Conferences, Vol. 208 (2019) p. 03002.
  * Hillas [1984] A. M. Hillas, The Origin of Ultra-High-Energy Cosmic Rays, araa 22, 425 (1984).
  * Kotera and Olinto [2011] K. Kotera and A. V. Olinto, The Astrophysics of Ultrahigh-Energy Cosmic Rays, araa 49, 119 (2011), arXiv:1101.4256 [astro-ph.HE] .
  * Aab  _et al._ [2017] A. Aab _et al._ (Pierre Auger), Combined fit of spectrum and composition data as measured by the Pierre Auger Observatory, JCAP 04, 038, [Erratum: JCAP 03, E02 (2018)], arXiv:1612.07155 [astro-ph.HE] .
  * Haungs, A.  _et al._ [2017] Haungs, A., Apel, W.D., Arteaga-Velázquez, J.C., Bekk, K., Bertaina, M., Blümer, J., Bozdog, H., Brancus, I.M., Cantoni, E., Chiavassa, A., Cossavella, F., Daumiller, K., Souza, V. de, Pierro, F. Di, Doll, P., Engel, R., Fuhrmann, D., Gherghel-Lascu, A., Gils, H.J., Glasstetter, R., Grupen, C., Heck, D., Hörandel, J.R., Huege, T., Kampert, K.-H., Kang, D., Klages, H.O., Link, K., Luczak, P., Mathes, H.J., Mayer, H.J., Milke, J., Mitrica, B., Morello, C., Oehlschläger, J., Ostapchenko, S., Pierog, T., Rebel, H., Roth, M., Schieler, H., Schoo, S., Schröder, F.G., Sima, O., Toma, G., Trinchero, G.C., Ulrich, H., Weindl, A., Wochele, J., and Zabierowski, J., Kascade-grande: Composition studies in the view of the post-lhc hadronic interaction models, EPJ Web Conf. 145, 13001 (2017).
  * Bergman [2020] D. Bergman (Telescope Array), Telescope Array: Latest Results and Expansion Plans, J. Phys. Conf. Ser. 1468, 012078 (2020).
  * Abraham and et al [2010] P. Abraham, J. and et al, Measurement of the energy spectrum of cosmic rays above 1018 eV using the Pierre Auger Observatory, Physics Letters B 685, 239 (2010), arXiv:1002.1975 [astro-ph.HE] .
  * Nagano  _et al._ [1992] M. Nagano, M. Teshima, Y. Matsubara, H. Y. Dai, T. Hara, N. Hayashida, M. Honda, H. Ohoka, and S. Yoshida, Energy spectrum of primary cosmic rays above $10^{17}$ ev determined from the extensive air shower experiment at akeno, Journal of Physics G: Nuclear and Particle Physics 18, 423 (1992).
  * Allard  _et al._ [2005] D. Allard, E. Parizot, A. V. Olinto, E. Khan, and S. Goriely, UHE nuclei propagation and the interpretation of the ankle in the cosmic-ray spectrum, aap 443, L29 (2005), arXiv:astro-ph/0505566 [astro-ph] .
  * Allard  _et al._ [2007] D. Allard, E. Parizot, and A. V. Olinto, On the transition from galactic to extragalactic cosmic-rays: Spectral and composition features from two opposite scenarios, Astroparticle Physics 27, 61 (2007), arXiv:astro-ph/0512345 [astro-ph] .
  * Berezinskii and Grigor’eva [1988] V. S. Berezinskii and S. I. Grigor’eva, A bump in the ultra-high energy cosmic ray spectrum, aap 199, 1 (1988).
  * Berezinsky  _et al._ [2006] V. Berezinsky, A. Gazizov, and S. Grigorieva, On astrophysical solution to ultrahigh energy cosmic rays, Phys. Rev. D 74, 043005 (2006), arXiv:hep-ph/0204357 [hep-ph] .
  * Aloisio  _et al._ [2014] R. Aloisio, V. Berezinsky, and P. Blasi, Ultra high energy cosmic rays: implications of Auger data for source spectra and chemical composition, jcap 2014, 020 (2014), arXiv:1312.7459 [astro-ph.HE] .
  * The Pierre Auger Collaboration and et al [2019] A. The Pierre Auger Collaboration, Aab and et al, The Pierre Auger Observatory: Contributions to the 36th International Cosmic Ray Conference (ICRC 2019), arXiv e-prints , arXiv:1909.09073 (2019), arXiv:1909.09073 [astro-ph.HE] .
  * Aab and et al [2020] A. Aab and et al (The Pierre Auger Collaboration), Features of the energy spectrum of cosmic rays above $2.5\times{}{10}^{18}\text{ }\text{ }\mathrm{eV}$ using the pierre auger observatory, Phys. Rev. Lett. 125, 121106 (2020).
  * di Matteo  _et al._ [2021] A. di Matteo, , and et al, UHECR arrival directions in the latest data from the original Auger and TA surface detectors and nearby galaxies, PoS ICRC2021, 308 (2021).
  * Biteau  _et al._ [2021] J. Biteau, , and et al, The ultra-high-energy cosmic-ray sky above 32 EeV viewed from the Pierre Auger Observatory, PoS ICRC2021, 307 (2021).
  * Abbasi and et al [2014] R. U. Abbasi and et al, Indications of Intermediate-scale Anisotropy of Cosmic Rays with Energy Greater Than 57 EeV in the Northern Sky Measured with the Surface Detector of the Telescope Array Experiment, apjl 790, L21 (2014), arXiv:1404.5890 [astro-ph.HE] .
  * Abbasi and et al [2020] R. U. Abbasi and et al, Evidence for a Supergalactic Structure of Magnetic Deflection Multiplets of Ultra-high-energy Cosmic Rays, apj 899, 86 (2020).
  * Jansson and Farrar [2012] R. Jansson and G. R. Farrar, The Galactic Magnetic Field, apjl 761, L11 (2012), arXiv:1210.7820 [astro-ph.GA] .
  * Tritsis  _et al._ [2019] A. Tritsis, C. Federrath, and V. Pavlidou, Magnetic Field Tomography in Two Clouds toward Ursa Major Using H I Fibers, apj 873, 38 (2019), arXiv:1810.00231 [astro-ph.GA] .
  * Puget  _et al._ [1976] J. L. Puget, F. W. Stecker, and J. H. Bredekamp, Photonuclear interactions of ultrahigh energy cosmic rays and their astrophysical consequences., apj 205, 638 (1976).
  * Stecker and Salamon [1999] F. W. Stecker and M. H. Salamon, Photodisintegration of ultrahigh-energy cosmic rays: A New determination, Astrophys. J. 512, 521 (1999), arXiv:astro-ph/9808110 .
  * Shaham and Piran [2013] N. Shaham and T. Piran, Implications of the penetration depth of ultrahigh-energy cosmic rays on physics at 100 tev, Phys. Rev. Lett. 110, 021101 (2013).
  * Aab and et al [2014] A. Aab and 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 (2014).
  * Aab and et al [2017] A. Aab and et al, Combined fit of spectrum and composition data as measured by the Pierre Auger Observatory, jcap 2017, 038 (2017), arXiv:1612.07155 [astro-ph.HE] .
  * Arsene [2021] N. Arsene, Mass Composition of UHECRs from Xmax Distributions Recorded by the Pierre Auger and Telescope Array Observatories, Universe 7, 321 (2021), arXiv:2109.03626 [hep-ph] .
  * Guido and et al [2022] T. Guido, E. and et al, Combined fit of the energy spectrum and mass composition across the ankle with the data measured at the Pierre Auger Observatory, in  _37th International Cosmic Ray Conference. 12-23 July 2021\. Berlin_ (2022) p. 311.
  * Unger  _et al._ [2015] M. Unger, G. R. Farrar, and L. A. Anchordoqui, Origin of the ankle in the ultrahigh energy cosmic ray spectrum, and of the extragalactic protons below it, Phys. Rev. D 92, 123001 (2015).
  * Globus  _et al._ [2015] N. Globus, D. Allard, and E. Parizot, A complete model of the cosmic ray spectrum and composition across the galactic to extragalactic transition, Phys. Rev. D 92, 021302 (2015).
  * Sánchez [2019] F. A. Sánchez, The muon component of extensive air showers above 10(17.5) eV measured with the Pierre Auger Observatory, in  _36th International Cosmic Ray Conference (ICRC2019)_ , International Cosmic Ray Conference, Vol. 36 (2019) p. 411.
  * Aab and et al [2016] A. Aab and 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 (2016).
  * Veberic [2017] D. Veberic, ed.,  _The Pierre Auger Observatory: Contributions to the 35th International Cosmic Ray Conference (ICRC 2017)_ (2017) arXiv:1708.06592 [astro-ph.HE] .
  * Aab  _et al._ [2018a] A. Aab _et al._ (Pierre Auger), An Indication of anisotropy in arrival directions of ultra-high-energy cosmic rays through comparison to the flux pattern of extragalactic gamma-ray sources, Astrophys. J. Lett. 853, L29 (2018a), arXiv:1801.06160 [astro-ph.HE] .
  * Farrar and Allen [2013] G. R. Farrar and J. D. Allen, A new physical phenomenon in ultra-high energy collisions, EPJ Web Conf. 53, 07007 (2013), arXiv:1307.2322 [hep-ph] .
  * Anchordoqui  _et al._ [2017] L. A. Anchordoqui, H. Goldberg, and T. J. Weiler, Strange fireball as an explanation of the muon excess in auger data, Phys. Rev. D 95, 063005 (2017).
  * Tomar [2017] G. Tomar, Lorentz invariance violation as an explanation of the muon excess in auger data, Phys. Rev. D 95, 095035 (2017).
  * Pavlidou and Tomaras [2019] V. Pavlidou and T. Tomaras, What do the highest-energy cosmic-ray data suggest about possible new physics around 50 TeV?, Phys. Rev. D 99, 123016 (2019), arXiv:1802.04806 [astro-ph.HE] .
  * Heck  _et al._ [1998] D. Heck, J. Knapp, J. N. Capdevielle, G. Schatz, and T. Thouw,  _CORSIKA: a Monte Carlo code to simulate extensive air showers._ (1998).
  * 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, Physical Review C 92, 10.1103/physrevc.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).
  * Ferrari  _et al._ [2005] A. Ferrari, P. R. Sala, M. /CERN /INFN, A. Fasso, /SLAC, J. Ranft, and S. U., Fluka: A multi-particle transport code 10.2172/877507 (2005).
  * Ulrich  _et al._ [2011] R. Ulrich, R. Engel, and M. Unger, Hadronic Multiparticle Production at Ultra-High Energies and Extensive Air Showers, Phys. Rev. D 83, 054026 (2011), arXiv:1010.4310 [hep-ph] .
  * Mohanty and Serreau [2005] B. Mohanty and J. Serreau, Disoriented chiral condensate: Theory and experiment [review article], physrep 414, 263 (2005), arXiv:hep-ph/0504154 [hep-ph] .
  * Cafarella  _et al._ [2005] A. Cafarella, C. Corianò, and T. N. Tomaras, Cosmic ray signals from mini black holes in models with extra dimensions: an analytical/Monte Carlo study, Journal of High Energy Physics 2005, 065 (2005), arXiv:hep-ph/0410358 [hep-ph] .
  * Mironov  _et al._ [2009] A. Mironov, A. Morozov, and T. N. Tomaras, Can Centauros or Chirons BE the First Observations of Evaporating Mini Black Holes?, International Journal of Modern Physics A 24, 4097 (2009), arXiv:hep-ph/0311318 [hep-ph] .
  * Patrignani  _et al._ [2016] C. Patrignani _et al._ (Particle Data Group), Review of Particle Physics, Chin. Phys. C 40, 100001 (2016).
  * Heitler [1936] W. Heitler,  _The quantum theory of radiation_ , International Series of Monographs on Physics, Vol. 5 (Oxford University Press, Oxford, 1936).
  * Matthews [2005] J. Matthews, A heitler model of extensive air showers, Astroparticle Physics 22, 387 (2005).
  * Aab  _et al._ [2015] A. Aab _et al._ (Pierre Auger), Muons in Air Showers at the Pierre Auger Observatory: Mean Number in Highly Inclined Events, Phys. Rev. D 91, 032003 (2015), [Erratum: Phys.Rev.D 91, 059901 (2015)], arXiv:1408.1421 [astro-ph.HE] .
  * Bergmann  _et al._ [2007] T. Bergmann, R. Engel, D. Heck, N. N. Kalmykov, S. Ostapchenko, T. Pierog, T. Thouw, and K. Werner, One-dimensional Hybrid Approach to Extensive Air Shower Simulation, Astropart. Phys. 26, 420 (2007), arXiv:astro-ph/0606564 .
  * Aab  _et al._ [2018b] A. Aab _et al._ (Pierre Auger), Large-scale cosmic-ray anisotropies above 4 EeV measured by the Pierre Auger Observatory, Astrophys. J. 868, 4 (2018b), arXiv:1808.03579 [astro-ph.HE] .
  * Abbasi  _et al._ [2020] R. U. Abbasi _et al._ (Telescope Array), Search for Large-scale Anisotropy on Arrival Directions of Ultra-high-energy Cosmic Rays Observed with the Telescope Array Experiment, Astrophys. J. Lett. 898, L28 (2020), arXiv:2007.00023 [astro-ph.HE] .
  * Dembinski  _et al._ [2021] H. Dembinski _et al._ , The Muon Puzzle in air showers and its connection to the LHC, PoS ICRC2021, 037 (2021).
  * Tassis  _et al._ [2018] K. Tassis, A. N. Ramaprakash, A. C. S. Readhead, S. B. Potter, I. K. Wehus, G. V. Panopoulou, D. Blinov, H. K. Eriksen, B. Hensley, A. Karakci, J. A. Kypriotakis, S. Maharana, E. Ntormousi, V. Pavlidou, T. J. Pearson, and R. Skalidis, PASIPHAE: A high-Galactic-latitude, high-accuracy optopolarimetric survey, arXiv e-prints , arXiv:1810.05652 (2018), arXiv:1810.05652 [astro-ph.IM] .
  * Tritsis  _et al._ [2018] A. Tritsis, C. Federrath, N. Schneider, and K. Tassis, A new method for probing magnetic field strengths from striations in the interstellar medium, mnras 481, 5275 (2018), arXiv:1810.01559 [astro-ph.GA] .
  * Clark and Hensley [2019] S. E. Clark and B. S. Hensley, Mapping the Magnetic Interstellar Medium in Three Dimensions over the Full Sky with Neutral Hydrogen, apj 887, 136 (2019), arXiv:1909.11673 [astro-ph.GA] .
  * Skalidis and Tassis [2021] R. Skalidis and K. Tassis, High-accuracy estimation of magnetic field strength in the interstellar medium from dust polarization, aap 647, A186 (2021), arXiv:2010.15141 [astro-ph.GA] .
  * Skalidis  _et al._ [2021] R. Skalidis, J. Sternberg, J. R. Beattie, V. Pavlidou, and K. Tassis, Why take the square root? An assessment of interstellar magnetic field strength estimation methods, aap 656, A118 (2021), arXiv:2109.10925 [astro-ph.GA] .
  * Magkos and Pavlidou [2019] G. Magkos and V. Pavlidou, Deflections of ultra-high energy cosmic rays by the Milky Way magnetic field: how well can they be corrected?, jcap 2019, 004 (2019), arXiv:1802.03409 [astro-ph.HE] .