# Diffuse flux of ultra-high energy photons from cosmic-ray interactions in
the disk of the Galaxy and implications for the search for decaying super-
heavy dark matter

Corinne Bérat Laboratoire de Physique Subatomique & Cosmologie, CNRS/IN2P3,
Université Grenoble Alpes, Grenoble, France  Carla Bleve Laboratoire de
Physique Subatomique & Cosmologie, CNRS/IN2P3, Université Grenoble Alpes,
Grenoble, France  Olivier Deligny Laboratoire de Physique des 2 Infinis Irène
Joliot-Curie, CNRS/IN2P3, Université Paris-Saclay, Orsay, France  François
Montanet Laboratoire de Physique Subatomique & Cosmologie, CNRS/IN2P3,
Université Grenoble Alpes, Grenoble, France  Pierpaolo Savina Laboratoire de
Physique des 2 Infinis Irène Joliot-Curie, CNRS/IN2P3, Université Paris-
Saclay, Orsay, France  University of Wisconsin-Madison, Department of Physics
and WIPAC, Madison, WI, U.S.A  Zoé Torrès Laboratoire de Physique Subatomique
& Cosmologie, CNRS/IN2P3, Université Grenoble Alpes, Grenoble, France

###### Abstract

An estimate of the expected photon flux above $10^{17}~{}$eV from the
interactions of ultra-high energy cosmic rays with the matter in the Galactic
disk is presented. Uncertainties arising from the distribution of the gas in
the disk, the absolute level of the cosmic ray flux, and the composition of
the cosmic rays are taken into account. Within these uncertainties, the
integrated photon flux above $10^{17}~{}$eV is, averaged out over Galactic
latitude less than $5^{\circ}$, between $\simeq
3.2{\times}10^{-2}~{}$km${}^{-2}~{}$yr${}^{-1}~{}$sr-1 and $\simeq
8.7{\times}10^{-2}~{}$km${}^{-2}~{}$yr${}^{-1}~{}$sr-1. The all-sky average
value amounts to $\simeq
1.1{\times}10^{-2}~{}$km${}^{-2}~{}$yr${}^{-1}~{}$sr-1 above $10^{17}~{}$eV
and decreases roughly as $E^{-2}$, making this diffuse flux the dominant one
from cosmic-ray interactions for energy thresholds between $10^{17}~{}$eV and
$10^{18}$ eV. Compared to the current sensitivities of detection techniques, a
gain between two and three orders of magnitude in exposure is required for a
detection below $\simeq 10^{18}$ eV. The implications for searches for photon
fluxes from the Galactic center that would be indicative of the decay of
super-heavy dark matter particles are discussed, as the photon flux presented
in this study can be considered as a floor below which other signals would be
overwhelmed.

astroparticle physics — cosmic rays — radiation mechanisms: non-thermal

††journal: ApJ

##  1 Introduction

Cosmic-ray electrons and protons in the $1$–$100$ GeV energy range can produce
gamma rays subsequent to their interactions with the interstellar matter in
the Galactic disk via high-energy electron Bremsstrahlung, nucleon-nucleon
processes, and, to a lesser extent, inverse-Compton interactions with low-
energy photons. The diffuse flux of ${>}100~{}$MeV and GeV gamma rays observed
in a narrow band along the Galactic plane (Narayanan et al., 2008; Fichtel et
al., 1972, 1975; Mayer-Hasselwander et al., 1980, 1982; Hunter et al., 1997;
Ackermann et al., 2012) is consistent with these expectations (Fichtel &
Kniffen, 1984; Bloemen et al., 1984; Bloemen, 1989). At higher energies, the
flux of cosmic-ray electrons is overwhelmed by that of protons and heavier
nuclei, so the dominant photon-producing process is expected to be the
nucleon-nucleon collision. The Milagro experiment has allowed the observation
of a diffuse flux at TeV energies from the Cygnus region in the Galactic disk
(Atkins et al., 2005), while the ARGO-YBJ observatory has allowed that of a
diffuse flux extended between $25^{\circ}$ and $100^{\circ}$ in Galactic
longitude for $|b|<5^{\circ}$ in Galactic latitude (Bartoli et al., 2015).
More recently, a similar diffuse flux at sub-PeV energies has been reported
from data collected at the Tibet-AS$\gamma$ observatory (Amenomori et al.,
2021). These observations are indeed in agreement with expectations of photon
emission from matter irradiated by cosmic rays, see e.g. Berezinsky et al.
(1993); Chardonnet et al. (1995); Evoli et al. (2008); Carceller & Masip
(2017); Lipari & Vernetto (2018). They could shed light on the spectra of
cosmic rays in distant regions of the Galaxy (Lipari & Vernetto, 2018).

Above $10^{17}~{}$eV, several campaigns are underway to search for ultra-high
energy (UHE) photons. On the one hand, the detection of a diffuse flux would
allow probing the expectations related to $\pi^{0}$ decays following the
interactions of ultra-high energy cosmic rays (UHECRs) with the background
photon fields or dust permeating the source environments or the extragalactic
space. On the other hand, the detection of localized fluxes would allow the
discovery of cosmic ray sources or, if they are in the direction of the
Galactic center, to highlight the presence of decaying super-heavy dark matter
produced in the early Universe. However, the interstellar matter of the
Galactic disk is irradiated by the UHECRs. Therefore, a diffuse flux of UHE
photons is expected to be produced by the same mechanism responsible for the
emission in the GeV–PeV range. An estimate of this flux, in addition to that
mentioned above originating from the interaction of UHECRs with the photon
fields permeating the Universe, is needed to know the “background” that could
hide the emission from sources in the Galactic disk or the Galactic center.

Because a few percent only of UHECRs interact at least once with the matter in
the Galactic disk, and because the Galaxy is essentially transparent to
photons with energies in excess of $10^{17}~{}$eV (see Appendix A),
estimations of the diffuse flux (per steradian) of UHE photons at energy $E$,
$\phi_{\gamma}(E,\mathbf{n})$, can be made in a straightforward manner by
integrating the position-dependent emission rate per unit volume and unit
energy along the line of sight,

$\phi_{\gamma}(E,\mathbf{n})=\frac{1}{4\pi}\int_{0}^{\infty}\mathrm{d}s~{}q_{\gamma}(E,\mathbf{x}_{\sun}+s\mathbf{n}).$
(1)

Here, $\mathbf{x}_{\sun}$ is the position of the Solar system in the Galaxy
and $\mathbf{n}\equiv\mathbf{n}(\ell,b)$ is a unit vector on the sphere
pointing to the longitude $\ell$ and latitude $b$, in Galactic coordinates.
The high degree of isotropy of UHECRs above $10^{17}~{}$eV, as observed on
Earth (Aab et al., 2020a), leads to an equally isotropic irradiation of gas
distribution in the Galactic disk, and therefore an isotropic emission that
explains the $1/4\pi$ factor. The emission of UHE photons stems from the
creation and decay of unstable mesons, mostly $\pi^{0}$ but also, to a lower
extent, $\eta$ and $\eta^{\prime}$, in the inelastic interactions of UHECRs
with interstellar gas,

$q_{\gamma}(E,\mathbf{x})=4\pi\sum_{i,j}n_{j}(\mathbf{x})\int_{E}^{\infty}\mathrm{d}E^{\prime}\phi_{i}(E^{\prime})\sigma_{ij}(E^{\prime})\frac{\mathrm{d}N_{ij}^{\gamma}}{\mathrm{d}E}(E^{\prime},E).$
(2)

The summation over $j$ accounts for the different gas elements (mainly
molecular and atomic hydrogen, and helium), while that over $i$ accounts for
the different UHECR-projectile elements. The $n_{j}(\mathbf{x})$ function
stands for the position-dependent number density of the $j$th element,
$\phi_{i}(E^{\prime})$ is the energy-differential flux per steradian (energy
spectrum) of the $i$th UHECR element with energy $E^{\prime}$,
$\sigma_{ij}(E^{\prime})$ is the inelastic cross section of the collision
between the projectile $i$ and the target $j$, and
$\mathrm{d}N_{ij}^{\gamma}/\mathrm{d}E(E^{\prime},E)$ is the inclusive
spectrum of photons generated in the $ij$ collision subsequently to the decay
of all unstable particles created in the collision. The integration is carried
out over all UHECR energies $E^{\prime}>E$ that allow for generating photons
with energy $E$. Finally, the $4\pi$ factor results from the integration of
the (isotropic) UHECR flux over solid angle.

The aim of this paper is thus to estimate $\phi_{\gamma}(E,\mathbf{n})$ by
means of Equations 1 and 2. To do so, the UHECR flux shaped by the various
contributing mass components is described in Section 2; a survey of the up-to-
date modeling of the interstellar gas density in the disk of the Milky Way is
presented in Section 3; and the photon production in cosmic ray-gas collisions
is detailed in Section 4. The diffuse flux of UHE photons is then obtained in
Section 5 for various energy thresholds, before implications for the searches
for super-heavy dark matter decaying into photons are drawn in Section 6.
Finally, the significance of the results is summarized in Section 7.

##  2 Ultra-high energy cosmic-ray flux

Figure 1: All-particle spectrum of UHECRs in thick grey, as reported in (Aab
et al., 2021). The various mass-discriminated energy spectra are obtained by
making use of the Sibyll2.3 model of hadronic interactions to interpret the
longitudinal profiles of extensive air showers observed at the Pierre Auger
Observatory (Aab et al., 2014; Bellido, 2018).

In the energy range of interest, knowledge of the energy spectrum and mass
composition of UHECRs is possible only through indirect measurements, via the
extensive air showers that they produce in the atmosphere. In particular, the
energy-dependent mass composition is accessible on a statistical basis only.
Hence, the various $\phi_{i}(E)$ components must be inferred by combining the
all-particle energy spectrum of UHECRs, $\phi(E)=\sum_{i}\phi_{i}(E)$, and the
abundances of the different elements as a function of energy.

The all-particle energy spectrum of UHECRs $\phi(E)$ has been measured between
$10^{17}~{}$and $10^{18}~{}$eV from several dedicated contemporary experiments
(Knurenko et al., 2013; Bertaina et al., 2016; Abbasi et al., 2018; Aartsen et
al., 2019; Budnev et al., 2020). At higher energy, above $10^{18}~{}$eV, the
fall of the flux required a leap in cumulative exposure to collect an
increased influx of events (Aab et al., 2015; Abu-Zayyad et al., 2012). The
measurement of the spectrum carried out at the Pierre Auger Observatory above
$10^{17}~{}$eV, reported in (Aab et al., 2021), is the one relying on the
largest cumulated exposure with a single detector type, avoiding the necessity
to combine measurements, the matching of which inevitably implies additional
systematic effects. In addition, the energy scale is set calorimetrically via
fluorescence telescopes, avoiding assumptions on the mass of the particles. We
use this measurement in the following, shown as the thick grey line in Fig. 1,
and make an explicit use of the reported covariance matrix
$\mathbf{\sigma_{\phi}}$ that accounts for both statistical and systematic
uncertainties to calculate (see Appendix B) the 2-sided 16% quantiles of the
underlying distribution of spectrum estimator, denoted as $\phi_{+}(E)$ and
$\phi_{-}(E)$ (shown as the thin grey lines in Fig. 1). This will allow us to
propagate in Section 5 the uncertainties in the spectrum measurement, coming
for the most part from those in the energy scale, in the calculation of the
photon flux.

The approximately power-law shape of the spectrum in $E^{-3}$ hides beneath a
complex intertwining of different astrophysical phenomena, which might be
exposed by looking at the spectrum of different primary elements. An
illustration of the contributions of protons, He, N and Fe nuclei to the total
spectrum is shown, by scaling the He (N) [Fe] components by 10 (100) [1000] so
as to distinguish each of them. They are obtained from the abundances of
elements reported in Aab et al. (2014); Bellido (2018), by weighting the total
spectrum with the fractions of the corresponding mass group as inferred from
the distributions of the depths of the shower maximum, $X_{\text{max}}$,
measured at the Pierre Auger Observatory and modelled according to, for this
specific illustration (Fig. 1), the hadronic model Sibyll2.3 (Riehn et al.,
2017). A linear interpolation (in $\log_{10}{E}$) is applied to smooth the
fluxes at the bin-center values. One can observe the steep fall-off of the Fe
component above $10^{17}$ eV, which is along the lines of the long-standing
scenario for Galactic CRs characterised by a rigidity-dependent maximum
acceleration energy (Peters, 1961). On the other hand, one can also observe
the presence of protons, He, and N nuclei, at least up to $\simeq 10^{19}$ eV,
at which point the proton flux begins to quench and that of the N hardens.
This is consistent with a rigidity-dependent maximum acceleration energy
scenario for extra-galactic sources, see e.g. Aloisio et al. (2014); Taylor et
al. (2015); Aab et al. (2017, 2020b). Within such a scenario, the presence of
intermediate He and N nuclei in the energy range between $\simeq 10^{17}$ and
$\simeq 10^{18.7}$ eV, where the Galactic Fe is declining and presumably the
extra-galactic protons are rising, could be fueled by a second Galactic “B
component”, as first suggested in Hillas (2005).

With caution over the interpretation of the $X_{\text{max}}$ measurements, we
thus use the data described above to model the various $\phi_{i}(E)$
components. As a proxy of the estimation of the uncertainties, we use the
three different results reported in Bellido (2018) for the fractions of mass
groups inferred from three different hadronic models, namely EPOS-LHC (Pierog
et al., 2015), QGSJetII-04 (Ostapchenko, 2013) and Sibyll2.3. The systematic
differences are on the order of, or larger than, the statistical
uncertainties, which, in the absence of an available covariance matrix
governing the correlations between energy bins, are not accounted for in this
study.

Our final goal is to calculate photon fluxes above $10^{17}$ eV, which are
mainly produced by cosmic rays of energy larger than $A~{}10^{18}$ eV, with
$A$ the atomic number of the particles. Hence mainly extra-galactic UHECRs,
the flux of which is assumed to illuminate the Galaxy uniformly in the
relevant energy range, can contribute to the photon fluxes sought for.
Subsequently, the various $\phi_{i}(E)$ components are assumed to fill
homogeneously the Galactic volume. This neglects flux-amplification effects
preferentially in the disk that would be indicative of larger residence times
for particles experiencing stronger values of the Galactic magnetic field.
This is reasonable as numerical simulations of particle trajectories in the
field show that the residence time of cosmic rays tends to equate their free
escape time for particle rigidities above $\simeq 1~{}$EV (Zirakashvili et
al., 1998; Kaapa et al., 2021).

##  3 Interstellar gas density in the Milky Way

The interstellar medium in the Galaxy is known to provide a mass of a few
$10^{9}~{}M_{\sun}$, which is distributed predominantly in the disk of the
Milky Way. Hydrogen ($\simeq 90\%$) and helium ($\simeq 10\%$) in gaseous
state make up most of this mass, while a small fraction is contributed by
silicate and graphite dust. Depending on its temperature, three forms of the
hydrogen gas are at play: atomic (H i), molecular (H2), and ionized (H ii)
hydrogen, while helium remains mostly neutral due to its much higher first
ionization potential. In this study, we use two models of the spatial
distribution of the gas recently developed from existing data. The differences
in the final photon emission from UHECR-gas interactions between the two
models are interpreted as contributing to the systematic uncertainties of
$\phi_{\gamma}(E,\mathbf{n})$. The first model, developed by Lipari & Vernetto
(2018) and dubbed as model A hereafter, aims at capturing the large-scale
properties of the gas based on an axially and up-down symmetric distribution
with a scale height increasing with the radial distance from the Galactic
center, hence without attempting to describe accurately fine details. It is
worth noting the diffuse gamma-ray fluxes observed around GeV energies and
sub-PeV ones could be reproduced using this model. The second model, developed
by Jóhannesson et al. (2018) and dubbed as model B, incorporates in addition
spiral arms and accounts for the warping of the disk so as to constrain models
of cosmic-ray propagation and to trace the structure of the Galaxy. We
describe below the main features of the models, without reproducing the
detailed parameterizations available in the original papers. A comprehensive
description of the data used to construct these models can be found in, e.g.,
Ferriere (2001).

The H i component represents a large fraction of the mass. The most important
probe of its distribution is the line at 21 cm observed in emission or
absorption. The density is inferred to be constant in the distance range from
4 to $\simeq 10~{}$kpc from the Galactic center and decreases steadily at
larger distances, but falls exponentially in the vertical direction with a
scale length of $\simeq 250$ pc for the warm component and $\simeq 130$ pc for
the cold one. The density of the cold component is highly non-uniform, as
illustrated by its filling factor of the volume of the thin disk less than
0.1.

Although almost equally abundant in the interstellar medium, the distribution
of H2 is less well known. This is because such a homonuclear diatomic molecule
does not possess a permanent electric dipole moment, and therefore has no
simple rotational transition. It can be however observed in the UV region
through electronic transitions, and indirectly inferred from the rotational
transition of the 12C16O molecule (CO hereafter for convenience), which is
mostly excited through collisions with H2 and which is the second most
abundant molecule in the interstellar medium. The integrated intensity of the
CO lines turns out to be almost linearly related to the column density of H2
so that the CO observations at radio wavelengths nicely complement the surveys
of the 21-cm line, the proportionality factor being taken as
$2{\times}10^{20}~{}$cm${}^{-2}~{}$(K km s-1)${}^{-1}~{}$(Bolatto et al.,
2013).

Figure 2: Surface density average over Galactic longitude as a function of
Galactic latitude for neutral (H i) and molecular (H2) hydrogen, as modeled by
Lipari & Vernetto (2018) (model A) and by Jóhannesson et al. (2018) (model B).
The surface density of H2 as obtained from the composite survey of Dame et al.
(2001) is shown as the orange line.

The surface densities for neutral and molecular hydrogen, as inferred from the
models A and B and averaged over Galactic longitude, are shown in Fig. 2. For
reference, the surface density of H2 as obtained from the composite survey of
Dame et al. (2001) is shown as the orange line. Once summing the H i and H2
components, the two models give rise to some differences by a factor of a few
in units of $10^{22}~{}$cm-2. We consider these differences as sources of
systematic uncertainties in the photon flux calculation.

The distribution of the helium contribution, $\simeq 10\%$ of that of
hydrogen, is assumed to follow closely the hydrogen one. No additional gas
element is considered in this study. In particular, the small contribution
from ionized hydrogen is neglected.

##  4 Photon production in cosmic ray-gas collisions

Figure 3: Inelastic cross section as a function of energy for the various
processes of interest. Figure 4: Photon yield per UHECR interaction for
different primaries and incoming energies.

Cosmic rays interact with the interstellar medium through the different
processes described in the introduction. These collisions create neutral pions
$\pi^{0}$ whose most probable mode of decay produces two photons:

$\pi^{0}\rightarrow 2\gamma.$ (3)

The production of secondary particles such as $K$, $\rho$, $\eta$ decaying
into neutral pions also occurs during these interactions, resulting in an
additional flux of photons.

To obtain the relevant inelastic cross sections $\sigma_{ij}(E^{\prime})$ as
well as the inclusive spectrum of photons
$\mathrm{d}N^{\gamma}_{ij}/\mathrm{d}E(E^{\prime},E)$, which corresponds to
the mean number of photons in the energy range $[E,E+\mathrm{d}E]$ produced in
a single interaction between the cosmic ray element $i$ with energy
$E^{\prime}$ and the gas element $j$ at rest, we use Cosmic Ray Monte Carlo
(CRMC) (Ulrich et al., 2021). This package provides different cosmic-ray event
generators that model the production of secondary particles occurring in
hadronic interactions. It also allows the follow-up of the decays of the
secondary particles that are likely to produce photons. In this study, the
hadronic interaction model EPOS-LHC (Pierog et al., 2015) is chosen.

Thus, for different values of the UHECR energies ranging from $10^{17}$ eV to
$10^{20}$ eV (in steps of $\Delta\log_{10}(E^{\prime})=0.5$), $100,000$
collisions are simulated for each couple of cosmic ray and target element
$(i,j)$ of interest here. The spectra of the photon energies are then produced
for each cosmic-ray energy, with a bin width of $\Delta\log_{10}(E)=0.05$.

The inelastic cross sections $\sigma_{ij}(E)$ obtained with CRMC are shown in
Fig. 3 for each cosmic-ray element as continuous lines for H targets and
dashed ones for He targets. They range typically in the hundred of millibarns,
reaching thousand ones for the heaviest collisions. The photon yields are
shown in Fig. 4 for the different primaries and three cosmic-ray energies,
$10^{9}~{}$GeV (dotted), $10^{10}~{}$GeV (dashed) and $10^{11}~{}$GeV
(continuous). For a fixed photon energy $E<E^{\prime}$, the expected increase
of the yield with the incident energy of cosmic rays is observed. On the other
hand, for a fixed cosmic-ray energy, the yield is also observed to depend
strongly on the cosmic-ray mass. Given that a single nucleus generally
intervenes in each interaction, this behavior is expected from the energy
available in each nucleus, which, compared to the total energy of the nucleus,
is reduced by the atomic number of the nucleus. As a result, the photon fluxes
are thus more copiously produced at higher energies by the lightest primaries.

##  5 Diffuse flux of ultra-high energy photons

On inserting the various ingredients detailed in Sections 2, 3, and 4 into
equation 2 and subsequently equation 1, the photon flux per energy unit sought
for is obtained. In Fig. 5 and Fig. 6 are shown the results by considering the
two models of gas distribution in the Galactic disk described in Section 3,
and by integrating over energy above $10^{17}~{}$eV. As expected, the flux is
concentrated a few degrees around the Galactic plane, amounting to $\simeq
8.7{\times}10^{-2}$ km-2yr-1sr-1 once averaged out over $|b|\leq 5^{\circ}$ in
the case of model A ($\simeq 3.2{\times}10^{-2}$ km-2yr-1sr-1 in the case of
model B). The pattern is brighter in the innermost region of the disk for both
models.

Figure 5: Photon flux in Galactic coordinates (Hammer projection) expected
from UHECR-gas interactions, integrated above $10^{17}~{}$eV. Model A is used
for the gas distribution in the Galactic disk. Figure 6: Same as Fig. 5, using
model B for the gas distribution in the Galactic disk.

The energy dependence of the photon flux is shaped by that of the cosmic-ray
flux, inheriting from its main features but shifted about a decade earlier.
Averaged over a 5∘-band around the Galactic plane, the photon flux is about
$10^{-5}$ that of the UHECRs for energy thresholds ranging from $10^{17}~{}$eV
to $10^{18}~{}$eV, before steepening to $10^{-6}$ above $10^{19}~{}$eV and
dropping sharply at higher thresholds.

Figure 7: Upper limits on diffuse photon fluxes reported by Fomin et al.
(2017) (pink triangles, EAS-MSU), Apel et al. (2017) (filled orange circles,
KASCADE-Grande), Abbasi et al. (2019) (open circles, Telescope Array), and
Savina & Pierre Auger Collaboration (2021) (dark blue squares, Auger) for
several energy thresholds. Expected fluxes from UHECRs interactions with the
gas distribution in the Galaxy are shown as the grey dashed line (model A) and
grey dashed-dotted line (model B), as well as from UHECRs interactions with
background photon fields (GZK mixed) estimated by Bobrikova et al. (2021).

The calculated photon flux can be compared to several existing upper limits
above $10^{17}~{}$eV from several observatories. Shown in Fig. 7 as pink
triangles (Fomin et al., 2017), filled orange circles (Apel et al., 2017),
dark blue squares (Savina & Pierre Auger Collaboration, 2021), and open
circles (Abbasi et al., 2019), they appear to be between $\simeq 2.5$ and
$\simeq 3$ orders of magnitude above our expectations for gas distribution
models A and B for energy thresholds between $10^{17}~{}$eV (amounting to
$\Phi_{\gamma}\simeq 1.1{\times}10^{-2}~{}$km${}^{-2}~{}$yr${}^{-1}~{}$sr-1
for model A) and $\simeq 10^{18.5}~{}$eV ($\Phi_{\gamma}\simeq
2.0{\times}10^{-6}~{}$km${}^{-2}~{}$yr${}^{-1}~{}$sr-1), and even higher for
larger thresholds. The bands correspond to the systematic uncertainties
arising from those in the all-particle energy spectrum and those of the
hadronic interaction model used to infer the mass composition of UHECRs. Note
that to be compared to the upper limits on the diffuse fluxes illuminating the
whole field of view of the observatories considered, the expected fluxes A and
B reported in Fig. 7 are “calibrated” on the whole sphere as follows:

$\Phi_{\gamma}(E)=\frac{1}{4\pi}\int_{E}^{\infty}\mathrm{d}E^{\prime}\int_{4\pi}\mathrm{d}\mathbf{n}~{}\phi_{\gamma}(E^{\prime},\mathbf{n}).$
(4)

For comparison with other astrophysical expectations, another photon flux,
resulting from interactions of extragalactic UHECRs with the background photon
fields permeating the Universe, is shown as the orange band (Bobrikova et al.,
2021). The hadrons that cause the creation of $\pi^{0}$ mesons, which then
decay into photons, may be primary proton cosmic rays, or secondary ones
produced by the photo-disintegration of nuclei interacting inelastically with
a cosmic background photon. Since the nucleons produced inherit the energy of
the fragmented nucleus divided by its atomic number, this photon flux
therefore depends on the nature of the UHECRs. Other dependencies come from
the maximum acceleration energy of the nuclei at the sources and the shape of
the energy spectrum of the accelerated particles. Different assumptions about
these ingredients have led to different estimates in the literature, see e.g.
Kampert et al. (2011). The recent estimate reported in Bobrikova et al. (2021)
that considers matching the energy spectrum measured at the Pierre Auger
Observatory simultaneously with the mass composition inferred from these data
is considered here. Considered as systematic uncertainties for the expected
photon flux, they give the orange band in Fig. 7.

Interestingly, among the two photon fluxes produced by cosmic-ray
interactions, the one explored in this study appears to be dominant in the
energy range between $10^{17}~{}$and $\simeq 10^{18}~{}$eV. Its observational
signature, that of an excess of a few degrees around the Galactic plane, makes
it unambiguous to identify once the required sensitivity is reached in the
future. At higher energy thresholds, the other one, expected to be isotropic,
takes over.

##  6 Implications for searches for super-heavy dark matter

Between $10^{17}~{}$and $\simeq 10^{18}~{}$eV, the flux of UHE photons
reported in this study can be seen as a floor that prevents probing sources in
the Galactic disk or patterns of super-heavy dark matter (SHDM) decaying
primarily from the Galactic center and producing a photon flux
$\phi_{\gamma}^{\mathrm{DM}}$ below the level of $\phi_{\gamma}$. This is
because the detection of such an overwhelmed flux would require an exposure
much larger than $1/\phi_{\gamma}^{\mathrm{DM}}$ so as to pick up a
significant excess above $\phi_{\gamma}$ integrated within some angular and
energy ranges. The production of super-heavy particles in the early Universe
indeed remains a possible solution to the dark matter puzzle because of the
high value of the instability energy scale in the Standard Model of particle
physics, which, according to current measurements of the Higgs boson mass and
the Yukawa coupling of the top quark, ranges between $10^{10}$ and
$10^{12}~{}$GeV (Degrassi et al., 2012; Bednyakov et al., 2015). The Standard
Model can therefore be extrapolated without encountering inconsistencies that
would make the electroweak vacuum unstable up to such energy scales (and even
to much higher ones given the slow evolution of the instability scale up to
the Planck mass (Degrassi et al., 2012)) where new physics could arise, giving
rise to a mass spectrum of super-heavy particles that could have been produced
during post-inflation reheating by various mechanisms, see e.g. Ellis et al.
(1990, 1992); Berezinsky et al. (1997); Chung et al. (1998); Garny et al.
(2016); Ellis et al. (2016); Dudas et al. (2017); Kaneta et al. (2019);
Mambrini & Olive (2021). The existence of a photon-flux floor in the direction
of the Galactic center has important implications for SHDM searches using UHE
photons from the Galactic center, and translates, depending on their mass
$M_{X}$, into a “ceiling” for the exploration of the particle lifetime that we
now discuss.

Figure 8: Allowed region of mass and lifetime of super-heavy dark matter
particles decaying into standard-model ones. The filled red region is excluded
from the upper limits in UHE photon fluxes while the hatched one corresponds
to the observational ceiling discussed in the text.

Secondary UHE photons produced by the decay of SHDM particles can be observed
today if these particles have a long enough lifetime. This usually requires a
fine tuning between the coupling constant of the particle with those produced
and the dimension of the operator governing the decay (de Vega & Sánchez,
2003). Alternatively, SHDM particles can be protected from decay in the
perturbative domain by the conservation of new quantum number(s) but
disintegrate by instanton non-perturbative effects (Kuzmin & Rubakov, 1998).
This mechanism offers the possibility of providing metastable particles with
lifetime $\tau_{X}$, which can produce detectable secondaries such as nucleons
and photons. In particular, the flux (per steradian) of photons produced can
then be obtained as the sum of a contribution from the Galaxy and of an
isotropic extragalactic one – see e.g. Berezinsky & Mikhailov (1999); Benson
et al. (1999); Medina-Tanco & Watson (1999); Aloisio et al. (2006); Siffert et
al. (2007); Kalashev & Kuznetsov (2016); Alcantara et al. (2019). The Galactic
contribution, which is the dominating one, is obtained by integrating the
position-dependent emission rate per unit volume and unit energy along the
path of the photons in the direction $\mathbf{n}$:

$\phi_{\gamma}^{\mathrm{DM}}(E,\mathbf{n})=\frac{1}{4\pi
M_{X}\tau_{X}}\frac{\mathrm{d}N}{\mathrm{d}E}\int_{0}^{\infty}\mathrm{d}s\rho_{\mathrm{DM}}(\mathbf{x}_{\sun}+s\mathbf{n}).$
(5)

Here, $\rho_{\mathrm{DM}}$ is the energy-density profile of dark matter in the
Galaxy. There are uncertainties in the determination of this profile (Guépin
et al., 2021). However, the typical $1$-degree angular resolution of UHECR-
Observatory capabilities cannot probe the uncertainties of the parametric
models Alcantara et al. (2019). We thus assume here the traditional Navarro-
Frenk-White profile (Navarro et al., 1996), normalized so that
$\rho(\mathbf{x}_{\sun})=0.3~{}$GeV cm-3 in the solar system. The other term,
$\mathrm{d}N/\mathrm{d}E$, stands for the energy spectrum of the secondary
photons, which depends on the hadronization process. Various computational
schemes, see e.g. Aloisio et al. (2004), show that the spectrum of the final
state particles (photons, nucleons, and neutrinos) varies as $E^{-1.9}$. We
use this widely-agreed energy dependence, normalizing the spectrum in the same
way as in Alcantara et al. (2019)

From the limits on $J_{\gamma}(E)$, constraints can be inferred in the plane
$(\tau_{X},M_{X})$ by requiring the all-direction flux calculated by averaging
Eqn. 5 over all directions to be less than the limits. For a specific upper
limit, a scan of the value of the mass $M_{X}$ is carried out so as to infer a
lower limit of the $\tau_{X}$ parameter. This defines a curve. By repeating
the procedure for each upper limit on $J_{\gamma}(E)$, a set of curves is
obtained, reflecting the sensitivity of a specific energy threshold to some
range of mass. The union of the excluded region is shown as the filled red
region in Fig. 8, which is similar to that of previous studies (Kalashev &
Kuznetsov, 2016; Ishiwata et al., 2020).

Although the directional dependence of $\phi_{\gamma}(E,\mathbf{n})$ (Eqn. 1)
is different from that of $\phi_{\gamma}^{\mathrm{DM}}(E,\mathbf{n})$ (a hot
spot falling steeply a few degrees away from the Galactic center), it shall
mask the presence of decaying SHDM for values of
$\phi_{\gamma}^{\mathrm{DM}}(E,\mathbf{n})$ in the direction of the Galactic
center below that of $\phi_{\gamma}(E,\mathbf{n})$ in the same direction. By
scanning $M_{X}$, we can also determine the corresponding value of $\tau_{X}$
beyond which the expected values of $\phi_{\gamma}^{\mathrm{DM}}(E)$ would be
overwhelmed by the $\phi_{\gamma}(E)$ averaged-directional fluxes. This
results in the ceiling region shown as the hatched region in Fig. 8. The
ceiling is observed to affect mass values up $M_{X}\simeq 10^{11}~{}$GeV,
reflecting the cut-off of the photon fluxes produced by cosmic-ray
interactions. Below $M_{X}\simeq 10^{10}~{}$GeV, the ceiling region is most
constrained by the photon flux inferred in this study. For larger masses, on
the other hand, the constraints obtained from the photon flux inferred in
Bobrikova et al. (2021) take over. The sharp increase just above
$M_{X}=10^{11}~{}$GeV reflects the sharp drop of the photon fluxes originating
from cosmic-ray interactions for energy thresholds above a few times
$10^{19}~{}$eV, leaving a large region of the phase space unaffected by the
ceiling.

##  7 Summary

The all-sky flux a few degrees around the Galactic plane of UHE photons
produced by UHECR interactions in the disk of the Galaxy amounts to,
integrated above $10^{17}~{}$eV,
$1.1{\times}10^{-2}~{}$km${}^{-2}~{}$yr${}^{-1}~{}$sr-1. It turns out to be
the dominating one from cosmic-ray interactions for energies between
$10^{17}~{}$and $\simeq 10^{18}~{}$eV. Although out of reach of the current
and planned observatories, it may be responsible for the first detection of
photons in this energy range once future experiments reach the required
exposure, unless brighter sources of photons are discovered while cumulating
this exposure. This flux will also be a floor below which other signals coming
from the same direction would be overwhelmed. This is relevant for searches
for photon fluxes from the Galactic center that would be indicative of the
decay of super-heavy dark matter particles, and this translates into a ceiling
region for the lifetime of the putative particles to be explored below
$M_{X}\simeq 10^{11}~{}$GeV.

Equally expected from UHECR-gas interactions, UHE fluxes of neutrinos are
produced by charged pions or neutrons decays. Quantitative estimates are left
for a future study.

## Acknowledgments

We thank Roger Clay for his comments. This work was made possible by with the
support of the Institut Pascal at Université Paris-Saclay during the Paris-
Saclay Astroparticle Symposium 2021, with the support of the P2IO Laboratory
of Excellence (program “Investissements d’avenir” ANR-11-IDEX-0003-01 Paris-
Saclay and ANR-10-LABX-0038), the P2I axis of the Graduate School Physics of
Université Paris-Saclay, as well as IJCLab, CEA, IPhT, APPEC, the IN2P3 master
projet UCMN and EuCAPT ANR-11-IDEX-0003-01 Paris-Saclay and ANR-10-LABX-0038).

##  Appendix A UHE photon absorption in the Galactic disk

The radiation length of UHE photons propagating in the matter of the Galactic
disk is of order $100~{}$g cm-2, much larger than the column density that is
traversed, which is at most a few tenths of g cm-2. UHE photon absorption by
the dust is thus negligible. The most important process is pair production
$\gamma\gamma\rightarrow e^{+}e^{-}$ with the radiation fields that permeate
the Galaxy. In the energy range considered in this study ($\geq
10^{17}~{}$eV), the radiation field that dominates by far is that of the
cosmic microwave background (Lipari & Vernetto, 2018). In this case, because
the density of the radiation field is uniform in space, the optical depth for
photons observed on Earth with energy $E$ in the direction $\mathbf{n}$,
$\tau(E,\mathbf{n},x)$, is proportional to the travelled distance $x$:

$\tau(E,\mathbf{n},x)=K(E,-\mathbf{n})x.$ (A1)

The function $K(E,\mathbf{n})$ is the interaction probability per unit length,
obtained by integrating over the energy and the (isotropic) angular
distribution of the target photons:

$K(E,\mathbf{n})=\int\mathrm{d}\mathbf{\kappa}(1-\cos{\theta_{\gamma\gamma}})n_{\gamma}(\epsilon)\sigma_{\gamma\gamma}(s).$
(A2)

In this expression, $\mathbf{\kappa}$ is the 3-momentum of the target photon
and $\epsilon=|\mathbf{\kappa}|$ its energy, $\theta_{\gamma\gamma}$ is the
angle between the interacting photons, and $\sigma_{\gamma\gamma}(s)$ is the
pair-production cross section expressed as a function of the square of the
center-of-mass energy $s=2E\epsilon(1-\cos{\theta_{\gamma\gamma}})$, which, in
the high-energy regime that is the relevant one here, behaves, modulo a
logarithmic correction, as $1/s$ (and thus as $1/E$):

$\sigma_{\gamma\gamma}(s)\simeq\frac{3\sigma_{\mathrm{T}}m_{\mathrm{e}}^{2}c^{4}}{2s}\left(\ln{\frac{s}{m_{\mathrm{e}}^{2}c^{4}}}-1\right),$
(A3)

with $\sigma_{\mathrm{T}}$ the Thomson cross section. In this way, the optical
depth averaged over angle behaves roughly as $\tau(E)\simeq 0.03/E$ (in
kpc-1). Plugging an absorption factor $\exp{(-\tau(E,\mathbf{n}))}$ into the
integrand in equation 1 leads however to small reductions, of order of 10% at
$10^{17}~{}$eV and much smaller at higher energies. This is much less than the
various systematic uncertainties considered in the study.

##  Appendix B Propagation of the uncertainties in the all-particle energy
spectrum

Denoting as a vector the set of measurements of the energy spectrum in each
energy bin, $\boldsymbol{\phi}=\\{\phi_{1},\phi_{2},...,\phi_{N}\\}$, the
$\boldsymbol{\phi}_{+}$ vector is defined as the set of values that satisfy

$\frac{1}{\sqrt{(2\pi)^{N}\mathrm{det}\boldsymbol{\sigma_{\phi}}}}\int_{\boldsymbol{\phi}_{+}}^{\infty}\mathrm{d}\boldsymbol{\phi}\exp{\left(-\frac{1}{2}\delta\boldsymbol{\phi}^{\mathrm{T}}\boldsymbol{\sigma^{-1}_{\phi}}\delta\boldsymbol{\phi}\right)}=0.84.$
(B1)

The notation $\delta\boldsymbol{\phi}$ stands for a random fluctuation around
the set of observed values. To solve equation B1 for the unknown
$\boldsymbol{\phi}_{+}$ and the corresponding one for $\boldsymbol{\phi}_{-}$,
we build the probability distribution function of $\boldsymbol{\phi}$ by
“whitening”, or “standardizing” the covariance matrix in the following way. We
first decompose the covariance matrix as
$\boldsymbol{\sigma_{\phi}}=\boldsymbol{R}\boldsymbol{\sigma_{\phi}^{\prime}}\boldsymbol{R}^{\mathrm{T}}$,
where the columns of the rotation matrix $\boldsymbol{R}$ are the eigenvectors
of $\boldsymbol{\sigma_{\phi}}$ and
$\boldsymbol{\sigma_{\phi}^{\prime}}=\mathrm{diag}(\sigma_{1}^{\prime
2},\dots,\sigma_{N}^{\prime 2})$, with $\sigma_{i}^{\prime 2}$ the eigenvalues
of $\boldsymbol{\sigma_{\phi}}$. Then, it is straightforward to build a
transformed flux $\delta\boldsymbol{\phi}^{\prime}$ for which the covariance
matrix is the identity one by considering the re-scaled flux
$\delta\boldsymbol{\phi}^{\prime}=\boldsymbol{\sigma_{\phi}^{\prime}}^{-1/2}\boldsymbol{R}^{\mathrm{T}}\delta\boldsymbol{\phi}$.
In this way, the directional properties of the initial covariance matrix are
eliminated in the new coordinate system so that random realizations of
$\delta\boldsymbol{\phi}^{\prime}$ can be drawn from a $N$D normal
distribution, with the set of values
$(\delta\phi^{\prime}_{1},\dots,\delta\phi^{\prime}_{N})$ independent from
each other. The corresponding values for $\delta\boldsymbol{\phi}$ are
recovered from the transformation
$\delta\boldsymbol{\phi}=\boldsymbol{R}\boldsymbol{\sigma_{\phi}^{\prime}}^{1/2}\delta\boldsymbol{\phi}^{\prime}$.
By repeating a large number of times the procedure, the 2-sided 16% quantiles
defining $\boldsymbol{\phi}_{+}$ and $\boldsymbol{\phi}_{-}$ are finally
estimated.

## References

  * Aab et al. (2014) Aab, A., et al. 2014, Phys. Rev. D, 90, 122006, doi: 10.1103/PhysRevD.90.122006
  * Aab et al. (2015) —. 2015, Nucl. Instrum. Meth. A, 798, 172, doi: 10.1016/j.nima.2015.06.058
  * Aab et al. (2017) —. 2017, JCAP, 04, 038, doi: 10.1088/1475-7516/2017/04/038
  * Aab et al. (2020a) —. 2020a, Astrophys. J., 891, 142, doi: 10.3847/1538-4357/ab7236
  * Aab et al. (2020b) —. 2020b, Phys. Rev. Lett., 125, 121106, doi: 10.1103/PhysRevLett.125.121106
  * Aab et al. (2021) —. 2021, Eur. Phys. J. C, 81, 966, doi: 10.1140/epjc/s10052-021-09700-w
  * Aartsen et al. (2019) Aartsen, M. G., et al. 2019, Phys. Rev. D, 100, 082002, doi: 10.1103/PhysRevD.100.082002
  * Abbasi et al. (2018) Abbasi, R. U., et al. 2018, Astrophys. J., 865, 74, doi: 10.3847/1538-4357/aada05
  * Abbasi et al. (2019) —. 2019, Astropart. Phys., 110, 8, doi: 10.1016/j.astropartphys.2019.03.003
  * Abu-Zayyad et al. (2012) Abu-Zayyad, T., et al. 2012, Nucl. Instrum. Meth. A, 689, 187, doi: 10.1016/j.nima.2012.05.079
  * Ackermann et al. (2012) Ackermann, M., et al. 2012, Astrophys. J., 750, 3, doi: 10.1088/0004-637X/750/1/3
  * Alcantara et al. (2019) Alcantara, E., Anchordoqui, L. A., & Soriano, J. F. 2019, Phys. Rev. D, 99, 103016, doi: 10.1103/PhysRevD.99.103016
  * Aloisio et al. (2014) Aloisio, R., Berezinsky, V., & Blasi, P. 2014, JCAP, 10, 020, doi: 10.1088/1475-7516/2014/10/020
  * Aloisio et al. (2004) Aloisio, R., Berezinsky, V., & Kachelriess, M. 2004, Nucl. Phys. B Proc. Suppl., 136, 319, doi: 10.1016/j.nuclphysbps.2004.10.023
  * Aloisio et al. (2006) —. 2006, Phys. Rev. D, 74, 023516, doi: 10.1103/PhysRevD.74.023516
  * Amenomori et al. (2021) Amenomori, M., et al. 2021, Phys. Rev. Lett., 126, 141101, doi: 10.1103/PhysRevLett.126.141101
  * Apel et al. (2017) Apel, W. D., et al. 2017, Astrophys. J., 848, 1, doi: 10.3847/1538-4357/aa8bb7
  * Atkins et al. (2005) Atkins, R., et al. 2005, Phys. Rev. Lett., 95, 251103, doi: 10.1103/PhysRevLett.95.251103
  * Bartoli et al. (2015) Bartoli, B., et al. 2015, Astrophys. J., 806, 20, doi: 10.1088/0004-637X/806/1/20
  * Bednyakov et al. (2015) Bednyakov, A. V., Kniehl, B. A., Pikelner, A. F., & Veretin, O. L. 2015, Phys. Rev. Lett., 115, 201802, doi: 10.1103/PhysRevLett.115.201802
  * Bellido (2018) Bellido, J. 2018, PoS, ICRC2017, 506, doi: 10.22323/1.301.0506
  * Benson et al. (1999) Benson, A., Wolfendale, A. W., & Smialkowski, A. 1999, Astropart. Phys., 10, 313, doi: 10.1016/S0927-6505(98)00065-6
  * Berezinsky et al. (1997) Berezinsky, V., Kachelrieß, M., & Vilenkin, A. 1997, Phys. Rev. Lett., 79, 4302, doi: 10.1103/PhysRevLett.79.4302
  * Berezinsky & Mikhailov (1999) Berezinsky, V., & Mikhailov, A. A. 1999, Phys. Lett. B, 449, 237, doi: 10.1016/S0370-2693(99)00053-2
  * Berezinsky et al. (1993) Berezinsky, V. S., Gaisser, T. K., Halzen, F., & Stanev, T. 1993, Astroparticle Physics, 1, 281, doi: 10.1016/0927-6505(93)90014-5
  * Bertaina et al. (2016) Bertaina, M., et al. 2016, PoS, ICRC2015, 359, doi: 10.22323/1.236.0359
  * Bloemen (1989) Bloemen, H. 1989, ARA&A, 27, 469, doi: 10.1146/annurev.aa.27.090189.002345
  * Bloemen et al. (1984) Bloemen, J. B. G. M., Blitz, L., & Hermsen, W. 1984, ApJ, 279, 136, doi: 10.1086/161873
  * Bobrikova et al. (2021) Bobrikova, A., Niechciol, M., Risse, M., & Ruehl, P. 2021, Proc. 37th International Cosmic Ray Conference (Berlin, Germany), PoS (ICRC2021), 449, doi: 10.22323/1.395.0449
  * Bolatto et al. (2013) Bolatto, A. D., Wolfire, M., & Leroy, A. K. 2013, Ann. Rev. Astron. Astrophys., 51, 207, doi: 10.1146/annurev-astro-082812-140944
  * Budnev et al. (2020) Budnev, N. M., et al. 2020, Astropart. Phys., 117, 102406, doi: 10.1016/j.astropartphys.2019.102406
  * Carceller & Masip (2017) Carceller, J. M., & Masip, M. 2017, J. Cosmology Astropart. Phys, 2017, 013, doi: 10.1088/1475-7516/2017/03/013
  * Chardonnet et al. (1995) Chardonnet, P., Salati, P., Silk, J., Grenier, I., & Smoot, G. 1995, ApJ, 454, 774, doi: 10.1086/176530
  * Chung et al. (1998) Chung, D. J., Kolb, E. W., & Riotto, A. 1998, Phys. Rev. D, 59, 023501, doi: 10.1103/PhysRevD.59.023501
  * Dame et al. (2001) Dame, T. M., Hartmann, D., & Thaddeus, P. 2001, Astrophys. J., 547, 792, doi: 10.1086/318388
  * de Vega & Sánchez (2003) de Vega, H. J., & Sánchez, N. G. 2003, Phys. Rev. D, 67, 125019, doi: 10.1103/PhysRevD.67.125019
  * Degrassi et al. (2012) Degrassi, G., Di Vita, S., Elias-Miro, J., et al. 2012, JHEP, 08, 098, doi: 10.1007/JHEP08(2012)098
  * Dudas et al. (2017) Dudas, E., Mambrini, Y., & Olive, K. 2017, Phys. Rev. Lett., 119, 051801, doi: 10.1103/PhysRevLett.119.051801
  * Ellis et al. (2016) Ellis, J., Garcia, M. A. G., Nanopoulos, D. V., Olive, K. A., & Peloso, M. 2016, JCAP, 03, 008, doi: 10.1088/1475-7516/2016/03/008
  * Ellis et al. (1992) Ellis, J. R., Gelmini, G., Lopez, J. L., Nanopoulos, D. V., & Sarkar, S. 1992, Nucl. Phys. B, 373, 399, doi: 10.1016/0550-3213(92)90438-H
  * Ellis et al. (1990) Ellis, J. R., Lopez, J. L., & Nanopoulos, D. V. 1990, Phys. Lett. B, 247, 257, doi: 10.1016/0370-2693(90)90893-B
  * Evoli et al. (2008) Evoli, C., Gaggero, D., Grasso, D., & Maccione, L. 2008, Nuclear Instruments and Methods in Physics Research A, 588, 22, doi: 10.1016/j.nima.2008.01.007
  * Ferriere (2001) Ferriere, K. M. 2001, Rev. Mod. Phys., 73, 1031, doi: 10.1103/RevModPhys.73.1031
  * Fichtel et al. (1972) Fichtel, C. E., Hartman, R. C., Kniffen, D. A., & Sommer, M. 1972, ApJ, 171, 31, doi: 10.1086/151254
  * Fichtel et al. (1975) Fichtel, C. E., Hartman, R. C., Kniffen, D. A., et al. 1975, ApJ, 198, 163, doi: 10.1086/153590
  * Fichtel & Kniffen (1984) Fichtel, C. E., & Kniffen, D. A. 1984, A&A, 134, 13
  * Fomin et al. (2017) Fomin, Y. A., Kalmykov, N. N., Karpikov, I. S., et al. 2017, Phys. Rev. D, 95, 123011, doi: 10.1103/PhysRevD.95.123011
  * Garny et al. (2016) Garny, M., Sandora, M., & Sloth, M. S. 2016, Phys. Rev. Lett., 116, 101302, doi: 10.1103/PhysRevLett.116.101302
  * Guépin et al. (2021) Guépin, C., Aloisio, R., Anchordoqui, L. A., et al. 2021, Phys. Rev. D, 104, 083002, doi: 10.1103/PhysRevD.104.083002
  * Hillas (2005) Hillas, A. M. 2005, J. Phys. G, 31, R95, doi: 10.1088/0954-3899/31/5/R02
  * Hunter et al. (1997) Hunter, S. D., et al. 1997, Astrophys. J., 481, 205, doi: 10.1086/304012
  * Ishiwata et al. (2020) Ishiwata, K., Macias, O., Ando, S., & Arimoto, M. 2020, JCAP, 01, 003, doi: 10.1088/1475-7516/2020/01/003
  * Jóhannesson et al. (2018) Jóhannesson, G., Porter, T. A., & Moskalenko, I. V. 2018, Astrophys. J., 856, 45, doi: 10.3847/1538-4357/aab26e
  * Kaapa et al. (2021) Kaapa, A., Kampert, K.-H., & Mayotte, E. 2021, PoS, ICRC2021, xxx, doi: 10.22323/1.395.0004
  * Kalashev & Kuznetsov (2016) Kalashev, O. K., & Kuznetsov, M. Y. 2016, Phys. Rev. D, 94, 063535, doi: 10.1103/PhysRevD.94.063535
  * Kampert et al. (2011) Kampert, K.-H., et al. 2011, Proc. 32nd International Cosmic Ray Conference (Beijing, China), doi: 10.7529/ICRC2011/V02/1087
  * Kaneta et al. (2019) Kaneta, K., Mambrini, Y., & Olive, K. A. 2019, Phys. Rev. D, 99, 063508, doi: 10.1103/PhysRevD.99.063508
  * Knurenko et al. (2013) Knurenko, S. P., Petrov, Z. E., Sidorov, R., et al. 2013, Proc. of 33rd ICRC.  https://arxiv.org/abs/1310.1978
  * Kuzmin & Rubakov (1998) Kuzmin, V., & Rubakov, V. 1998, Phys. Atom. Nucl., 61, 1028.  https://arxiv.org/abs/astro-ph/9709187
  * Lipari & Vernetto (2018) Lipari, P., & Vernetto, S. 2018, Phys. Rev. D, 98, 043003, doi: 10.1103/PhysRevD.98.043003
  * Mambrini & Olive (2021) Mambrini, Y., & Olive, K. A. 2021, Phys. Rev. D, 103, 115009, doi: 10.1103/PhysRevD.103.115009
  * Mayer-Hasselwander et al. (1980) Mayer-Hasselwander, H. A., Buccheri, R., Kanbach, G., et al. 1980, in Ninth Texas Symposium on Relativistic Astrophysics, Vol. 336, 211–222, doi: 10.1111/j.1749-6632.1980.tb15931.x
  * Mayer-Hasselwander et al. (1982) Mayer-Hasselwander, H. A., Bennett, K., Bignami, G. F., et al. 1982, A&A, 105, 164
  * Medina-Tanco & Watson (1999) Medina-Tanco, G. A., & Watson, A. A. 1999.  https://arxiv.org/abs/astro-ph/9905240
  * Narayanan et al. (2008) Narayanan, G., Heyer, M. H., Brunt, C., et al. 2008, Astrophys. J. Suppl., 177, 341, doi: 10.1086/587786
  * Navarro et al. (1996) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, Astrophys. J., 462, 563, doi: 10.1086/177173
  * Ostapchenko (2013) Ostapchenko, S. 2013, in EPJ Web of Conferences, Vol. 52, EDP Sciences, 02001
  * Peters (1961) Peters, B. 1961, Il Nuovo Cimento, 22, 800, doi: 10.1007/BF02783106
  * Pierog et al. (2015) Pierog, T., Karpenko, I., Katzy, J. M., Yatsenko, E., & Werner, K. 2015, Phys. Rev. C, 92, 034906, doi: 10.1103/PhysRevC.92.034906
  * Riehn et al. (2017) Riehn, F., Dembinski, H. P., Engel, R., et al. 2017, PoS, ICRC2017, 301, doi: 10.22323/1.301.0301
  * Savina & Pierre Auger Collaboration (2021) Savina, P., & Pierre Auger Collaboration. 2021, Proc. 37th International Cosmic Ray Conference (Berlin, Germany), PoS (ICRC2021), 373, doi: 10.22323/1.395.0373
  * Siffert et al. (2007) Siffert, B. B., Lazarotto, B., de Mello Neto, J. R. T., & Olinto, A. 2007, Braz. J. Phys., 37, 48, doi: 10.1590/S0103-97332007000100016
  * Taylor et al. (2015) Taylor, A. M., Ahlers, M., & Hooper, D. 2015, Phys. Rev. D, 92, 063011, doi: 10.1103/PhysRevD.92.063011
  * Ulrich et al. (2021) Ulrich, R., Pierog, T., & Baus, C. 2021, doi: 10.5281/zenodo.4558706
  * Zirakashvili et al. (1998) Zirakashvili, V. N., Pochepkin, D. N., Ptuskin, V. S., & Rogovaya, S. I. 1998, Astronomy Letters, 24, 139