# Search for photons above $\text{10}^{\text{19}}$ eV with the surface
detector of the Pierre Auger Observatory

P. Abreu  M. Aglietta  I. Allekotte  K. Almeida Cheminant  A. Almela  J.
Alvarez-Muñiz  J. Ammerman Yebra  G.A. Anastasi  L. Anchordoqui  B. Andrada
S. Andringa  C. Aramo  P.R. Araújo Ferreira  E. Arnone  J. C. Arteaga
Velázquez  H. Asorey  P. Assis  G. Avila  E. Avocone  A.M. Badescu  A.
Bakalova  A. Balaceanu  F. Barbato  J.A. Bellido  C. Berat  M.E. Bertaina  G.
Bhatta  P.L. Biermann  V. Binet  K. Bismark  T. Bister  J. Biteau  J. Blazek
C. Bleve  J. Blümer  M. Boháčová  D. Boncioli  C. Bonifazi  L. Bonneau
Arbeletche  N. Borodai  J. Brack  T. Bretz  P.G. Brichetto Orchera  F.L.
Briechle  P. Buchholz  A. Bueno  S. Buitink  M. Buscemi  M. Büsken  A. Bwembya
K.S. Caballero-Mora  L. Caccianiga  I. Caracas  R. Caruso  A. Castellina  F.
Catalani  G. Cataldi  L. Cazon  M. Cerda  J.A. Chinellato  J. Chudoba  L.
Chytka  R.W. Clay  A.C. Cobos Cerutti  R. Colalillo  A. Coleman  M.R. Coluccia
R. Conceição  A. Condorelli  G. Consolati  F. Contreras  F. Convenga  D.
Correia dos Santos  C.E. Covault  M. Cristinziani  S. Dasso  K. Daumiller
B.R. Dawson  R.M. de Almeida  J. de Jesús  S.J. de Jong  J.R.T. de Mello Neto
I. De Mitri  J. de Oliveira  D. de Oliveira Franco  F. de Palma  V. de Souza
E. De Vito  A. Del Popolo  O. Deligny  L. Deval  A. di Matteo  M. Dobre  C.
Dobrigkeit  J.C. D’Olivo  L.M. Domingues Mendes  R.C. dos Anjos  J. Ebr  M.
Eman  R. Engel  I. Epicoco  M. Erdmann  C.O. Escobar  A. Etchegoyen  H. Falcke
J. Farmer  G. Farrar  A.C. Fauth  N. Fazzini  F. Feldbusch  F. Fenu  B. Fick
J.M. Figueira  A. Filipčič  T. Fitoussi  T. Fodran  T. Fujii  A. Fuster  C.
Galea  C. Galelli  B. García  H. Gemmeke  F. Gesualdi  A. Gherghel-Lascu  P.L.
Ghia  U. Giaccari  M. Giammarchi  J. Glombitza  F. Gobbi  F. Gollan  G. Golup
M. Gómez Berisso  P.F. Gómez Vitale  J.P. Gongora  J.M. González  N. González
I. Goos  D. Góra  A. Gorgi  M. Gottowik  T.D. Grubb  F. Guarino  G.P. Guedes
E. Guido  S. Hahn  P. Hamal  M.R. Hampel  P. Hansen  D. Harari  V.M. Harvey
A. Haungs  T. Hebbeker  D. Heck  C. Hojvat  J.R. Hörandel  P. Horvath  M.
Hrabovský  T. Huege  A. Insolia  P.G. Isar  P. Janecek  J.A. Johnsen  J.
Jurysek  A. Kääpä  K.H. Kampert  B. Keilhauer  A. Khakurdikar  V.V. Kizakke
Covilakam  H.O. Klages  M. Kleifges  J. Kleinfeller  F. Knapp  N. Krohm  N.
Kunka  B.L. Lago  N. Langner  M.A. Leigui de Oliveira  V. Lenok  A. Letessier-
Selvon  I. Lhenry-Yvon  D. Lo Presti  L. Lopes  R. López  L. Lu  Q. Luce  J.P.
Lundquist  A. Machado Payeras  G. Mancarella  D. Mandat  B.C. Manning  J.
Manshanden  P. Mantsch  S. Marafico  F.M. Mariani  A.G. Mariazzi  I.C. Mariş
G. Marsella  D. Martello  S. Martinelli  O. Martínez Bravo  M.A. Martins  M.
Mastrodicasa  H.J. Mathes  J. Matthews  G. Matthiae  E. Mayotte  S. Mayotte
P.O. Mazur  G. Medina-Tanco  D. Melo  A. Menshikov  S. Michal  M.I. Micheletti
L. Miramonti  S. Mollerach  F. Montanet  L. Morejon  C. Morello  A.L. Müller
K. Mulrey  R. Mussa  M. Muzio  W.M. Namasaka  A. Nasr-Esfahani  L. Nellen  G.
Nicora  M. Niculescu-Oglinzanu  M. Niechciol  D. Nitz  I. Norwood  D. Nosek
V. Novotny  L. Nožka  A Nucita  L.A. Núñez  C. Oliveira  M. Palatka  J.
Pallotta  P. Papenbreer  G. Parente  A. Parra  J. Pawlowsky  M. Pech  J.
Pȩkala  R. Pelayo  E.E. Pereira Martins  J. Perez Armand  C. Pérez Bertolli
L. Perrone  S. Petrera  C. Petrucci  T. Pierog  M. Pimenta  M. Platino  B.
Pont  M. Pothast  M. Pourmohammad Shavar  P. Privitera  M. Prouza  A. Puyleart
S. Querchfeld  J. Rautenberg  D. Ravignani  M. Reininghaus  J. Ridky  F. Riehn
M. Risse  V. Rizi  W. Rodrigues de Carvalho  J. Rodriguez Rojo  M.J. Roncoroni
S. Rossoni  M. Roth  E. Roulet  A.C. Rovero  P. Ruehl  A. Saftoiu  M. Saharan
F. Salamida  H. Salazar  G. Salina  J.D. Sanabria Gomez  F. Sánchez  E.M.
Santos  E. Santos  F. Sarazin  R. Sarmento  R. Sato  P. Savina  C.M. Schäfer
V. Scherini  H. Schieler  M. Schimassek  M. Schimp  F. Schlüter  D. Schmidt
O. Scholten  H. Schoorlemmer  P. Schovánek  F.G. Schröder  J. Schulte  T.
Schulz  S.J. Sciutto  M. Scornavacche  A. Segreto  S. Sehgal  S.U.
Shivashankara  G. Sigl  G. Silli  O. Sima  R. Smau  R. Šmída  P. Sommers  J.F.
Soriano  R. Squartini  M. Stadelmaier  D. Stanca  S. Stanič  J. Stasielak  P.
Stassi  M. Straub  A. Streich  M. Suárez-Durán  T. Sudholz  T. Suomijärvi
A.D. Supanitsky  Z. Szadkowski  A. Tapia  C. Taricco  C. Timmermans  O.
Tkachenko  P. Tobiska  C.J. Todero Peixoto  B. Tomé  Z. Torrès  A. Travaini
P. Travnicek  C. Trimarelli  M. Tueros  R. Ulrich  M. Unger  L. Vaclavek  M.
Vacula  J.F. Valdés Galicia  L. Valore  E. Varela  A. Vásquez-Ramírez  D.
Veberič  C. Ventura  I.D. Vergara Quispe  V. Verzi  J. Vicha  J. Vink  S.
Vorobiov  C. Watanabe  A.A. Watson  A. Weindl  L. Wiencke  H. Wilczyński  D.
Wittkowski  B. Wundheiler  A. Yushkov  O. Zapparrata  E. Zas  D. Zavrtanik  M.
Zavrtanik  and L. Zehrer

###### Abstract

We use the surface detector of the Pierre Auger Observatory to search for air
showers initiated by photons with an energy above $10^{19}$ eV. Photons in the
zenith angle range from 30∘ to 60∘ can be identified in the overwhelming
background of showers initiated by charged cosmic rays through the broader
time structure of the signals induced in the water-Cherenkov detectors of the
array and the steeper lateral distribution of shower particles reaching
ground. Applying the search method to data collected between January 2004 and
June 2020, upper limits at 95% CL are set to an $E^{-2}$ diffuse flux of
ultra-high energy photons above $10^{19}$ eV, $2{\times}10^{19}$ eV and
$4{\times}10^{19}$ eV amounting to $2.11{\times}10^{-3}$,
$3.12{\times}10^{-4}$ and $1.72{\times}10^{-4}$ km-2 sr-1 yr-1, respectively.
While the sensitivity of the present search around $2\times 10^{19}$ eV
approaches expectations of cosmogenic photon fluxes in the case of a pure-
proton composition, it is one order of magnitude above those from more
realistic mixed-composition models. The inferred limits have also implications
for the search of super-heavy dark matter that are discussed and illustrated.

##  1 Introduction

Photons with energies above $10^{19}$ eV can be produced by $\pi^{0}$ decays
subsequent to the interactions of ultra-high energy cosmic rays (UHECRs) with
the photon fields or the dust permeating the source environments, or the
background photon fields in the extragalactic space. The resulting fluxes are
attenuated over distances of ${\sim}$10 Mpc by $e^{\pm}$ pair production in
interactions of those photons with the cosmic-background ones, see e.g. [1,
2], so that the detectable volume of photon sources encompasses only the local
Universe. A diffuse photon flux is expected from the cosmogenic photons
produced by UHECR interactions while they are propagating within the
detectable volume.

The cosmogenic photon flux depends on the nature of the UHECRs. The hadrons
that cause the creation of the $\pi^{0}$ mesons, through the process of
resonant photopion reaction, must have energies typically ten times higher
than the secondary photons. Such hadrons can be primary proton CRs, or
secondary ones produced from the photo-disintegration of nuclei interacting
inelastically with a cosmic-background photon, which leads to the production
of nucleons inheriting the energy of the fragmented nucleus divided by its
atomic number. Given the steepening of the UHECR intensity with energy,
photons are thus more efficiently produced above
${10}^{19}\text{\,}\mathrm{e}\mathrm{V}$ by UHE protons. Several mass-
sensitive observables are however providing evidence that the mass composition
of UHECRs is gradually getting heavier above $10^{18.3}$ eV and is, in
particular, not compatible with a pure-proton composition [3, 4, 5, 6]. This
is in line with the absence of copious fluxes of cosmogenic photons and
neutrinos with energies ranging from GeV to EeV, as reported in [7] from the
extragalactic gamma-ray flux at GeV energies and in [8, 9] from neutrino
searches above a hundred of PeV. These results provide important constrains on
the sources of protons at UHE [10, 11, 12, 13, 14]. Even though scenarios
based on a mixed composition of UHECRs are more demanding in terms of exposure
to photons, the search for these emblematic messengers is thus complementary
to complete the multi-messenger approach aimed at understanding the non-
thermal processes producing UHECRs in the Universe.

While the search for photons above ${10}^{19}\text{\,}\mathrm{e}\mathrm{V}$ is
of primary importance to decipher further the origin of UHECRs, the detection
of photons of even higher energies, above
${10}^{20}\text{\,}\mathrm{e}\mathrm{V}$, would open an unexpected window,
revealing either new physics such as Lorentz invariance violation [15, 16, 17,
18, 19, 20, 21] or signatures of axion mixing models [22], or some new
particle acceleration never seen or imagined until now. At the same time, the
detection of a flux of UHE photons could be compelling evidence for dark
matter (DM) composed of super-heavy particles. In cosmological models with an
inflation phase, such particles never at thermal equilibrium could have been
produced at reheating after inflation through mechanisms involving gravitation
[23, 24, 25]. Despite being metastable particles, they can decay through non-
perturbative effects into standard-model particles [26, 23], and hence produce
detectable secondaries such as nucleons and photons. Of particular interest
would thus be the detection of UHE photons from regions of denser DM density
such as the center of our Galaxy. The limits on such photon fluxes translates
into constrains on the lifetimes and masses of DM particles [27, 28, 29, 30,
31].

In this work we update the search for UHE photons above $10^{19}$ eV using the
surface detector (SD) array of the Pierre Auger Observatory [32]. Compared to
previous analyses [33, 34, 35, 36], this work benefits from the increased
exposure cumulated from January 2004 to June 2020, as well as from a refined
search method and data selection. The search for photons presented here
complements and extends previous searches using data from the Pierre Auger
Observatory at lower energies [37, 38]. The paper is organized as follows. In
Section 2, general features of photon induced extensive air showers are
presented focusing on the differences expected with respect to the bulk of
showers initiated by nuclei that constitute the background for the search. The
Pierre Auger Observatory is also briefly described, with a more specific
emphasis on the SD array used in this analysis. The discriminating variables
aiming at identifying photon showers are introduced in Section 3. The Monte
Carlo simulations of photons used in the analysis, the photon energy scale,
and the analysis method combining the discriminating variables to extract
photons from the bulk of events are detailed in Section 4. The results of the
photon search and the upper limits to the diffuse flux of UHE photons are
presented in Section 5. Finally, some astrophysical implications of the
results are discussed in Section 6.

##  2 Photon showers at the Pierre Auger Observatory

The Pierre Auger Observatory is a ground-based instrument designed to detect
the extensive air showers (EAS) induced in the atmosphere by UHECRs. We
briefly discuss here the main features of the photon-induced showers compared
to nucleus-induced ones and how the showers are detected and reconstructed at
the Pierre Auger Observatory. A more detailed description can be found
elsewhere for extensive air showers [39, 40], the Observatory [32] and the
reconstruction of EAS [41].

###  2.1 Main features of photon showers

Figure 1: Main features of photon- and nucleus-induced showers.

Air showers initiated by high-energy photons in the atmosphere differ
significantly from those from nuclei. For a photon-induced shower, the first
interactions and generations are purely electromagnetic, since the radiation
length is more than two orders of magnitude smaller than the mean free path
for photo-nuclear interactions. Yet, the development of the shower is delayed
by the typically small multiplicity of electromagnetic interactions. Thus the
maximum development of the shower is reached at a slant atmospheric depth
$X_{\rm max}$ larger for photon primaries than for nuclei, with a difference
of $\simeq 200\ \rm{g\ cm^{-2}}$ between photons and protons at $10^{19}$ eV
and even larger between photons and heavy nuclei.

The lateral distribution of secondary particles at a given stage of
development is governed by the moderate transverse momentum of the processes
in the cascade and by the mean free path of the particles. Overall, the
steepness of the lateral distribution decreases with the slant depth $X$ so as
to get flatter through the shower development, and the fall-off with the
distance to the axis of the shower depends on the primary mass of the cosmic
rays. At ground level, the steepness is thus relevant to distinguish between
nucleus-induced showers and photon-induced ones.

Since the mean free path for photo-nuclear interactions is much larger than
the radiation length, the transfer of energy to the hadron and muon channels
is reduced hence only a small fraction of the electromagnetic component in a
photon-induced shower is injected into the hadronic cascade. Showers induced
by photons are thus characterized by a lower content of muons: on average,
simulations show that photon showers have nearly one order of magnitude less
muons than proton showers of the same energy.

These main features of photon showers, depicted in Fig. 1, are amplified by
the Landau-Pomeranchuk-Migdal (LPM) effect [42, 43] resulting in a suppression
of the bremsstrahlung and pair-production cross sections.

The picture of UHE photon showers is supplemented by accounting for the
influence of the magnetic field of the Earth, which can allow for the
conversion of photons into an $e^{\pm}$ pair before they enter the upper
atmosphere (“preshowering” effect [44]). The resulting showers are a
superposition of cascades initiated by lower energy electrons and photons,
giving rise to smaller $X_{\rm max}$ values on average.

###  2.2 The Pierre Auger Observatory and the events collected with the
surface detector

The Pierre Auger Observatory is located in the province of Mendoza, Argentina,
at 1400 m a.s.l. – corresponding to an atmospheric overburden of $\simeq$
$875\text{\,}\mathrm{g}\mathrm{/}\mathrm{c}\mathrm{m}^{2}$. It is designed as
a hybrid cosmic-ray detector using two proven techniques to measure the
properties of EAS by observing their longitudinal development in the
atmosphere with a fluorescence detector (FD) and their lateral spread at
ground level with a surface detector (SD) array. The search for photons
presented here makes use of data collected only with the SD, which operates
with $\simeq$100% duty cycle.

The SD array consists of a triangular grid of about 1600 water-Cherenkov
detectors (WCDs), with a spacing of 1500 m, covering a total area of 3000 km2.
Each WCD is a cylinder with a surface area of 10 m2 and a height of 1.2 m,
holding 12 tonnes of ultrapure water viewed by three $9^{\prime\prime}$
photomultipliers (PMTs). These detect the Cherenkov light emitted in water by
charged particles and the $e^{\pm}$ pairs produced in water by secondary
photons reaching the ground. The signals from the PMTs are digitized using 40
MHz 10-bit flash analog-to-digital converters (FADCs). The signals are
normalized to the signal obtained for a vertical muon and expressed in
vertical equivalent muons (VEM). Data are collected in real time by searching
for temporal and spatial coincidences at a minimum of three WCD locations to
build the event triggers. When this occurs, the pulse amplitude and time of
detection of signals are obtained from the FADC data of the PMTs. The data
quality is guaranteed by both an on-line and a long-term continuous monitoring
of the detectors.

The arrival direction of the primary particle initiating the EAS is
reconstructed using the start time of the signals recorded in individual
detectors and is determined with a resolution of $\simeq 1^{\circ}$. To
estimate the energy of the primary particle, the total signal amplitude,
integrated in time, of each triggered detector, $S_{i}$, is used.

The lateral extension of the showers, also known as the lateral distribution
function (LDF), is measured at the ground level. The signal deposit $S_{i}$ in
each WCD is adjusted by scaling the normalisation of an average LDF in a
fitting procedure so as to best reproduce the observed signal amplitudes [41].
The energy estimator is then the signal interpolated at 1000 m from the shower
axis, $S(1000)$. The attenuation of $S(1000)$ in the atmosphere for showers
with same energy but different zenith angle is accounted for with the constant
intensity cut method [45]. The final value of $S(1000)$ is calibrated with the
quasi-calorimetric measurement of the primary energy provided by the
fluorescence detector for a subset of hybrid events [46].

According to the procedure described, the energy $E_{\rm had}$ assigned to
each event is a function $f_{\rm had}(S(1000),\theta)$, whose parameters are
calibrated with data and represents the correct energy scale for nuclear
primaries. Since photons are characterised by an $X_{\rm max}$ and a muonic
content significantly different from the bulk of data, $E_{\rm had}$ is not
suitable to provide an accurate estimate of the energy of photon showers,
resulting indeed in an overestimate of more than a factor two. An alternative
procedure is therefore required to assign the correct photon energy scale, as
explained in Section 4.2.

##  3 Observables from the Surface Detector

Figure 2: A shower induced by a simulated photon of 39 EeV with a zenith angle
of 44∘: (top) the lateral distribution of signals in the WCDs is steeper than
the LDF obtained from data (solid line) and (bottom) the risetime of the
signal in the WCDs is larger than the average data benchmark (solid line).

The main characteristics of an air shower for the identification of the nature
of the primary particle are its $X_{\rm max}$ value and the muonic content
(Section 2.1). While $X_{\rm max}$ can be measured directly with the FD, this
is not possible with the SD that detects the secondary particles of the air
shower reaching the ground. The muon content, as well, cannot be estimated
directly with WCDs, which record an aggregate signal from muons and other
electromagnetic shower components. We use instead two robust data-driven
variables describing the spread in time of the shower front and the steepness
of the lateral distribution of time-integrated signals, sensitive both to the
depth of the shower maximum and the muonic content of the EAS. The strategy to
search for photons is to identify, in the bulk of events detected, showers
that depart significantly from the average behaviour of data in the direction
expected for photon primaries, i.e. events with a larger spread in the arrival
times of the secondaries and a steeper LDF.

###  3.1 Signal risetime

The spread in arrival times of secondary particles in individual WCDs can be
measured through the risetime defined as the time at which the integrated
signal in the FADC time trace rises from 10% to 50% of its total value. The
risetime is increased by a larger contribution of the electromagnetic
component as secondary photons and electrons undergo scattering and
attenuation in the atmosphere, unlike muons, which are concentrated in time
close to the shower front. It also increases when the difference in depth
between $X_{\rm max}$ and the observation level becomes smaller, for
geometrical reasons [47]. Being sensitive to both the deeper $X_{\rm max}$ and
the limited muon content of photon showers, the risetime is a suitable
variable for the search of photons.

A “Data Benchmark” is produced to describe the average risetime of data as a
function of the zenith angle and of the distance to the shower axis, following
a procedure similar to the one described in detail in [47]. A correction for
asymmetries in the observed risetime is obtained from data, accounting for the
fact that for non-vertical air showers, the ground level observation
corresponds to different stages of development (or “age”) of the shower for
geometrical reasons. In the following, we will denote the asymmetry corrected
risetime as $t_{1\mspace{-2.0mu}/\mspace{-1.0mu}2}$. Sampling fluctuations,
$\sigma_{t_{1\mspace{-2.0mu}/\mspace{-1.0mu}2}}$, are also estimated from the
data, using the difference between the measurements of SD doublets (a WCD in
the regular SD grid plus a second one off-grid deployed close to it) or SD
pairs (two WCDs in the same event with similar distance from the axis and
total signal).

The $t_{1\mspace{-2.0mu}/\mspace{-1.0mu}2}^{\rm bench}$ and
$\sigma_{t_{1\mspace{-2.0mu}/\mspace{-1.0mu}2}{}}$ describe the average
thickness of the local shower disk for the bulk of cosmic rays detected with
the SD. For each triggering detector in an event, the quantity

$\delta_{i}=\frac{(t_{1\mspace{-2.0mu}/\mspace{-1.0mu}2}^{i}-t_{1\mspace{-2.0mu}/\mspace{-1.0mu}2}^{\rm
bench})}{\sigma^{i}_{t_{1\mspace{-2.0mu}/\mspace{-1.0mu}2}}}$ (3.1)

is then providing a measurement of the deviation of the risetime in the $i$-th
WCD from the data benchmark in units of expected standard deviation. An
example of $\delta_{i}$ values for one simulated photon shower event is shown
in Fig. 2, where large departures from the benchmark curve are observed.

We can therefore define an observable that measures the departure of an
individual event from the data-averaged behaviour of air-showers as

$\Delta=\frac{1}{N}\sum_{i=1}^{N}{\delta_{i}},$ (3.2)

where $N$ is the number of the triggered detectors in the event. $\Delta$ is
expected to average to zero for data by construction and to be significantly
positive for air showers initiated by photons.

###  3.2 Steepness of the lateral distribution of signals

The reduced muon content of photon showers with respect to data produces, as
already mentioned, a steeper LDF of the signals in the detectors at ground
level. At large distances from the axis, photon showers thus produce typically
smaller signals than expected from the data LDF. This is illustrated in the
top panel of Fig. 2 where the LDF function is

$f_{\mathrm{LDF}}(r)=S(1000)\ \left(\frac{r}{r_{\mathrm{opt}}}\right)^{\beta}\
\left(\frac{r+r_{s}}{r_{\mathrm{opt}}+r_{s}}\right)^{\gamma+\beta},$ (3.3)

with $r_{\rm opt}=1000$ m, $r_{s}=700$ m and $\gamma$ and $\beta$ are
parameterized as a function of $S(1000)$ and $\theta$ to describe the average
behaviour of data [41].

We define an observable $L_{\rm LDF}$ measuring the departure of the observed
signals from the average data LDF as the logarithm of the average ratio
between the SD signals and $f_{\rm LDF}(r)$:

$L_{\rm LDF}=\log_{10}\left(\frac{1}{N}\sum_{i=1}^{N}\frac{S_{i}}{f_{\rm
LDF}(r_{i})}\right),$ (3.4)

where $S_{i}$ is the total signal of the $i$-th detector and $i$ runs over the
$N$ detectors with radial distance from the shower axis $r_{i}>$1000 m, where
the signal for photon showers is expected to be lower than $f_{\rm LDF}$.
$L_{\rm LDF}$ is expected to be close to 0 for data, as by construction the
LDF function is built to describe the data average behaviour, and negative for
photons.

##  4 Analysis method

The data collected with the SD of the Pierre Auger Observatory between January
$1^{\text{st}}$ 2004 and June $30^{\text{th}}$ 2020 are used for the analysis
described here. The expected physical differences between signal and
background are quantified using simulations of air showers initiated by
photons. The background of showers induced by nuclei is not simulated; instead
a fraction of the data set is used as a burn sample to define the selection
for candidate events. The search sample will consist of the remaining events.

###  4.1 Monte Carlo simulations of photon showers

A set of 27,000 photon showers is simulated using CORSIKA [48] with EPOS-LHC
[49] as the high-energy generator of hadronic interactions and FLUKA [50, 51,
52] as the low-energy interaction model. The number of muons at ground level
in simulated showers relies on hadronic interaction properties at very high
energies and in phase-space regions not well covered by accelerator
experiments. The EPOS-LHC generator used here is one of the leading LHC-tuned
high-energy hadronic interaction model predicting in general the largest
number of muons.

The energy distribution of simulated showers follows a $\Phi_{\rm
gen}(E)\propto E^{-1}$ law in the energy range $10^{18.5}-10^{20.5}\,$eV; the
arrival directions are distributed in zenith angle $\theta$ between
$20^{\circ}$ and $70^{\circ}$ according to a $\cos^{2}(\theta)$ distribution
to simulate an isotropic flux impinging on a flat surface at ground. Any
arbitrary energy spectrum other than the one simulated can be reproduced by
weighting the simulated events with $w(E)=\Phi(E)/\Phi_{\rm gen}(E)$.

Both the LPM effect and geomagnetic cascading (preshowering), described in
Section 2.1, are considered in the simulations [53, 54].

The response of the SD array is simulated using the $\overline{\textrm{Off}}$
$\underline{\textrm{line}}$ package [55] providing as output simulated events
in the same format as real events. Each CORSIKA shower is used five times,
placing randomly the shower footprint on the area of the SD array.

###  4.2 Photon energy scale

To account for the different energy scale of photon showers (see Section 2.2),
we retain the standard reconstruction of the shower direction and energy
estimator $S(1000)$ and replace $E_{\rm had}$ with a new function of $S(1000)$
and $\theta$ calibrated with photon simulations.

Taking advantage of the large statistics of simulated events, a look-up table
is built with the mean value of the logarithm of the Monte Carlo true energy
of photons, $L_{E}={\log_{10}(E/{\rm eV}})$, in bins of $S(1000)$ and
$\theta$. Only non-preshowering simulated events triggering the SD array are
used and the photon simulations are weighted to a reference spectrum $\propto
E^{-2}$. The table serves as the desired function of $S(1000)$ and $\theta$:
for each event the assigned photon energy $E_{\gamma}$ is $10^{L_{E}}$ eV,
where $L_{E}$ is the tabulated value for the bin to which the event belongs.

###  4.3 Event selection

To ensure events with a good reconstruction and optimize the photon search, we
select only events in which:

  * •

the detector with the highest signal is surrounded by a hexagon of six
stations that are fully operational,

  * •

the reconstructed zenith is in the range $30^{\circ}-60^{\circ}$,

  * •

the reconstructed energy in hadronic scale is required to be $E_{\rm
had}>10^{18.5}$ eV.

The second condition ensures that the majority of selected photon-induced
showers reach their maximum development before being detected, as $\langle
X_{\rm max}\rangle$ exceeds the vertical atmospheric depth at the Observatory
site already at $E_{\gamma}>10^{19}$ eV and increases with increasing
$E_{\gamma}$. The latter condition reduces the background from low energy
showers initiated by nuclei while not impacting the selection efficiency for
photons above 10 EeV when combined with the observable-related cuts.

The timing and shape of the signal acquired by each PMT of the WCDs is of
paramount importance for the risetime observable $\Delta$ defined in Eq. 3.2.
A special attention is given to the quality of the PMT outputs to avoid
distorted signals from features shown by single PMTs, which may significantly
affect the discriminant observables and could result in mimicking photon-like
signals. Besides the quality criteria applied routinely to the SD events,
three new filters have been designed specifically for the photon searches,
rejecting PMTs with significant afterpulses, an oscillating baseline pattern
or a nonphysically slow decrease of their signal. Since only the individual
PMTs are removed from the events, the effect on the number of events selected
is negligible.

With the selection criteria described above, the set of air showers detected
with the SD in the time period of interest consists of 105,064 events.

To guarantee a good quality of the observables defined in Section 3,
additional criteria are applied:

  * •

only WCDs with non-saturated signal $>6$ VEM and radial distance in the range
$600-2000$ m are considered in the calculation of $\Delta$, the event being
selected when the resulting number of selected stations $N\geq 4$,

  * •

a minimum of one WCD above radial distance $1000\text{\,}\mathrm{m}$ is
required to compute $L_{\rm LDF}$.

The search is restricted to $E_{\gamma}>10^{19}$ eV, corresponding to the
energy at which photons are detected by the SD with a trigger efficiency close
to 100%.

These criteria reduce the set of selected data to 48,947 events.

The resolution on $E_{\gamma}$ for simulated photon events fulfilling the same
selection criteria is almost constant over the energy range considered, being
$\simeq 30\%$ with no bias for non-preshowering photons from a $E^{-2}$
spectrum. Showers undergoing preshowering in the geomagnetic field, on the
other hand, are characterized by a shallower depth of the shower maximum with
respect to non-preshowering ones of the same primary energy. This results in
an underestimate of their energy of 30% using the look-up table constructed as
described in Section 4.2.

###  4.4 Selection of photon candidates

To combine the information contained in the two discriminating variables and
define a criterion for the identification of photon candidate events, a burn
sample is extracted from the set in a way that guarantees a time distribution
of events representing a fair sample of the instantaneous exposure over the
time period considered. The burn sample consists of 886 events, corresponding
to $\sim$1.8% of the total selected events. The use of such a subset avoids
reliance on simulations of showers initiated by nuclei, which constitute the
background for the photon search, and the related uncertainties stemming from
the assumptions on the mass composition and the modeling of hadronic
interactions.

Figure 3: Distribution of $\Delta$ (left) and $L_{\rm LDF}$ (right) as a
function of the photon energy $E_{\gamma}$ for the burn sample and photon
simulations. The bands represent one standard deviation of the photon
distributions.

For photon showers, the distribution of the SD observables $\Delta$ and
$L_{\rm LDF}$ described in Section 3 is dependent on the energy and zenith
angle of the primary particle (see, e.g., Fig. 3). To define a single
selection criterion for photon candidates as independent as possible from
direction and energy, the mean and standard deviation of the distribution of
each initial variable are computed, for non-preshowering photons, in 30 bins
of roughly equal statistics in the ($S(1000)$, $\theta$) space (five in
$S(1000)$ and six in $\theta$). The reference spectrum is $\propto E^{-2}$
hence simulations are weighted accordingly. The variables $\tilde{\Delta}$ and
$\tilde{L}_{\rm LDF}$ are then defined as linear transformations of the
initial ones centered around $0$ and expressed in units of the standard
deviation of the corresponding distributions of non-preshowering photons:

$\displaystyle\tilde{\Delta}$ $\displaystyle=$
$\displaystyle\frac{\Delta-\langle\Delta\rangle^{i}}{\sigma_{\Delta}^{i}},$
(4.1) $\displaystyle\tilde{L}_{\rm LDF}$ $\displaystyle=$
$\displaystyle\frac{L_{\rm LDF}-\langle L_{\rm LDF}\rangle^{i}}{\sigma_{L_{\rm
LDF}}^{i}},$ (4.2)

where $i$ is the index of the bin corresponding to the specific event. The
distributions of $\tilde{\Delta}$ and $\tilde{L}_{\rm LDF}$ are shown in Fig.
4.

Figure 4: Distributions of the variables $\tilde{L}_{\rm LDF}$ and
$\tilde{\Delta}$ of the burn sample (points) and non-preshowering photons
weighted to an $E^{-2}$ spectrum (contours). The contour levels encompass
respectively 10-30-50-70-90% of the distribution. The threshold photon energy
is $10^{19}\leavevmode\nobreak\ $eV.

The two variables are combined using a Fisher discriminant analysis with the
burn sample representing the background and photon simulations the signal. The
transformation is normalized as to be equivalent to a rotation in the
$(\tilde{\Delta},\tilde{L}_{\rm LDF})$ plane. The resulting axis is shown in
Fig. 4.

The photon candidate cut is chosen a priori as the median of the photon sample
of non-preshowering events weighted to a $E^{-2}$ spectrum. This cut value
constitutes a good compromise between efficiency and purity. Any event falling
above this cut, shown as a dashed line in Fig. 4, will be considered as a
photon candidate.

##  5 Results of the photon search

Excluding the burn sample from the final analysis, the search sample consists
of 48,061 selected events. Application of the photon search method yields the
summary plots shown in Fig. 5 for $E_{\gamma}\geq 10^{19}\leavevmode\nobreak\
$eV. Analyzing the data in the $(\tilde{\Delta}$, $\tilde{L}_{\rm LDF})$ plane
results in the red points displayed in the left panel, on which are drawn the
same contour levels as in Fig. 4 of the distribution for photons as well as
the Fisher axis and the candidate-cut Fisher value. In the right panel, the
corresponding distributions of the Fisher discriminant value are shown as
normalized histograms for the burn sample, the search sample, as well as the
simulated photon sample separated in non-preshowering and preshowering. For
reference, the candidate cut is shown as the vertical line, while the result
of an exponential fit to the 5% of events from the burn sample with the
largest Fisher values is drawn to guide the eye in the interpretation of the
tail of the Fisher distribution of the search sample.

We find 16 (1) [0] photon candidates above $10^{19}\leavevmode\nobreak\ $eV
($2{\times}10^{19}\leavevmode\nobreak\ $eV)
[$4{\times}10^{19}\leavevmode\nobreak\ $eV]. The number of observed candidates
is in statistical agreement with what is expected from the exponential fit to
the burn sample, with a difference of -0.3 standard deviations. In addition,
no peak-like features above the selection cuts that would indicate the
presence of a photon population are observed above the fall-off of the
distribution. Overall, therefore, the Fisher distribution of photon candidates
is consistent with the expectations of a background of UHECR events.

To search for further imprints that would be indicative of the presence of
photon events, we have checked that no candidates are coincident in time. We
have also searched for small-scale clustering in arrival directions that would
be indicative of repeaters and thus of point-like source of photons. No such
clustering is observed, and the arrival directions of the candidates are
distributed in accordance with the directional exposure of the cosmic-ray
background events.

Figure 5: (Left) Distribution in the $(\tilde{\Delta}$,$\tilde{L}_{\rm LDF})$
plane of the events of the search sample compared with the $E^{-2}$-weighted
distribution for non-preshowering photons (same contour levels as in Fig. 4).
(Right) Distribution of the Fisher discriminant value for the burn sample
(grey), search sample (red) and $E^{-2}$-weighted photon simulations (non-
preshowering in blue and preshowering in light blue). The search sample and
photon distributions are scaled as to have the same integral as the burn
sample one. The vertical line stands for the candidate cut. The dashed line
shows the result of the exponential fit to the 5% of events in the burn sample
with larger Fisher values.

From the absence of photon flux measurements, upper limits are derived from
the number of candidates above a minimum energy $E_{0}$. The signal efficiency
of the analysis, $\epsilon_{\gamma}$, is estimated by applying the quality
cuts and the selection of candidates to a Monte-Carlo set of photon air
showers using the same procedure as for the data. The (weighted) ratio of the
number of selected candidates (preshowering or not) reconstructed in the
selected zenithal range and with assigned photon energy above the threshold
$E_{0}$ to the total number of simulated events with generated zenith and
energy in the desired ranges is computed. The weight is assigned so as to
obtain, from the generated $E^{-1}$ spectrum, the results corresponding to the
spectral shape assumed for the photon flux $\Phi(E)$ we want to bind, as
explained in Section 4.1.

$\epsilon_{\gamma}(>E_{0})=\frac{\int_{\Delta\Omega}\int_{E_{\gamma}>E_{0}}\Phi_{\rm
obs}(E_{\gamma},\Omega_{\rm rec})\,d\Omega_{\rm
rec}\,dE_{\gamma}}{\int_{\Delta\Omega}\int_{E>E_{0}}\Phi_{\rm
gen}(E)\,w(E)\,d\Omega\,dE}$ (5.1)

where $\Delta\Omega$ is the solid angle corresponding to the analysis
selection,

$\Phi_{\rm obs}(E_{\gamma},\Omega_{\rm rec})=\int_{\Delta\Omega_{\rm
gen}}\int_{\Delta E_{\rm gen}}p^{\rm cand}(E_{\gamma},\Omega_{\rm
rec}|E,\Omega)\,\Phi_{\rm gen}(E)\,w(E)\,d\Omega\,dE$ (5.2)

$p^{\rm cand}(E_{\gamma},\Omega_{\rm rec}|E,\Omega)$ represents the
probability that a shower of true Monte Carlo energy $E$ generated in the
solid angle $\Omega$ is selected as candidate photon event with photon energy
$E_{\gamma}$ and in the reconstructed solid angle $\Omega_{\rm rec}$. Here the
integral is extended to the whole range of simulations in solid angle
$\Delta\Omega_{\rm gen}$ and energy range $\Delta E_{\rm gen}$. In other
words, the efficiency represents the ratio between the integral spectrum
folded into the detector response and the true underlying one, accounting for
any effects due to inaccuracies in the photon energy and direction assignment
(and possible biases for model fluxes different from $E^{-2}$).

Upper limits to the integrated flux at CL confidence level are derived using
the Feldman-Cousins method [56]:

$\mathcal{F}^{\text{CL}}(E_{\gamma}\geq E_{0})=\frac{N_{\gamma}^{\rm
CL}}{A(1-\eta_{B})\epsilon_{\gamma}},$ (5.3)

with $A=46,900\leavevmode\nobreak\ $km2 sr yr the total exposure in the time
period 1 January 2004 to 30 June 2020 and in the zenithal range considered,
$\eta_{B}$ the exposure fraction corresponding to the burn sample of data used
to set the candidate cut and excluded from the search, and
$N_{\gamma}^{\text{CL}}$ the Feldman-Cousins upper limit to the number of
photon events computed at confidence level CL in the conservative hypothesis
of no background event expected.

While the reconstruction and candidate selection procedure is independent from
the model spectrum, upper limits depend on the spectrum through the
efficiency. The values obtained for both efficiencies and upper limits at 95%
CL are listed in Table 1 in the case of a power law spectrum $E^{-\alpha}$ for
three different values of the exponent.

$E_{0}$ | 10 EeV | 20 EeV | 40 EeV | 10 EeV | 20 EeV | 40 EeV
---|---|---|---|---|---|---
$\alpha$ | Efficiency ($\epsilon_{\gamma}$) | Flux U.L. ($10^{-3}\,{\rm km}^{-2}\,{\rm sr}^{-1}\,{\rm yr}^{-1}$)
1.0 | 0.37 | 0.41 | 0.42 | 1.48 | 0.273 | 0.162
2.0 | 0.26 | 0.36 | 0.39 | 2.11 | 0.312 | 0.172
3.0 | 0.19 | 0.34 | 0.41 | 2.86 | 0.332 | 0.166
Table 1: Efficiencies and upper limits to the integral flux of photons above
$E_{0}$ for different values of the exponent $\alpha$ in the case of a power-
law spectrum. The limits (95% CL) are listed on the right.

The 95% CL upper limits to the integral photon flux with a spectral index
$\alpha=2$ are shown in Fig. 6, together with the results of other photon
searches using other components of the Observatory for lower energy thresholds
[37] and those reported from Telescope Array [57]. In the energy region above
$10^{19}\leavevmode\nobreak\ $eV, the limits obtained in this study are the
most stringent currently available in literature. Their astrophysical
implications are briefly discussed in Section 6.

Figure 6: Upper limits (at 95 % CL) on the integral photon flux above
$10^{19}$ eV determined in this study (red squares). Shown are also previous
upper limits reported in [37] (Auger Hybrid, blue circles), and Telescope
Array [57] (grey triangles). The grey band stands for the range of expected
photon fluxes under the assumption of a pure-proton scenario [58], while the
green one for the case of a mixed composition [59]. In addition, several
expected photon fluxes from the decay of SHDM particles are shown (lines). The
dashed violet [60] and dot-dashed orange lines [29] have been obtained through
private communication with one of the authors.

##  6 Discussion

As mentioned in Section 1, the interactions of UHECRs with the background
photon fields permeating the Universe, most notably the cosmic microwave
background, guarantee the existence of a cosmogenic diffuse flux of UHE
photons. However, the short photon horizon compared to the cosmic-ray one
makes the photon flux lower than the UHECR one by several orders of magnitude.
Moreover, the exact photon flux is further reduced as the mass composition of
UHECRs gets heavier, because the photons produced from primary heavy nuclei
are of lower energies than those from lighter ones. This is evidenced in Fig.
6 through the grey band standing for the photon flux expected from proton
interactions [58], which is observed one order of magnitude (as well as to
extend to higher energies) above that expected from a mixed composition shown
in green [59]. The size of the bands accounts for several hypotheses necessary
for the modeling of the UHECR flux, such as the maximum acceleration energy of
the nuclei at the sources and the shape of the energy spectrum of the
accelerated particles. Note that the mixed-composition case is designed to
reproduce the mass composition and the energy spectrum inferred from the data
collected at the Observatory. The sensitivity of the present search around
$2\times 10^{19}$ eV approaches the most optimistic expectations of the
cosmogenic photon flux from protons, while it is one order of magnitude above
those from the mixed-composition model. In the future, it will thus be
possible to probe the proton expectations with improved sensitivity to
photons, and thus to confirm indirectly the presence of a mixed composition at
the highest energies.

Our upper limits have also consequences for models of super-heavy dark matter
(SHDM). Although the absence of a detectable flux of UHE photons at a level
comparable to that of cosmic rays excludes that SHDM particles could explain
the bulk of UHECRs, they can still contribute in a subdominant way to the all-
particle spectrum. For a $E^{-2}$ photon spectrum, and using the flux of
cosmic rays measured using the Pierre Auger Observatory [46], the upper limits
(at 95% CL) to the integral photon fraction are: 1.6%, 1.2% and 3.2% above 10
EeV, 20 EeV and 40 EeV respectively.

There are several theoretical motives to search for SHDM related to particle
physics in the early Universe. A recent exploration of the constraints
provided by our photon-flux limits on some models is given elsewhere [31].
Here, we restrict ourselves to illustrate in Fig. 6 the discovery potential
with searches for UHE photons: we show as the dashed violet line and the dot-
dashed red line the expected photon fluxes in the case of hadronic [60] and
leptonic [29] decay channels, for a mass of the SHDM particles of
$10^{10}\leavevmode\nobreak\ $GeV and a lifetime of
$3{\times}10^{21}\leavevmode\nobreak\ $yr. In addition, we also show as the
black dot-dashed line an allowed flux extending to higher energy for a mass of
$10^{12}\leavevmode\nobreak\ $GeV and a lifetime of
$10^{23}\leavevmode\nobreak\ $yr (hadronic decay channel). As the sensitivity
of current photon searches increases, it will be possible to further constrain
the allowed values of mass and lifetime.

In summary, a search for photons using the surface detector of the Pierre
Auger Observatory has been performed, leading to the most stringent limits on
the diffuse photon flux above $10^{19}$ eV and complementing the limits
obtained in the lower energy decade. The current limits are not able to
challenge the flux of cosmogenic photons expected for the mixed-composition
model but start to probe the most optimistic predictions of the pure-proton
scenarios. Future data will enable us to constrain further the different
mechanisms expected to produce UHE photons. In particular, the completion of
the upgraded Observatory [61] will allow us to increase the sensitivity of the
different discriminating analyses. Finally, the analysis presented in this
study can also serve for directional searches from specific targets or
searches in time coincidence with observations from other messengers.

## Acknowledgments

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); 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 – Grant No. MSMT CR LTT18004,
LM2015038, LM2018102, 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;
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’Istruzione, dell’Universitá e della Ricerca (MIUR); CETEMPS Center of
Excellence; Ministero degli Affari Esteri (MAE); 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, grant No. DIR/WK/2018/11; National Science
Centre, Grants No. 2016/22/M/ST9/00198, 2016/23/B/ST9/01635, and
2020/39/B/ST9/01398; 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/CCCDI UEFISCDI, grant no. PN19150201/16N/2019 and
PN1906010 within the National Nucleus Program, and projects number TE128, PN-
III-P1-1.1-TE-2021-0924/TE57/2022 and PED289, 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

  * [1] P. S. Coppi and F. A. Aharonian, Constraints on the VHE emissivity of the universe from the diffuse GeV gamma-ray background, Astrophys. J. Lett. 487 (1997) L9–L12, [astro-ph/9610176].
  * [2] C. Heiter, D. Kuempel, D. Walz, and M. Erdmann, Production and propagation of ultra-high energy photons using CRPropa 3, Astroparticle Physics 102 (2018) 39–50, [astro-ph/1710.11406].
  * [3] Pierre Auger Collaboration, A. Aab et al., Depth of maximum of air-shower profiles at the Pierre Auger Observatory. I. Measurements at energies above $10^{17.8}$ eV, Phys.Rev. D90 (2014), no. 12 122005, [arXiv:1409.4809].
  * [4] Pierre Auger Collaboration, A. Aab et al., Depth of maximum of air-shower profiles at the Pierre Auger Observatory. II. Composition implications, Phys. Rev. D 90 (2014), no. 12 122006, [arXiv:1409.4809].
  * [5] Pierre Auger Collaboration, A. Aab et al., Evidence for a mixed mass composition at the ‘ankle’ in the cosmic-ray spectrum, Phys. Lett. B762 (2016) 288–295, [arXiv:1609.08567].
  * [6] Pierre Auger Collaboration, A. Aab et al., Inferences on mass composition and tests of hadronic interactions from 0.3 to 100 EeV using the water-Cherenkov detectors of the Pierre Auger Observatory, Phys. Rev. D96 (2017), no. 12 122003, [arXiv:1710.07249].
  * [7] Fermi-LAT Collaboration, A. A. Abdo et al., The Spectrum of the Isotropic Diffuse Gamma-Ray Emission Derived From First-Year Fermi Large Area Telescope Data, Phys. Rev. Lett. 104 (2010) 101101, [arXiv:1002.3603].
  * [8] IceCube Collaboration, M. G. Aartsen et al., Constraints on Ultrahigh-Energy Cosmic-Ray Sources from a Search for Neutrinos above 10 PeV with IceCube, Phys. Rev. Lett. 117 (2016), no. 24 241101, [arXiv:1607.05886]. [Erratum: Phys. Rev. Lett.119,no.25,259902(2017)].
  * [9] Pierre Auger Collaboration, A. Aab et al., Probing the origin of ultra-high-energy cosmic rays with neutrinos in the EeV energy range using the Pierre Auger Observatory, JCAP 1910 (2019), no. 10 022, [arXiv:1906.07422].
  * [10] V. Berezinsky, A. Gazizov, and O. Kalashev, Cascade photons as test of protons in UHECR, Astropart. Phys. 84 (2016) 52–61, [arXiv:1606.09293].
  * [11] J. Heinze, D. Boncioli, M. Bustamante, and W. Winter, Cosmogenic Neutrinos Challenge the Cosmic Ray Proton Dip Model, Astrophys. J. 825 (2016), no. 2 122, [arXiv:1512.05988].
  * [12] A. D. Supanitsky, Implications of gamma-ray observations on proton models of ultrahigh energy cosmic rays, Phys. Rev. D94 (2016), no. 6 063002, [arXiv:1607.00290].
  * [13] R. Alves Batista, R. M. de Almeida, B. Lago, and K. Kotera, Cosmogenic photon and neutrino fluxes in the Auger era, JCAP 1901 (2019) 002, [arXiv:1806.10879].
  * [14] M. S. Muzio, M. Unger, and G. R. Farrar, Progress towards characterizing ultrahigh energy cosmic ray sources, Phys. Rev. D100 (2019), no. 10 103008, [arXiv:1906.06233].
  * [15] S. Coleman and S. L. Glashow, Cosmic ray and neutrino tests of special relativity, Physics Letters B 405 (1997), no. 3 249 – 252.
  * [16] P. Bhattacharjee and G. Sigl, Origin and propagation of extremely high-energy cosmic rays, Physics Reports 327 (2000), no. 3 109 – 247.
  * [17] W. Bietenholz, Cosmic rays and the search for a Lorentz Invariance Violation, Physics Reports 505 (2011), no. 5 145 – 185.
  * [18] P. W. Gorham, A. Connolly, P. Allison, J. J. Beatty, K. Belov, D. Z. Besson, W. R. Binns, P. Chen, J. M. Clem, S. Hoover, M. H. Israel, J. Nam, D. Saltzberg, G. S. Varner, and A. G. Vieregg, Implications of ultrahigh energy neutrino flux constraints for Lorentz-invariance violating cosmogenic neutrinos, Phys. Rev. D 86 (Nov, 2012) 103006.
  * [19] G. Rubtsov, P. Satunin, and S. Sibiryakov, Prospective constraints on Lorentz violation from ultrahigh-energy photon detection, Phys. Rev. D 89 (Jun, 2014) 123011.
  * [20] J. D. Tasson, What do we know about Lorentz invariance?, Reports on Progress in Physics 77 (2014), no. 6 062901.
  * [21] G. Rubtsov, P. Satunin, and S. Sibiryakov, Constraints on violation of Lorentz invariance from atmospheric showers initiated by multi-TeV photons, Journal of Cosmology and Astroparticle Physics 2017 (2017), no. 05 049.
  * [22] A. Nicolaidis, Axionic shortcuts for high energy photons, Journal of Cosmology and Astroparticle Physics 2010 (2010), no. 04 013.
  * [23] V. Berezinsky, M. Kachelriess, and A. Vilenkin, Ultrahigh-energy cosmic rays without GZK cutoff, Phys. Rev. Lett. 79 (1997) 4302–4305, [astro-ph/9708217].
  * [24] D. J. Chung, E. W. Kolb, and A. Riotto, Superheavy dark matter, Phys. Rev. D 59 (1998) 023501, [hep-ph/9802238].
  * [25] M. Garny, M. Sandora, and M. S. Sloth, Planckian Interacting Massive Particles as Dark Matter, Phys. Rev. Lett. 116 (2016), no. 10 101302, [arXiv:1511.03278].
  * [26] V. Kuzmin and V. Rubakov, Ultrahigh-energy cosmic rays: A Window to postinflationary reheating epoch of the universe?, Phys. Atom. Nucl. 61 (1998) 1028, [astro-ph/9709187].
  * [27] R. Aloisio, S. Matarrese, and A. V. Olinto, Super Heavy Dark Matter in light of BICEP2, Planck and Ultra High Energy Cosmic Rays Observations, JCAP 1508 (2015), no. 08 024, [arXiv:1504.01319].
  * [28] E. Alcantara, L. A. Anchordoqui, and J. F. Soriano, Hunting for superheavy dark matter with the highest-energy cosmic rays, Phys. Rev. D99 (2019), no. 10 103016, [arXiv:1903.05429].
  * [29] M. Kachelriess, O. E. Kalashev, and M. Yu. Kuznetsov, Heavy decaying dark matter and IceCube high energy neutrinos, Phys. Rev. D98 (2018), no. 8 083016, [arXiv:1805.04500].
  * [30] K. Ishiwata, O. Macias, S. Ando, and M. Arimoto, Probing heavy dark matter decays with multi-messenger astrophysical data, JCAP 2001 (2020), no. 01 003, [arXiv:1907.11671].
  * [31] Pierre Auger Collaboration, P. Abreu et al., Limits to gauge coupling in the dark sector set by the non-observation of instanton-induced decay of Super-Heavy Dark Matter in the Pierre Auger Observatory data, arXiv:2203.08854.
  * [32] Pierre Auger Collaboration, A. Aab et al., The Pierre Auger Cosmic Ray Observatory, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 798 (2015) 172 – 213.
  * [33] Pierre Auger Collaboration, J. Abraham et al., An upper limit to the photon fraction in cosmic rays above $10^{19}\leavevmode\nobreak\ $eV from the Pierre Auger Observatory, Astroparticle Physics 27 (2007), no. 2 155 – 168\.
  * [34] Pierre Auger Collaboration, J. Abraham et al., Upper limit on the cosmic-ray photon flux above $10^{19}\leavevmode\nobreak\ $eV using the surface detector of the Pierre Auger Observatory, Astroparticle Physics 29 (2008), no. 4 243 – 256.
  * [35] Pierre Auger Collaboration, C. Bleve, Update of the neutrino and photon limits from the Pierre Auger Observatory, Proc. 34th ICRC, The Hague (2015).
  * [36] Pierre Auger Collaboration, J. Rautenberg, Limits on ultra-high energy photons with the Pierre Auger Observatory, PoS ICRC2019 (2019) 398.
  * [37] P. Savina and Pierre Auger Collaboration, A search for ultra-high-energy photons at the Pierre Auger Observatory exploiting air-shower universality, Proc. 37th International Cosmic Ray Conference (Berlin, Germany) PoS (ICRC2021) (2021) 373.
  * [38] Pierre Auger Collaboration, P. Abreu et al., A search for photons with energies above $2\times 10^{17}$ eV using hybrid data from the low-energy extensions of the Pierre Auger Observatory, accepted for publication in Astrophys. J. (2022) [arXiv:2205.14864].
  * [39] R. Engel, D. Heck, and T. Pierog, Extensive air showers and hadronic interactions at high energy, Annual review of nuclear and particle science 61 (2011) 467–489.
  * [40] M. Risse and P. Homola, Search for ultra-high energy photons using air showers, Modern Physics Letters A 22 (2007), no. 11 749–766, [0702632].
  * [41] Pierre Auger Collaboration, A. Aab et al., Reconstruction of events recorded with the surface detector of the pierre auger observatory, Journal of Instrumentation 15 (oct, 2020) P10021–P10021.
  * [42] L. Landau and I. Pomeranchuk, Limits of applicability of the theory of bremsstrahlung electrons and pair production at high-energies, Dokl. Akad. Nauk Ser. Fiz. 92 (1953) 535–536.
  * [43] A. Migdal, Bremsstrahlung and pair production in condensed media at high-energies, Phys. Rev. 103 (1956) 1811–1820.
  * [44] B. McBreen and C. J. Lambert, Interactions of High-energy ($E>5\times 10^{19}$-eV) Photons in the Earth’s Magnetic Field, Phys. Rev. D24 (1981) 2536–2538.
  * [45] Gaisser, Thomas K and Hillas, A Michael, Reliability of the method of constant intensity cuts for reconstructing the average development of vertical showers, .
  * [46] Pierre Auger Collaboration, A. Aab et al., Measurement of the cosmic-ray energy spectrum above $2.5{\times}10^{18}$ eV using the Pierre Auger Observatory, Phys. Rev. D 102 (2020), no. 6 062005, [arXiv:2008.06486].
  * [47] Pierre Auger Collaboration, A. Aab et al., Inferences on mass composition and tests of hadronic interactions from 0.3 to 100 EeV using the water-Cherenkov detectors of the Pierre Auger Observatory, Phys. Rev. D 96 (2017), no. 12 122003, [arXiv:1710.07249].
  * [48] D. Heck, J. Knapp, J. N. Capdevielle, G. Schatz, and T. Thouw, CORSIKA: a Monte Carlo code to simulate extensive air showers. 1998\.
  * [49] 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 (Sep, 2015) 034906.
  * [50] “Fluka website.” https://fluka.cern.
  * [51] C. Ahdida et al., New Capabilities of the FLUKA Multi-Purpose Code, Frontiers in Physics (2022), no. 9 788253.
  * [52] G. Battistoni, T. Boehlen, F. Cerutti, P. Chin, L. Esposito, A. Fassò, A. Ferrari, A. Lechner, A. Empl, A. Mairani, A. Mereghetti, P. Garcia Ortega, J. Ranft, S. Roesler, P. Sala, V. Vlachoudis, and G. Smirnov, Overview of the FLUKA code, Annals of Nuclear Energy (2015), no. 82 10.
  * [53] P. Homola, M. Risse, R. Engel, D. Gora, J. Pekala, B. Wilczynska, and H. Wilczynski, Characteristics of geomagnetic cascading of ultrahigh energy photons at the southern and northern sites of the Pierre Auger Observatory, Astropart. Phys. 27 (2007) 174–184, [astro-ph/0608101].
  * [54] P. Homola, R. Engel, A. Pysz, and H. Wilczyński, Simulation of ultra-high energy photon propagation with PRESHOWER 2.0, Computer Physics Communications 184 (2013) 1468.
  * [55] S. Argirò et al., The Offline Software Framework of the Pierre Auger Observatory, Nuclear Instruments and Methods in Physics Research Section A 580 (2007) 1485.
  * [56] G. J. Feldman and R. D. Cousins, Unified approach to the classical statistical analysis of small signals, Physical Review D 57 (Apr, 1998) 3873–3889.
  * [57] Telescope Array Collaboration, R. U. Abbasi et al., Constraints on the diffuse photon flux with energies above $10^{18}$ eV using the surface detector of the Telescope Array experiment, Astropart. Phys. 110 (2019) 8–14, [arXiv:1811.03920].
  * [58] B. Sarkar, K.-H. Kampert, and J. Kulbartz, Ultra-High Energy Photon and Neutrino Fluxes in Realistic Astrophysical Scenarios, in Proceedings, 32nd International Cosmic Ray Conference (ICRC 2011): Beijing, China, August 11-18, 2011, vol. 2, p. 198, 2011.
  * [59] A. Bobrikova, M. Niechciol, M. Risse, and P. Ruehl, Predicting the UHE photon flux from GZK-interactions of hadronic cosmic rays using CRPropa 3, Proc. 37th International Cosmic Ray Conference (Berlin, Germany) PoS (ICRC2021) (2021) 449.
  * [60] O. K. Kalashev and M. Yu. Kuznetsov, Constraining heavy decaying dark matter with the high energy gamma-ray limits, Phys. Rev. D94 (2016), no. 6 063535, [arXiv:1606.07354].
  * [61] Pierre Auger Collaboration, A. Aab et al., The Pierre Auger Observatory Upgrade - Preliminary Design Report, arXiv:1604.03637.