# Where are Milky Way’s Hadronic PeVatrons?

Takahiro Sudoh<EMAIL_ADDRESS>Center for Cosmology and AstroParticle Physics
(CCAPP), Ohio State University, Columbus, OH 43210, USA  Department of
Physics, Ohio State University, Columbus, OH 43210, USA  Department of
Astronomy, Ohio State University, Columbus, OH 43210, USA  John F. Beacom
<EMAIL_ADDRESS>Center for Cosmology and AstroParticle Physics (CCAPP), Ohio
State University, Columbus, OH 43210, USA  Department of Physics, Ohio State
University, Columbus, OH 43210, USA  Department of Astronomy, Ohio State
University, Columbus, OH 43210, USA

###### Abstract

Observations of the Milky Way at TeV–PeV energies reveal a bright diffuse flux
of hadronic cosmic rays and also bright point sources of gamma rays. If the
gamma-ray sources are hadronic cosmic-ray accelerators, then they must also be
neutrino sources. However, no neutrino sources have been detected. Where are
they? We introduce a new population-based approach to probe Milky Way hadronic
PeVatrons, demanding consistency between diffuse and point-source PeV-range
data on cosmic rays, gamma rays, and neutrinos. For the PeVatrons, two extreme
scenarios are allowed: (1) the hadronic cosmic-ray accelerators and the gamma-
ray sources are the same objects, so that bright neutrino sources exist and
improved telescopes can detect them, versus (2) the hadronic cosmic-ray
accelerators and the gamma-ray sources are distinct, so that there are no
detectable neutrino sources. The latter case is possible if hadronic
accelerators have sufficiently thin column densities. We quantify present
constraints and future prospects, showing how to reveal the nature of the
hadronic PeVatrons.

##  I Introduction

The Milky Way is blazing in the GeV–PeV diffuse emission of hadronic cosmic
rays (CRs) Ginzburg and Syrovatskii (1964); Berezinskii  _et al._ (1990);
Gaisser (1990); Gaisser  _et al._ (2016). This requires the existence of
powerful accelerators, though their locations are obscured by magnetic
deflections during CR propagation. While many source classes seem able to
accelerate hadrons to GeV energies, it remains mysterious which ones can reach
PeV energies Lagage and Cesarsky (1983); Bell (2004); Antoni  _et al._ (2005);
Bell  _et al._ (2013); Blasi (2013); Anchordoqui  _et al._ (2014); Cristofari
(2021) — i.e., the hadronic PeVatrons, which are our focus. In principle,
these should emit bright fluxes of gamma rays and neutrinos due to the pion-
producing CR interactions with source matter and photons Stecker (1971);
Aharonian (2004); Dermer and Menon (2009).

The Milky Way is also blazing in the GeV–PeV point-source emission of gamma
rays Abdo  _et al._ (2007); Abeysekara  _et al._ (2017a); Abdalla  _et al._
(2018a); Abdollahi  _et al._ (2020); Cao  _et al._ (2021a). As the detector
energy range is increased, sources become rarer, but they are still found,
which indicates their powerful emission at high energies. Recently, the Large
High Altitude Air Shower Observatory (LHAASO) detected twelve northern-sky
sources up to 1.4 PeV Cao  _et al._ (2021a). It is often assumed that the
highest-energy gamma-ray sources are hadronic CR accelerators. This may be
true. But it might instead be true that these sources are only leptonic CR
accelerators, which produce gamma rays — but not neutrinos — through the
inverse-Compton scattering of CR electrons with source photons Atoyan and
Aharonian (1996).

How can we identify the Milky Way’s hadronic PeVatrons? Neutrino emission
would be a smoking gun. However, no sources have been found in more than a
decade of searches Ahlers and Halzen (2018); Aartsen  _et al._ (2019a, 2020);
Albert  _et al._ (2020a); Kheirandish (2020) despite a series of optimistic
predictions Kistler and Beacom (2006); Beacom and Kistler (2007); Halzen  _et
al._ (2008); Gonzalez-Garcia  _et al._ (2009). Figure 1 schematically
illustrates the range of possibilities. In the optimistic case, the highest-
energy gamma-ray sources are hadronic CR accelerators that produce neutrinos
and electrons, making them exciting multi-messenger sources. Then bright
neutrino sources exist and will be found with future detectors Aartsen  _et
al._ (2021); Adrián-Martínez  _et al._ (2016); Avrorin  _et al._ (2018);
Agostini  _et al._ (2020); Ye  _et al._ (2022). In the pessimistic case, the
highest-energy gamma ray sources are all leptonic CR accelerators. However,
this seems to be at odds with the requirement that some sources must
accelerate hadronic CRs. Here we point out a viable possibility that has
received little attention: the hadronic CR accelerators may be so _thin_ in
column density that the CRs escape without interaction. Then bright neutrino
sources would not exist and identifying the hadronic CR accelerators would be
difficult.

Figure 1: Possible scenarios for PeV sources. Each diamond is a source,
emitting the noted particles. Bottom left shows the “optimistic” case, where
all gamma-ray sources are hadronic PeVatrons with neutrino emission. Bottom
right shows the “pessimistic” case, where all gamma-ray sources are leptonic
and no hadronic PeVatrons produce detectable gamma-ray or neutrino fluxes.

We distinguish the last possibility from so-called hidden sources, which are
hadronic CR accelerators that are so _thick_ that gamma rays cannot escape,
though neutrinos can. But then these sources are not what we call hadronic
PeVatrons, because in many cases the CRs cannot escape either, leaving
unanswered where the Milky Way hadronic CRs come from. Hidden sources are
interesting, but they are not our focus.

In this paper, we present a new theoretical framework to systematically study
Milky Way hadronic PeVatrons, introducing three innovations. First, contrary
to previous studies on individual candidate sources, here we take a
population-studies approach. Second, we also take a multi-messenger approach,
demanding consistency with diffuse and point-source data on CRs, gamma rays,
and neutrinos. Third, we quantify the properties of PeVatrons in a semi-model-
independent plane of source gas density ($n_{\rm src}^{\rm gas}$) and CR
escape time ($\tau_{\rm src}^{\rm esc}$). As the inputs are uncertain over
orders of magnitude, we aim for a precision of a factor of a few, which we
show is adequate given present flux sensitivities. Our goal is to guide the
next steps as observations and theory improve. Our new theoretical framework
is a valuable complement to other approaches.

In Sec. II, we review very high energy astronomy. In Sec. III, we describe our
models for the hadronic source population in a simplified way, followed in
Sec. IV by the full details. In Sec. V, we calculate constraints based on
current observations. In Sec. VI, we calculate prospects for future neutrino
telescopes and discuss ways to make progress. In Sec. VII, we summarize key
new insights.

##  II Overview of Very High Energy Emission

In this section, we review observations and theory for very high energy
astrophysics. We start with diffuse hadronic CR emission, then point sources
of gamma rays and neutrinos, and finally diffuse gamma-ray and neutrino
emission.

The CR spectrum is dominated by protons. The observed CR proton intensity
Panov  _et al._ (2009); Adriani  _et al._ (2011); Aguilar  _et al._ (2015);
Yoon  _et al._ (2017); Atkin  _et al._ (2018); An  _et al._ (2019); Adriani
_et al._ (2019) is

$E_{p}^{2}\Phi_{p}\simeq 7\times 10^{-5}\left(\frac{E_{p}}{\rm
PeV}\right)^{-0.7}{\rm GeV~{}cm^{-2}~{}s^{-1}~{}sr^{-1}}.$ (1)

Near the PeV range, helium and other heavy components are also important
Höorandel (2006); Albert  _et al._ (2022); their contribution to hadronic
gamma-ray emission is discussed below. The all-particle CR spectrum has a
break at about 3 PeV, commonly defined as the “knee” Amenomori  _et al._
(2008); Apel  _et al._ (2013); Aartsen  _et al._ (2019b); Lipari and Vernetto
(2020); CRs below the knee are believed to be of a Galactic origin. Our
baseline scenario for the CR sources is supernova remnants (SNRs), which are
the most promising candidate up to the knee and possibly beyond Ptuskin  _et
al._ (2010), although decisive conclusions are yet to be reached Ahnen  _et
al._ (2017); Abeysekara  _et al._ (2020a). For population approaches based on
specific SNR models, see Refs. Cristofari  _et al._ (2013, 2018). We strive to
make our constraints general, so that they also cover other scenarios Bykov
(2014); Fujita  _et al._ (2017); Ioka  _et al._ (2017); Guépin  _et al._
(2018); Ohira  _et al._ (2018); Kimura  _et al._ (2018); Sudoh  _et al._
(2020); Cooper  _et al._ (2020); Kimura  _et al._ (2020); Guépin  _et al._
(2020); Morlino  _et al._ (2021); Kimura  _et al._ (2021).

An important observable is the grammage that CRs accumulate before escaping
the Milky Way, which is constrained by data on the secondary-to-primary ratio
Engelmann  _et al._ (1990); Blum  _et al._ (2013):

$X\simeq 8.7~{}{\rm g~{}cm^{-2}}\left(\frac{E_{p}}{\rm
10~{}GeV}\right)^{-\delta}.$ (2)

Above 65 GeV/nucleon, observations are consistent with a power law of
$\delta=0.33$ (at lower energies, the index is somewhat larger), but detailed
constraints are available only below 3 TeV/nucleon, inducing uncertainties
Aguilar  _et al._ (2016). Nevertheless, we assume that the scaling
$X\propto(E_{p})^{-0.33}$ holds up to the knee, as is expected for CRs that
scatter with external turbulence (as opposed to turbulence generated by CRs
themselves) Blasi  _et al._ (2012); Aloisio and Blasi (2013).

The Milky Way CR production rate can be estimated from the observed flux as

$E_{p}^{2}\frac{d^{2}N_{p}}{dE_{p}dt}\biggr{\rvert}_{\rm
MW}=\frac{4\pi}{c}E_{p}^{2}\Phi_{p}\frac{V_{\rm CR}}{\tau^{\rm esc}(E_{p})},$
(3)

where $V_{\rm CR}$ is the Galactic volume that CRs are confined within and
$\tau^{\rm esc}$ is the escape time from it. (Note that here we focus on the
Milky Way as a whole, whereas in most of the paper we focus on sources.) This
equation is derived assuming steady-state conditions, in which case the CR
energy density is the production rate on the left-hand side times $\tau^{\rm
esc}/V_{\rm CR}$. In the standard paradigm, it is assumed that the grammage,
at least in the GeV range, is dominated by the diffuse interstellar medium
(ISM) Moskalenko and Strong (1998); Strong  _et al._ (2007), (see, however,
Refs. Cowsik and Burch (2010); Lipari (2017, 2019); Gabici  _et al._ (2019)).
Then, the factor $V_{\rm CR}/\tau^{\rm esc}$ is calculated via $X=\mu
m_{p}n_{\rm ISM}^{\rm gas}c\tau^{\rm esc}=(M_{\rm gas}/V_{\rm CR})c\tau^{\rm
esc}$, where $m_{p}$ is the proton mass, $\mu\simeq 1.4$ accounts for the
composition of the gas, $n_{\rm ISM}^{\rm gas}$ is the gas number density, and
$M_{\rm gas}$ is the Milky Way’s gas mass (set to $10^{10}M_{\odot}$ Murase
and Fukugita (2019)). The energy-dependent CR production rate is then

$E_{p}^{2}\frac{d^{2}N_{p}}{dE_{p}dt}\biggr{\rvert}_{\rm MW}\sim 1.3\times
10^{38}\left(\frac{E_{p}}{\rm PeV}\right)^{-0.37}~{}{\rm erg~{}s^{-1}},$ (4)

which defines the energy budget of hadronic PeV sources Murase and Fukugita
(2019). Compared to widely assumed $E^{-2}$ (or $E^{-2.2}$), the above
spectral index is softer, yielding less power in the PeV range for a fixed
total energy. We focus on the flux in a narrow energy range near 1 PeV, which
makes our results robust to spectrum uncertainties. The CR spectrum from
PeVatrons might be harder than -2.37; the observed power-law spectrum can be
still obtained with additional soft-spectrum sources with lower maximum
energies. The typical spectral index for PeVatrons is uncertain, but
unimportant here, as our analysis is independent of the lower-energy emission.

The major uncertainties in Eq. (4) are the observed intensity and CR grammage.
First, the observed proton intensities differ between experiments by a factor
of about 1.6 at 1 PeV, which increases to 2.5 at 3 PeV Lipari and Vernetto
(2020). Eq. (1) uses the upper end of the range. The actual production rate
near 1 PeV might be smaller by a factor of 3. At higher energies,
discrepancies between measurements are even larger (a factor of 6 at 10 PeV).
Although such high energies are not our focus, in Sec. VI.2 we comment on how
future LHAASO gamma-ray data might help to resolve this tension. Second, we
extrapolate $X$ from 3 TeV to 1 PeV assuming $\delta=0.33$, which is predicted
by the Kolmogorov theory. Had we used $\delta=0.5$, as predicted by Kraichnan
theory, the value of $X$ would be a factor of (1 PeV / 3 TeV)0.17 smaller,
which means that the actual production rate might be larger by a factor of 3.
For the overall uncertainties on the production rate near $\sim$1 PeV, we thus
expect a factor less than $\sim 3$.

At the sources or during propagation, CR hadrons collide with gas, producing
roughly equal numbers of $\pi^{+}$, $\pi^{-}$, and $\pi^{0}$. The decay of a
$\pi^{0}$ meson generates two gamma rays, each with energy $E_{\gamma}\sim
0.1E_{p}$, while the decays of $\pi^{+}$ and $\pi^{-}$ (and the subsequent
decays of $\mu^{+}$ and $\mu^{-}$) produce three neutrinos, each with an
energy of $E_{\nu}\sim 0.05E_{p}$ (see Refs. Kelner  _et al._ (2006); Kafexhiu
_et al._ (2014); Celli  _et al._ (2020) for details). Hadronic gamma rays must
thus be accompanied by neutrinos with comparable numbers and energies.

Gamma-ray observatories have detected a large number and large variety of
Galactic sources in the GeV–PeV range. The search for the sources of
$E_{p}\sim 1$ PeV protons is most effective for secondaries around
$E_{\gamma}\sim 100~{}$TeV. Gamma-ray sources at such high energies are now
detected. Tibet AS$\gamma$ and the High-Altitude Water Cherenkov Observatory
(HAWC) have detected several gamma-ray sources emitting above 100 TeV
Amenomori  _et al._ (2019); Abeysekara  _et al._ (2020b); Albert  _et al._
(2020b, 2021a); Amenomori  _et al._ (2021a, b); Abeysekara  _et al._ (2021);
Amenomori  _et al._ (2022). LHAASO has expanded the source count above 100
TeV, detecting twelve uniformly selected sources in Ref. Cao  _et al._
(2021a), which we use. (They also reported one more in Ref. Cao  _et al._
(2021b).) Incredibly, some are observed at energies beyond 1 PeV, which is in
some tension with data on the CR knee, as we discuss in Sec. VI.2. These
telescopes, combined with observations by imaging atmospheric Cherenkov
telescopes (IACTs), have discovered particularly strong candidates for
hadronic PeVatrons Albert  _et al._ (2020b, 2021a); Amenomori  _et al._
(2021a); Abdalla  _et al._ (2021); Xin  _et al._ (2019); Fang  _et al._
(2022).

On the contrary, no Galactic neutrino sources have been discovered. As
discussed below, this can be interpreted as meaning that the sensitivities are
insufficient. However, it might also indicate that a significant fraction of
100-TeV gamma-ray sources are leptonic and that the accelerators of PeV
hadrons are so thin that they produce little gamma-ray and neutrino emission.
This suspicion is strengthened by observations that indicate that young and
middle-aged pulsars have high efficiencies for producing TeV gamma rays via
leptonic processes Tanaka and Takahara (2011); Torres  _et al._ (2014);
Abeysekara  _et al._ (2017b); Linden  _et al._ (2017); Abdalla  _et al._
(2018b); Sudoh  _et al._ (2019); Di Mauro  _et al._ (2021); Albert  _et al._
(2021b); Sudoh  _et al._ (2021); Breuhaus  _et al._ (2022a); de Oña Wilhelmi
_et al._ (2022); Joshi  _et al._ (2022); Fiori  _et al._ (2022); Martin  _et
al._ (2022).

Measurements of diffuse gamma-ray and neutrino emission are also important for
understanding the hadronic source population. Recently, Tibet AS$\gamma$ made
the first observations of diffuse gamma rays near the 100-TeV range Amenomori
_et al._ (2021c). The flux may be originating from the ISM Lipari and Vernetto
(2018); De la Torre Luque  _et al._ (2022), where CRs lose a small fraction,
$f_{\pi,\rm ISM}$, of energy to pions, where

$f_{\pi,\rm ISM}\sim 10^{-2}\left(\frac{X_{\rm ISM}}{1~{}\rm
g~{}cm^{-2}}\right),$ (5)

which is energy-dependent. This is derived using $f_{\pi,\rm ISM}\sim\tau_{\rm
ISM}^{\rm esc}/\tau^{pp}_{\rm ISM}$ and $\tau^{pp}_{\rm
ISM}=(\kappa_{pp}n_{\rm ISM}^{\rm gas}\sigma_{pp}c)^{-1}$, where
$\kappa_{pp}=0.5$ is the inelasticity and $\sigma_{pp}$ is the $pp$ cross
section, for which we use 50 mb Kelner  _et al._ (2006). We use $X_{\rm ISM}$
to explicitly note the contribution from the ISM. If we extrapolate Eq. (2),
the grammage for PeV particles is only $\sim 0.2$ g cm-2, implying that such
CRs lose only a small fraction of their energy, $\sim 2\times 10^{-3}$, to
pion production in the ISM.

Unresolved sources (both leptonic and hadronic) can contribute to the diffuse
flux Berezhko and Völk (2000); Prodanović  _et al._ (2007); Casanova and
Dingus (2008); Ahlers and Murase (2014); Ahlers  _et al._ (2016); Linden and
Buckman (2018); Fang and Murase (2021); Liu and Wang (2021); Vecchiotti  _et
al._ (2022). The grammage acquired in a source is

$X_{\rm src}=0.1\left(\frac{n_{\rm src}^{\rm gas}}{5~{}\rm
cm^{-3}}\right)\left(\frac{\tau_{\rm src}^{\rm esc}}{\rm
10~{}kyr}\right)\rm~{}g~{}cm^{-2},$ (6)

where $n_{\rm src}^{\rm gas}$ is the gas density at the source and $\tau_{\rm
src}^{\rm esc}$ is the escape time (i.e., how long CRs are confined by the
source.) The fraction of proton energy lost in sources can be estimated from
$X_{\rm src}$, similar to Eq. (5).

Because source spectra are often harder than the Galactic spectrum, the total
source emission might be expected to dominate the diffuse emission at the
highest energies. For source emission to be important, $X_{\rm src}$ must be
comparable to $X_{\rm ISM}$. Such a possibility has been discussed in the
GeV–TeV range, but not been constrained in the PeV range. Below, we show that
this is constrained by source counts.

##  III Models of the Hadronic Source Population: Simplified Overview

In this section, we introduce a model of the hadronic source population in
which the results are _explicitly normalized by the energy-dependent hadronic
CR energy budget in Eq. ( 4)_. We predict the expected source counts and their
gamma-ray and neutrino properties (e.g., luminosities, positions, and fluxes)
as well as the total emission from sources. Here we clarify the basic ideas
for the $n_{\rm src}^{\rm gas}$–$\tau_{\rm src}^{\rm esc}$ plane (hereafter
“$n$–$\tau$ plane”) introduced in this paper.

###  III.1 Impulsive Injection

We first consider the injection of CRs by impulsive events, as expected for
SNRs; in the subsequent subsection, we discuss the case where injection is
continuous. The total energy produced in CRs per event follows from Eq. (4)
and the rate of events that produce them (“CR source rate”), $\Gamma_{\rm
CR}$, as

$E_{p}^{2}\frac{dN_{p}}{dE_{p}}\biggr{\rvert}_{I}=\frac{1}{\Gamma_{\rm
CR}}E_{p}^{2}\frac{d^{2}N_{p}}{dE_{p}dt}\biggr{\rvert}_{\rm MW},$ (7)

where $I$ stands for impulsive injection. In the baseline scenario,
$\Gamma_{\rm CR}$ is set by the SN rate (0.03 yr-1), though it might be that
only a subclass of SN explosions produces CRs that reach PeV energies Murase
_et al._ (2014); Cristofari  _et al._ (2020). We discuss other choices below.

The gamma-ray and neutrino emission from these sources depends on how dense
the target is ($n_{\rm src}^{\rm gas}$) and how long CRs are confined
($\tau_{\rm src}^{\rm esc}$). We assume that the CR energy is impulsively
injected, i.e., in a timescale $\Delta t$ shorter than $\tau_{\rm src}^{\rm
esc}$ (although Eq. (7) is a good approximation even if $\Delta t\sim\tau_{\rm
src}^{\rm esc}$). We thus take the confined CR energy to be constant, as
described by Eq. (7), during the time $\tau_{\rm src}^{\rm esc}$ and zero
afterwards, when the source CRs mix in with Galactic CRs. The total number of
sources is then $\Gamma_{\rm CR}\tau_{\rm src}^{\rm esc}$.

The gamma-ray and neutrino ($\nu+\bar{\nu}$, per flavor) luminosities of a
source, normalized to CR data as per Eq. (7), are

$E_{\gamma}^{2}\frac{d^{2}N_{\gamma}}{dE_{\gamma}dt}=\frac{\epsilon_{A}}{3\tau^{pp}_{\rm
src}}E_{p}^{2}\frac{dN_{p}}{dE_{p}}\biggr{\rvert}_{I,\ E_{p}=10E_{\gamma}},$
(8)
$E_{\nu}^{2}\frac{d^{2}N_{\nu_{i}}}{dE_{\nu}dt}=\frac{\epsilon_{A}}{6\tau^{pp}_{\rm
src}}E_{p}^{2}\frac{dN_{p}}{dE_{p}}\biggr{\rvert}_{I,\ E_{p}=20E_{\nu}},$ (9)

where $\epsilon_{A}$ takes into account the fact that both the CR and the
target material contain heavy nuclei (mostly helium) Stecker (1970); Gaisser
and Schaefer (1992); Mori (2009); Kachelriess  _et al._ (2014). We use
$\epsilon_{A}=2.6$, as found by Ref. Peron and Aharonian (2022) for the PeV
range, which depends on the uncertain CR composition; see also Ref. Breuhaus
_et al._ (2022b). We refer to them as luminosities _near_ specific energies,
which means integrated over bins with $\Delta\ln E\sim 1$. We do not include
gamma-ray attenuation; at 100 TeV, the effect is at most a factor of $\sim$2
and certainly smaller for most sources Moskalenko  _et al._ (2006), but note
that this becomes increasingly important at higher energies.

Figure 2: Schematic constraints in the $n$–$\tau$ plane for PeV sources. Panel
A considers hypothetical gamma-ray source data in light of real CR data. If
all gamma-ray sources are hadronic, then the sources must lie within the blue
band; if none are, then the sources must lie below the blue dotted line.
Regardless of the interpretation of the nature of the gamma-ray sources, the
upper right of the plane is excluded. Panels B and C compare hypothetical
neutrino source data (green bands/lines, with similar interpretations). The
dashed gray line indicates where the column density encountered by CRs in the
source ($n_{\rm src}^{\rm gas}c\tau_{\rm src}^{\rm esc}$) equals that
encountered in the ISM before escape from the Milky Way ($X_{\rm ISM}$), a
relevant scale for comparison. Stars indicate possible source classes. For
this figure, we assume a CR source rate of $0.03$ yr-1. This new, semi-model-
independent plane allows consistent comparisons of CR, gamma-ray, and neutrino
data. See text for details.

Quantitatively, we find the total number of gamma-ray and neutrino sources in
the Milky Way to be :

$\mathcal{N}_{\gamma}=\mathcal{N}_{\nu}\sim 300\left(\frac{\Gamma_{\rm
CR}}{0.03~{}\rm yr^{-1}}\right)\left(\frac{\tau_{\rm src}^{\rm esc}}{10~{}\rm
kyr}\right).$ (10)

These sources have luminosities

$\begin{split}E_{\gamma}^{2}\frac{d^{2}N_{\gamma}}{dE_{\gamma}dt}\biggr{\rvert}_{E_{\gamma}=100~{}\rm
TeV}&\sim 10^{32}~{}{\rm erg~{}s^{-1}}\\\ &\times\left(\frac{\Gamma_{\rm
CR}}{0.03~{}\rm yr^{-1}}\right)^{-1}\left(\frac{n_{\rm src}^{\rm gas}}{\rm
cm^{-3}}\right),\end{split}$ (11)
$\begin{split}E_{\nu}^{2}\frac{d^{2}N_{\nu_{i}}}{dE_{\nu}dt}\biggr{\rvert}_{E_{\nu}=50~{}\rm
TeV}&\sim 5\times 10^{31}~{}{\rm erg~{}s^{-1}}\\\
&\times\left(\frac{\Gamma_{\rm CR}}{0.03~{}\rm
yr^{-1}}\right)^{-1}\left(\frac{n_{\rm src}^{\rm gas}}{\rm
cm^{-3}}\right).\end{split}$ (12)

The details of the spectrum shape are not important here, as we only focus on
a range of energies near $E_{p}\sim 1$ PeV. For comparison, the gamma-ray
luminosity of the Crab Nebula at 100 TeV is $\sim 5\times 10^{32}$ erg s-1 Cao
_et al._ (2021c).

The detection horizon depends on the experimental sensitivity, $F_{\rm lim}$,
as

$d_{\rm lim}=\sqrt{\frac{1}{4\pi F_{\rm
lim}}E_{i}^{2}\frac{d^{2}N_{i}}{dE_{i}dt}},$ (13)

where $i=\gamma$ or $\nu$. A fraction $\sim\xi_{\rm f}(d_{\rm lim}/R_{\rm
gal})^{2}$ of the total source population can be detected, where we assume
that we are at the center of a thin disk and the factor $\xi_{\rm f}$ accounts
for the limited fraction of the sky (we take $\xi_{\rm f}\sim 1/3$). The
expected gamma-ray source count at 100 TeV is then

$\mathcal{N}_{\rm det}\sim 3\left(\frac{\tau_{\rm src}^{\rm esc}\,n_{\rm
src}^{\rm gas}}{10~{}\rm kyr~{}cm^{-3}}\right)\left(\frac{F_{\rm
lim}}{10^{-13}~{}\rm erg~{}cm^{-2}~{}s^{-1}}\right)^{-1},$ (14)

where we take the Galactic radius to be $R_{\rm gal}=$15 kpc.

Figure 2 shows the concepts behind using the $n$–$\tau$ plane to constrain
hadronic PeVatrons. Here we illustrate a simple version of the figure with
hypothetical data, explaining the implications one by one. In the next
section, we use real data. We start by focusing on $n_{\rm src}^{\rm gas}$ and
$\tau_{\rm src}^{\rm esc}$, and then discuss the effects of changing the CR
source rate, $\Gamma_{\rm CR}$, and the source size, $R_{\rm src}$.

  * •

_Cosmic-ray data:_ By construction, any source class in this plane is
consistent with the energy-dependent CR production rate in the PeV range, as
in Eq. (4). Hadronic sources that do not reach PeV energies, or those that
accelerate PeV particles but do not allow escape (CR hidden sources), are not
considered. To derive the CR yield per source, as in Eq. (7), we assume a CR
source rate of $\Gamma_{\rm CR}=0.03$ yr-1 (alternatives are discussed below).

  * •

_Gamma-ray data:_ Consider a hypothetical gamma-ray observation that attains
$F_{\rm lim}=10^{-13}~{}\rm erg~{}cm^{-2}~{}s^{-1}$ and observes 10 sources at
$E_{\gamma}=100$ TeV. If all of them are hadronic PeVatrons, Eq. (14)
indicates that $\tau_{\rm src}^{\rm esc}\,n_{\rm src}^{\rm gas}\sim 30~{}\rm
kyr~{}cm^{-3}$, with smaller values of $\tau_{\rm src}^{\rm esc}\,n_{\rm
src}^{\rm gas}$ allowed if there is a smaller fraction of hadronic sources
(and hence a larger fraction of leptonic sources). In panel A, we show two
extreme cases, where 100$\%$ or 0$\%$ of them are hadronic, corresponding to
the optimistic and pessimistic cases in Fig. 1. Gamma-ray data alone are
usually insufficient to decisively determine if a source is hadronic or not.

If all of the observed gamma-ray sources are hadronic, then the parameters
$n_{\rm src}^{\rm gas}$ and $\tau_{\rm src}^{\rm esc}$ must be within the
shaded area enclosed by the solid band labeled “Required if 100$\%$ hadronic.”
This area is defined by $\mathcal{N}_{\rm det}=10$, with the width indicating
statistical fluctuations. If none of the observed sources are hadronic,
$n_{\rm src}^{\rm gas}$ and $\tau_{\rm src}^{\rm esc}$ must be below the
dotted line labeled “Upper limit if 0$\%$ hadronic,” obtained by setting
$\mathcal{N}_{\rm det}<1$ (ignoring statistical fluctuations). One could also
display parameter spaces that are consistent with intermediate cases (e.g.,
half of them being hadronic).

Regardless of the hadronic fraction of gamma-ray sources, the upper right side
of the plane is _always_ ruled out, due to predicting more sources than
observed.

While parameters within the band predict the same _source counts_ , the
_luminosities_ of gamma-ray sources vary along the band. In the upper left,
sources are common but less luminous, making the observable sources nearby.
Towards the middle, sources become rarer but more luminous, making the
observable sources more distant. Eventually, $\tau_{\rm src}^{\rm esc}$
becomes smaller than $10/(\xi_{f}\Gamma_{\rm CR})$, at which point the gamma-
ray source number in the field of view becomes fewer than ten. Then it becomes
impossible to attain $\mathcal{N}_{\rm det}=10$ for any $n_{\rm src}^{\rm
gas}$ by continuing the diagonal, which causes a flattening of the band. For
the same reason, the dotted lines flatten at $1/(\xi_{f}\Gamma_{\rm CR})$,
corresponding to one source.

  * •

_Neutrino data:_ Consider two hypothetical neutrino experiments, with
sensitivities of $F_{\rm lim}=10^{-11}$ and $2\times 10^{-13}~{}\rm
erg~{}cm^{-2}~{}s^{-1}$ at $E_{\nu}=50$ TeV, and suppose they detect no
sources. Non-detections define dotted lines, and $n_{\rm src}^{\rm gas}$ and
$\tau_{\rm src}^{\rm esc}$ must be below them, similar to the case of zero
hadronic sources for the gamma-ray data. Those cases are marked as “Upper
limit by $\nu$ (Weak)” and “Upper limit by $\nu$ (Strong)” in panel B of Fig.
2. As with the gamma-ray case, the upper right of the plane is ruled out.

Next, consider a hypothetical neutrino experiment with sensitivity of
$10^{-13}~{}\rm erg~{}cm^{-2}~{}s^{-1}$ at $E_{\nu}=50$ TeV, and suppose that
it finds one source, which must then be hadronic, defining a robust allowed
parameter space. This case is shown in panel C of Fig. 2, marked as “Required
by $\nu$,” where we set the neutrino source count to $\mathcal{N}_{\rm
det}=1$. We expect the area (statistical fluctuations) to be larger than the
gamma-ray case, due to the smaller source number. And as with the gamma-ray
case, the parameters must be in the band, meaning that the lower left of the
plane is also ruled out.

  * •

_Gamma-ray data vs. Neutrino data:_ Comparing gamma-ray and neutrino
constraints provides information on the fraction of hadronic sources in the
gamma-ray observations. In panel B, the “Neutrino (Weak)” limit would be
consistent with a scenario where 100$\%$ of gamma-ray sources are hadronic,
except for the bottom right region, where this is excluded. In contrast, the
“Neutrino (Strong)” limit would exclude all scenarios where 100$\%$ of gamma-
ray sources are hadronic.

In panel C, there is a narrow parameter space that is consistent with the
scenario where 100$\%$ of gamma-ray sources are hadronic, which is where the
two bands overlap. The neutrino data rule out the scenario where 0$\%$ of
gamma-ray sources are hadronic, as expected, because the dotted blue line is
outside the band allowed by the neutrino source detection. This could be true
even if the gamma-ray and neutrino telescopes viewed disjoint regions of the
sky.

  * •

_Data vs. Theory:_ With observational constraints on $n_{\rm src}^{\rm gas}$
and $\tau_{\rm src}^{\rm esc}$, we can compare these with expectations from
specific models of hadronic PeVatrons. In panels B and C, we show three
representative cases with star symbols (“theory points”), each explained
below.

If all gamma-ray sources are hadronic, then the theory points should be within
the allowed band. If none of the sources are hadronic, then the theory points
should be below the blue dotted line. Note that even if a theory point is in
the “100$\%$ hadronic” band, it does not always mean that this source class
must explain 100$\%$ of the gamma-ray sources, due to statistical
fluctuations.

  * •

_Source vs. Diffuse:_ The gray dashed lines show the grammage that CRs
accumulate in the ISM, or equivalently the fraction of energy lost, as per Eq.
(5), obtained by extrapolating Eq. (2) to 1 PeV. For points above this line,
CRs lose more energy at the source than in the ISM before they escape the
Galaxy, i.e., $X_{\rm src}>X_{\rm ISM}$. Above this line, the total Galactic
emission in gamma rays and neutrinos would be dominated by sources as opposed
to diffuse emission. Our full population models (next section) consistently
calculate the resolved source counts and also the contribution of sources to
the total Galactic emission of gamma rays and neutrinos.

The underlying idea of the $n$–$\tau$ plane is allow a consistent comparison
of different results — CR, gamma-ray, and neutrino data, for both source and
diffuse emission — to each other and to theoretical expectations.

As noted, we _assume_ a CR source rate of $\Gamma_{\rm CR}=0.03$ yr-1, and
then determine the CR energy per source from the observed CR data. If we adopt
different values of $\Gamma_{\rm CR}$, most constraints are unchanged, as the
source count, as in Eq. (14), does not depend on it. However, the flattening
of the band is shifted, proportionally to $(\Gamma_{\rm CR})^{-1}$; we show
this case below and in the Appendix.

For convenience, we convert the source number into the distance of the nearest
source. Assuming a two-dimensional geometry as above, it is $d_{\rm
near}=R_{\rm MW}/\sqrt{\mathcal{N}}$, or

$d_{\rm near}\sim 0.9~{}{\rm kpc}\left(\frac{\Gamma_{\rm CR}}{0.03~{}\rm
yr^{-1}}\right)^{-1/2}\left(\frac{\tau_{\rm src}^{\rm esc}}{10~{}\rm
kyr}\right)^{-1/2}.$ (15)

This is to be contrasted with the maximum distance telescopes can see, $d_{\rm
lim}$, which we note for each experiment in the next section.

An important effect not yet discussed is the source size, $R_{\rm src}$. In
the rest of the paper, we adopt a source size of 10 pc unless otherwise noted.
The sensitivities of gamma-ray and neutrino telescopes are degraded for
sources that are more extended than the size of the point-spread function
(PSF; $\theta_{\rm PSF}$) by a factor of $\sim\theta_{\rm src}/\theta_{\rm
PSF}$ Hinton and Hofmann (2009), where $\theta_{\rm src}$ is the source
angular extension. This quantitatively changes the band when the point-source
detection horizon in Eq. (13) is smaller than the distance at which
$\theta_{\rm src}=\theta_{\rm PSF}$ (which we denote as $d_{\rm size}$). The
detection horizon in this case shrinks to

$d_{\rm lim}=\frac{L}{4\pi F_{\rm lim}^{\rm PS}d_{\rm size}},$ (16)

where $F_{\rm lim}^{\rm PS}$ is the sensitivity for point sources. The
dependence of horizon distance on luminosities changes to $d_{\rm lim}\propto
L$, resulting in $\mathcal{N}_{\rm det}\propto(d_{\rm lim})^{2}\tau_{\rm
src}^{\rm esc}\propto L^{2}\tau_{\rm src}^{\rm esc}\propto(n_{\rm src}^{\rm
gas})^{2}\tau_{\rm src}^{\rm esc}$, so the source count does not depend
linearly on $X_{\rm src}$.

Having introduced the basics of this figure, the following questions should be
addressed:

  * •

_What sets the range of parameters plotted?_ While both $n_{\rm src}^{\rm
gas}$ and $\tau_{\rm src}^{\rm esc}$ can vary over many orders of magnitude,
we only show limited parameter ranges. We encourage new work to sharpen our
choices.

There are considerations that set the largest values that the parameters can
take. Sources can have $n_{\rm src}^{\rm gas}$ larger than displayed in Fig.
2. However, we show below that $n_{\rm src}^{\rm gas}\gtrsim 1000$ cm-3 is
firmly ruled out by existing IceCube data, except for very short escape time
($\lesssim$0.5 kyr). For $\tau_{\rm src}^{\rm esc}$, an upper bound comes from
a physical argument. Particle diffusion is slowest in the Bohm diffusion
regime, where the coefficient is $cR_{L}/3$, which is $3\times 10^{28}$ cm2
s-1 for a 1-PeV particle in a $\mu$G field, where $R_{L}$ is the Larmor
radius. If such a slow diffusion is sustained over 100-pc scale (note that
this is extremely optimistic), we obtain $\tau_{\rm src}^{\rm esc}\sim 100$
kyr. Such a large value, although cannot be excluded from first principles, is
certainly too high, though we display it in the figure to show present
constraints without theoretical bias.

There are also considerations that set the smallest values that the parameters
can take. No immediate physical arguments prohibit sources from having small
values of $n_{\rm src}^{\rm gas}$, though some matter is required to support
the magnetic fields that accelerate (and confine) CRs. The gas densities can
be much smaller than 1 cm-3 if, for example, a SN occurs in a cavity region
where material has been blown out due to stellar winds or SN explosions. For
$\tau_{\rm src}^{\rm esc}$, the absolute minimum would be the time needed to
accelerate particles to energies beyond 1 PeV, which in principle could be as
small as $\sim R_{L}/c$, which is a few years for a $\mu$G magnetic field and
is smaller for a stronger field Aharonian  _et al._ (2002). These small values
are not shown in the figure, but may be physically plausible.

  * •

_How should the complexity of sources be incorporated into phenomenological
descriptions with $n_{\rm src}^{\rm gas}$ and $\tau_{\rm src}^{\rm esc}$?_ The
environments of CR sources and the production sites of gamma rays should be
quite complicated. In the case of SNRs, gamma rays may be produced by either
diffuse matter or gas clumps inside the shell or in their close vicinity, with
all of these components likely being highly inhomogeneous Aharonian  _et al._
(1994); Bykov  _et al._ (2000); Uchiyama  _et al._ (2010); Zirakashvili and
Aharonian (2010); Caprioli (2011); Inoue  _et al._ (2012); Aharonian (2013);
Celli  _et al._ (2019a); Baade and Zwicky (1934); Sano  _et al._ (2020); Sano
and Fukui (2021); Yasuda  _et al._ (2021, 2022). Moreover, the accelerator of
particles and the emitter of gamma rays might be physically distinct. For
example, an SNR may accelerate protons, which can interact with nearby
molecular clouds (MCs), emitting gamma rays and neutrinos Aharonian and Atoyan
(1996); Gabici and Aharonian (2007); Aharonian  _et al._ (2008); Gabici  _et
al._ (2009); Fujita  _et al._ (2009); Ohira  _et al._ (2011); Uchiyama  _et
al._ (2012); Mitchell  _et al._ (2021).

For a single source class, we can separately consider multiple “emitting
regions,” each characterized by different combinations of $n_{\rm src}^{\rm
gas}$ and $\tau_{\rm src}^{\rm esc}$. In the case of SNRs, one can separately
consider SNR shells, wind bubbles, wind cavities, and nearby molecular clouds,
for example. Due to the likely inhomogeneities of the densities, the parameter
$n_{\rm src}^{\rm gas}$ should be regarded as an average gas density that PeV
CRs encounter. Formally, it should be calculated by convolving the gas density
with the spatial distribution of CRs, which requires a detailed understanding
of CR propagation and the gas distribution. More simply, we can separate the
source into sub-volumes of density $n_{v}$ and filling fraction $f_{v}$ and
expect $n_{\rm src}^{\rm gas}\sim\sum_{v}f_{v}n_{v}$, with $\sum_{v}f_{v}=1$.

  * •

_Where do realistic models lie in this plane?_ Below, we discuss three
concrete cases, each marked in Fig. 2. Although we focus on SNRs as
illustrative examples, our framework can be applied to more general source
classes, for which we encourage further work.

    1. 1.

SNR: The most plausible candidate for the source of PeV hadrons are SNRs.
Indeed, both GeV and TeV data support the scenario where SNRs accelerate
protons to TeV scales Morlino and Caprioli (2012); Ackermann  _et al._ (2013);
Jogler and Funk (2016); Abdalla  _et al._ (2018c), although it remains unknown
if the maximum energy can reach the PeV range. If produced, PeV protons are
expected to escape in the very early phase of the SNR evolution, $\tau_{\rm
src}^{\rm esc}\sim$ 1 kyr, comparable to when the Sedov-Taylor phase starts,
although many details are uncertain Drury (2011); D’Angelo  _et al._ (2016);
Nava  _et al._ (2016); Gaggero  _et al._ (2018); Nava  _et al._ (2019); Celli
_et al._ (2019b); Brose  _et al._ (2020, 2021); Kamijima and Ohira (2021);
Suzuki  _et al._ (2022). Acceleration to the PeV scale might take place on
timescales much shorter than kiloyears Ptuskin and Zirakashvili (2003); Schure
and Bell (2013); Cardillo  _et al._ (2015); Marcowith  _et al._ (2018); Inoue
_et al._ (2021) and $\tau_{\rm src}^{\rm esc}$ includes the time particles are
in the vicinity of the accelerator. This choice of $\tau_{\rm src}^{\rm esc}$
is likely optimistic, given the escape of PeV particles from the shock might
be as short as $\sim$10 yrs Schure and Bell (2013) and the escape from the
larger surroundings are highly uncertain. If PeV protons interact with the
average gas densities in the ISM, then $n_{\rm src}^{\rm gas}\sim$ 1 cm-3.

    2. 2.

SNR (Dense): For a handful of shell-type SNRs, TeV gamma rays are spatially
coincident with gas clouds, supporting a hadronic origin for the gamma rays:
RX 1713.7-3946 Aharonian  _et al._ (2007); Fukui  _et al._ (2012, 2021), Vela
Jr. Fukui  _et al._ (2017), HESS J1731-347 Fukuda  _et al._ (2014), and RCW86
Sano  _et al._ (2019). They are young ($\simeq 1$$-$$5$ kyr) and have high
target gas densities of $\gtrsim$10–100 cm-3, although the latter
significantly depends on the volume filling factor of dense gas, which is
usually highly uncertain. Here, we take RX J1713.7-3946 as an example case.
Multi-wavelength modeling of this SNR suggests high-density ($2.5\times
10^{4}$ cm-3) clumps with a volume filling factor $10^{-2}$, embedded in low-
density gas ($\sim 10^{-2}$ cm-3), for an average density of $n_{\rm src}^{\rm
gas}\sim$ 250 cm-3 in the 10-pc shell Fujita  _et al._ (2022). The escape of
PeV particles is model-dependent. Given the lack of $>$ 10 TeV gamma-ray
emission from this object, $\tau_{\rm src}^{\rm esc}$ is likely smaller than
its age (1.4 kyr). We optimistically take $\tau_{\rm src}^{\rm esc}=1$ kyr, as
in the previous case.

    3. 3.

MC near SNR: Emission may be produced by CRs that escape from the accelerators
and diffuse around them, interacting with a massive gas cloud or clouds. The
duration is determined by the local propagation of CRs. As a reference, we
consider a molecular cloud with a mass of $M_{\rm cl}=10^{5}~{}M_{\odot}$ and
a size of $R_{\rm cl}=20$ pc at a distance of $d_{\rm cl}=$ 50 pc from an SNR.
With a diffusion coefficient of $D=10^{29}~{}$cm2 s-1 (ten times smaller than
the ISM average), the diffusion time is $\tau_{\rm src}^{\rm esc}\sim(d_{\rm
cl})^{2}/D\sim$ 5 kyr. The gas density of this MC is very high, $\simeq
100~{}{\rm cm^{-3}}$, but the volume filling fraction of this is $\sim(R_{\rm
cl}/d_{\rm cl})^{3}\sim 0.06$, resulting in a modest value of $n_{\rm
src}^{\rm gas}\sim$ 6 cm-3.

  * •

_How should the variation among sources be treated?_ Both $n_{\rm src}^{\rm
gas}$ and $\tau_{\rm src}^{\rm esc}$ should vary among sources. In the case of
SNRs, the measured gas densities are known to vary by more than an order of
magnitude among different sources Dermer and Powale (2013); Zeng  _et al._
(2019); Suzuki  _et al._ (2020). There are two theoretical possibilities to
account for this source-to-source variation.

First, among all varieties of a source class, only a subclass with specific
$n_{\rm src}^{\rm gas}$ and/or $\tau_{\rm src}^{\rm esc}$ values might be able
to produce PeV hadronic CRs. In this case, an appropriate value of
$\Gamma_{\rm CR}$ should be chosen to take into account the lower rate. A
smaller value of $\Gamma_{\rm CR}$ does not move the “theory points,” but does
change the allowed bands and limit lines, as discussed above.

Second, sources with a wide range of $n_{\rm src}^{\rm gas}$ and/or $\tau_{\rm
src}^{\rm esc}$ values may produce PeV hadronic CRs, and only those with
sufficiently large values of these parameters may be detectable in gamma rays
and neutrinos. Formally, one should introduce a probability density
$d^{2}\mathcal{P}/dn_{\rm src}^{\rm gas}d\tau_{\rm src}^{\rm esc}$. Simply, we
may use the $n_{\rm src}^{\rm gas}$ and $\tau_{\rm src}^{\rm esc}$ values that
maximally contribute to the source counts. We can instead imagine a situation
where the distributions are approximated as bimodal. Such a scenario can be
studied by using Fig. 2 as a guide. For example, a fraction $f_{\rm dense}$ of
SNR may be in dense regions, like in the “SNR (Dense)” point, while the
remainder are in a normal environment, like in the “SNR” point. Figure 2 shows
that the “SNR (Dense)” and “SNR” models predict that $\sim 100\%$ and $\sim
0\%$ of gamma-ray sources are hadronic, respectively. Then, in total we expect
that $\sim 100f_{\rm dense}\%$ are hadronic. This consideration can be
extended to more general distributions. The key point is that the
applicability of Fig. 2 is wide, despite the fact that it assumes a universal
$n_{\rm src}^{\rm gas}$ and $\tau_{\rm src}^{\rm esc}$.

By construction, the $n$–$\tau$ plane aims to constrain the source class that
dominates PeV CR energy budget. It is quite plausible that other source
classes produce gamma-ray (or even neutrino) fluxes but make subdominant
contributions to the hadronic CR flux. For a source class that produce CRs at
a rate a factor of $\xi$ smaller than in Eq. (4), all the constrains would
approximately be shifted rightwards by a factor of $\xi$.

###  III.2 Continuous Injection

The case of continuous production is straightforwardly obtained from the
results above. The luminosity of each source follows from Eq. (4) and the
number of objects producing CRs, $\mathcal{N}_{\rm CR}$, as

$E_{p}^{2}\frac{d^{2}N_{p}}{dE_{p}dt}\biggr{\rvert}_{C}=\frac{1}{\mathcal{N}_{\rm
CR}}E_{p}^{2}\frac{d^{2}N_{p}}{dE_{p}dt}\biggr{\rvert}_{\rm MW},$ (17)

where $C$ stands for continuous injection. The CR energy spectrum is

$E_{p}^{2}\frac{dN_{p}}{dE_{p}}\biggr{\rvert}_{C}=\frac{\tau_{\rm src}^{\rm
esc}}{\mathcal{N}_{\rm
CR}}E_{p}^{2}\frac{d^{2}N_{p}}{dE_{p}dt}\biggr{\rvert}_{\rm MW}.$ (18)

Here, the sources are assumed to be in steady-state conditions, which
necessitates the loss time $\tau_{\rm src}^{\rm esc}$ to be less than the age
of the source (more than Myr-scales for the case of star clusters). Then Eqs.
(8) and (9) can be used to calculate luminosities, with $I$ replaced with $C$.

In the continuous case, the total Galactic source number $\mathcal{N}$ is
fixed, and larger values of $n_{\rm src}^{\rm gas}$ and $\tau_{\rm src}^{\rm
esc}$ both lead to larger luminosities, as $L_{\gamma}\propto n_{\rm src}^{\rm
gas}\tau_{\rm src}^{\rm esc}$. In the plane of $n_{\rm src}^{\rm
gas}$–$\tau_{\rm src}^{\rm esc}$, the upper right region corresponds to more
luminous sources. Because the detection area increases with $d_{\rm
lim}^{2}\propto L_{\gamma}$, the detectable source count increases accordingly
as $\mathcal{N}_{\rm det}\propto\mathcal{N}L_{\gamma}\propto n_{\rm src}^{\rm
gas}\tau_{\rm src}^{\rm esc}$. This is in contrast to the impulsive case,
where the total Galactic source number is $\mathcal{N}\propto\tau_{\rm
src}^{\rm esc}$ and their luminosities are $L_{\gamma}\propto n_{\rm src}^{\rm
gas}$, although this case also results in $\mathcal{N}_{\rm
det}\propto\mathcal{N}L_{\gamma}\propto n_{\rm src}^{\rm gas}\tau_{\rm
src}^{\rm esc}$. Additionally, the band does not flatten in the continuous
case.

An example case is stellar winds of massive stars in star clusters Bykov  _et
al._ (2020). Gamma-ray observations support the idea that this source class is
continuously injecting CRs into the ISM Aharonian  _et al._ (2019). The extent
of the gamma-ray emitting region is often greater than the size of the cluster
itself and can exceed $\sim$50 pc. The diffusion in this region is slow.
Assuming $D=10^{29}$ cm2 s-1 (ten times smaller than the ISM average),
$\tau_{\rm src}^{\rm esc}\sim 5$ kyr. The gas densities are suggested to be
high $n_{\rm src}^{\rm gas}\sim 10~{}$cm-3 Yang  _et al._ (2018); Sun  _et
al._ (2020); Abeysekara  _et al._ (2021). These points suggest that star
clusters are promising PeVatron candidates that can be probed by gamma rays
and neutrinos. (In compact star clusters or more loose young star
associations, multiple SNRs can also be produced, which might efficiently
accelerate hadrons compared to isolated SNRs Higdon and Lingenfelter (2005);
Gupta  _et al._ (2020); Bhadra  _et al._ (2022); Vieu  _et al._ (2022a, b),
which should be viewed as “impulsive” injection.)

##  IV Models of the Hadronic Source Population: Technical Methods

In this section, we describe the methods for our Monte-Carlo simulations of
the population of hadronic gamma-ray and neutrino sources, with the results
and interpretation given in the next section. We continue to use the
phenomenological parameters $\tau_{\rm src}^{\rm esc}$ and $n_{\rm src}^{\rm
gas}$ as in Eqs. (10), (11), and (12).

For each parameter set $\\{\tau_{\rm src}^{\rm esc},n_{\rm src}^{\rm
gas},R_{\rm src}\\}$, we run the simulation $10^{4}$ times, each time sampling
$\Gamma_{\rm CR}\tau_{\rm src}^{\rm esc}$ sources and calculating the
distributions of observables. We place sources randomly in the Milky Way
plane, following the radial distribution of initial positions of neutron stars
used in Ref. Sartore  _et al._ (2010) (the “F06 model” Faucher-Giguère and
Kaspi (2006)), and assuming a Gaussian distribution for the height above the
plane, with a standard deviation of 30 pc Faucher-Giguère and Kaspi (2006). We
then calculate the source luminosities with Eqs. (11) and (12). We include a
variety of corrections, including how the search sensitivities depend on
declination and the extent of the sources, as discussed above. We set the
maximum detectable angular size to be $\theta_{\rm max}=5^{\circ}$, which
excludes sources closer than 100 $(R_{\rm src}/10~{}\rm pc)$ pc.

In the next section, we compare several quantities with present data, while
here we only describe what we calculate and how. For predictions to agree with
observations, we require that the observed values be between the 10th and 90th
percentiles. To disagree, we require that the observed value be less than the
10th percentile. This is conservative because we consider a parameter set to
be ruled out only if more than 90$\%$ of the predicted range disagrees with
observations.

  * •

_Source counts from LHAASO at 100 TeV:_ LHAASO reports the detection of 12
Galactic sources Cao  _et al._ (2021a). This search covers a sky region of
$-15^{\circ}<\delta<75^{\circ}$. Because the details of the sensitivities are
not reported, we assume a sensitivity of 0.4 Crab at 100 TeV at best ($4\times
10^{-13}$ erg cm-2 s-1), which is comparable to the minimal flux reported in
Ref. Cao  _et al._ (2021a). We assume a PSF size of $\theta_{\rm
PSF}=0.3^{\circ}$ Aharonian  _et al._ (2021). The sensitivities are degraded
closer to the edge of observable sky region; due to the lack of detailed
information, we assume a declination dependence of the sensitivity as in the
2HWC survey for a $E^{-2.5}$ spectrum Abeysekara  _et al._ (2017a) . This flux
sensitivity corresponds to a distance range of 1.5 kpc for $\Gamma_{\rm
CR}=0.03$ yr-1 and $n_{\rm src}^{\rm gas}=1$ cm-3, as per Eq. (11). This
distance being moderately small indicates that LHAASO and future experiments
have many more sources to discover.

  * •

_Source counts from IceCube:_ The latest results from full-sky point-source
searches are reported in Ref. Aartsen  _et al._ (2020), which uses ten years
of track-like events data and looks for the clustering of neutrino events over
background. Additional searches with seven years of cascade events data are
presented in Ref. Aartsen  _et al._ (2019a), which is more sensitive to the
sources located in the southern sky. Both searches reported the non-detection
of any Galactic sources. These papers report “5$\sigma$ discovery potential”
flux levels as a function of source declination, for two assumed spectrum
slopes, $-2.0$ and $-3.0$. We use the former to be conservative, as we are
fixing the flux near 50 TeV. We consider sources to be detected when the
source flux at 50 TeV exceeds 5$\sigma$ discovery potential of either track or
cascade searches; tracks are more sensitive for a large portion of the sky,
while cascade are slightly better for sources at $\sin\delta\lesssim-0.3$. The
IceCube analysis method integrates the significance over energy, putting more
weight on high energies Aartsen  _et al._ (2019a, 2020). Based on our
estimates, the most important energy range is somewhat below the 50 TeV we
use, but the sensitivity calculations can be improved in future work. The best
detectable flux level reaches $10^{-12}$ erg cm-2 s-1 at 50 TeV. We assume
$\theta_{\rm PSF}=0.3^{\circ}$ for track and $\theta_{\rm PSF}=10^{\circ}$ for
cascade events. This flux sensitivity corresponds to a distance horizon of 0.6
kpc for $\Gamma_{\rm CR}=0.03$ yr-1 and $n_{\rm src}^{\rm gas}=1$ cm-3 (Eq.
12). This distance being very small indicates that IceCube has not yet
significantly probed Milky Way sources.

In addition to source counts, we calculate the total emission from the Milky
Way plane:

  * •

_Total gamma-ray flux from sources:_ Tibet AS$\gamma$ reported the detection
of emission from the Milky Way plane Amenomori  _et al._ (2021c). We use the
measured diffuse intensity, $5\times 10^{-8}$ GeV cm-2 s-1 sr-1 at
$E_{\gamma}=$100 TeV for the region of $25^{\circ}<l<100^{\circ}$ and
$|b|<5^{\circ}$, which includes the truly diffuse ISM emission and
contributions from _unresolved_ sources, but excludes the contribution from
_resolved_ (known) TeV sources. The total emission from known sources is on
the order of 10$\%$ level Amenomori  _et al._ (2021c); Lipari and Vernetto
(2018), so we do not correct for this. We calculate the sum of 100-TeV gamma-
ray fluxes from simulated sources in the region observed by Tibet, and require
that it does not exceed the measured flux.

  * •

_Total neutrino flux from sources:_ IceCube has searched for the diffuse
emission of neutrinos from the Milky Way plane Aartsen  _et al._ (2019a),
which is the sum of the truly diffuse flux plus the contributions from
sources. Unlike Tibet, IceCube has no sources to be removed. IceCube set an
upper limit at $E_{\nu}=50$ TeV ($2\times 10^{-7}$ GeV cm-2 s-1,
$\nu+\bar{\nu}$, all flavor) Aartsen  _et al._ (2017) obtained using templates
for diffuse neutrino emission from Refs. Gaggero  _et al._ (2015, 2017). Note
that the latest cascade search reported a 2$\sigma$ level detection, with the
best-fit flux about half of this upper limit Aartsen  _et al._ (2019a). We
calculate the sum of 50-TeV neutrino fluxes from the Milky Way plane
($|b|<5^{\circ}$) sources, and require that it does not exceed the reported
upper limit.

By using the observations noted above, we connect the CR energy budget near
$E_{p}=$ 1 PeV to gamma-ray data near $E_{\gamma}=$ 100 TeV and neutrino data
near $E_{\nu}=$ 50 TeV, focusing on narrow energy bins that are connected by
typical kinematic relations. To extend our discussions to somewhat lower
energies, we make an additional calculation:

  * •

_Source counts from HAWC at $E_{\gamma}=7$ TeV:_ We calculate the expected
source counts for the latest survey, 3HWC Albert  _et al._ (2020c), which
contains 63 sources (removing two sources that are known to be extragalactic).
We calculate the source flux by extrapolating Eq. (11) with
$d^{2}N_{\gamma}/(dE_{\gamma}dt)\propto(E_{\gamma})^{-2.37}$, motivated by the
CR production rate, although the caveats discussed below Eq. (4) should be
kept in mind. We adopt the differential sensitivity near 7 TeV quoted in Ref.
Albert  _et al._ (2020c) for a spectrum index of $-2.5$ and assume
$\theta_{\rm PSF}=0.3^{\circ}$. The best sensitivity reaches 3 $\times
10^{-13}$ erg cm-2 s-1 (at 7 TeV).

We also calculate expectations for future data.

  * •

_Expected source counts for IceCube-Gen2:_ The sensitivity depends on the
details of the detector and search strategies, which are uncertain, but an
overall factor of five improvement is expected Aartsen  _et al._ (2021). We
estimate the sensitivity of Gen2 by scaling the 5$\sigma$ discovery-level
fluxes from the existing IceCube searches by this factor. We use the same
$\theta_{\rm PSF}$ as for IceCube, although some improvements are expected.
The detection horizon is 1.3 kpc $(\Gamma_{\rm CR}/0.03$ yr${}^{-1})^{-1/2}$
$(n_{\rm src}^{\rm gas}/1$ cm${}^{-3})^{1/2}$.

  * •

_Expected source counts for other neutrino experiments:_ Several neutrino
telescopes, all in the northern hemisphere, are planning to use water instead
of ice. Those include KM3NeT, P-ONE, Baikal-GVD, and TRIDENT Adrián-Martínez
_et al._ (2016); Avrorin  _et al._ (2018); Agostini  _et al._ (2020); Ye  _et
al._ (2022). Those experiments are expected to have better angular resolution
than IceCube, improving the power for finding sources, and also covering a
large portion of sky outside IceCube’s best range. Although the details are
still uncertain, analyses by KM3NeT suggests a source flux sensitivity of
$\simeq$ (0.8–1.6) $\times 10^{-12}$ erg cm-2 s-1, comparable to IceCube, but
over a different range of sky ($-1<\sin\delta<0.8$) Aiello  _et al._ (2019).
We calculate the source-detection expectations for these combined experiments
by assuming a uniform sensitivity over the entire sky (combining with IceCube
or Gen2), using two values of 1.2 (“KM3”) and 0.16 (“KM3$\times 5$”) in units
of $10^{-12}$ erg cm-2 s-1 per flavor at 50 TeV. As a reference, the neutrino
sources detectable with “KM3$\times 5$” sensitivity should have a hadronic
gamma-ray flux of 0.3 Crab at 100 TeV. The detection horizon for the
“KM3$\times 5$” case is still small; 1.5 kpc for $\Gamma_{\rm CR}=0.03$ yr-1
and $n_{\rm src}^{\rm gas}=1$ cm-3. In all cases, we assume $\theta_{\rm
PSF}=0.3^{\circ}$, which is conservative for KM3NeT Adrián-Martínez  _et al._
(2016).

While gamma-ray telescopes typically have better flux sensitivity, neutrino
telescopes have advantages beyond being able to decisively indicated hadronic
sources. First, they cover larger regions of the sky. Second, sources are
typically small compared to the angular resolution, so the flux sensitivity is
less subject to being degraded due to the source extent.

##  V Present Constraints on the Hadronic Source Population

Figure 3: Constraints on PeVatron models from point-source observations, with
the models calibrated to hadronic CR data in the PeV range. Gamma-ray
constraints are based on LHAASO source counts (blue band or dotted line) and
neutrino constraints are based on IceCube non-detections of sources (green
dotted line). Other aspects follow Fig. 2. Left Panel: Cases where 0$\%$ or
100$\%$ of LHAASO sources are hadronic. Right Panel: Same, but for cases with
10$\%$, 50$\%$, or 100$\%$. _Key takeaways:_ _(Left) A wide range of models
remain viable, with present IceCube constraints allowing between 0% and 100%
of the LHAASO sources to be hadronic accelerators, although the latter
requires a large $n_{\rm src}^{\rm gas}\tau_{\rm src}^{\rm esc}$, so that the
source grammage dominates over that of the ISM. (Right) LHAASO is beginning to
probe scenarios where SNRs are in high density environments or have a nearby
molecular cloud, but not yet the ordinary SNR case, despite our being
optimistic on $\tau_{\rm src}^{\rm esc}$._ Figure 4: Same as Fig. 3 (left
panel), adding constraints on sources from the Tibet AS$\gamma$ (a
measurement) and IceCube (a limit) observations of the Milky Way plane. The
alternate axes also apply to the other main text figures. _Key takeaway:_
_Constraints from the Milky Way plane emission are important._ Figure 5:
Constraints from HAWC source measurements (at $E_{\gamma}=7$ TeV). $X_{\rm
ISM}$ is scaled to the value for $E_{p}=70$ TeV. _Key takeaway:_ _Large values
of $n_{\rm src}^{\rm gas}$ and $\tau_{\rm src}^{\rm esc}$ are needed to have
all HAWC sources to be hadronic, suggesting that leptonic sources may be
dominant._

In this section, we develop new constraints on hadronic sources using present
data on CRs, gamma rays, and neutrinos, as discussed above. We stress that the
$n$–$\tau$ plane is constructed from the energy budget of PeV protons, and is
thus consistent with CR observations.

In short, the LHAASO highest-energy sources are consistent with being either
all hadronic or all leptonic. The IceCube non-detections still allow a wide
range of PeVatron models, indicating the need for upgraded neutrino
telescopes. The source grammage in the PeV range is likely subdominant
compared to that in the ISM. Next, we quantify these points.

Figure 3 (left panel) shows important results on gamma-ray point-source
observations. The LHAASO sources may all be hadronic. We show by a star one
theoretical model that is nominally consistent with this, but we note that
these models are optimistic and in conflict with IceCube data (below). It is
also possible that none of the LHAASO sources are hadronic. This is hinted at
by the disagreement of the LHAASO band with the fiducial hadronic model of PeV
protons from SNRs interacting with gas of density 1 cm-3, characteristics of
the ISM. It may thus be that hadronic CR accelerators are so thin (i.e., have
small $n_{\rm src}^{\rm gas}\tau_{\rm src}^{\rm esc}$) that they are not
observable multi-messenger sources.

Figure 3 (left panel) shows important results on neutrino point-source
observations. The IceCube non-detections exclude scenarios with very high
average gas densities (and part of the LHAASO “100$\%$” band on the right),
and that is quite useful. However, many other scenarios remain allowed,
including those where all the LHAASO sources are hadronic. More sensitive
neutrino observations are needed.

Figure 3 (right panel) shows further details about the gamma-ray data. We show
bands assuming that 100%, 50%, and 10% of the LHAASO sources are hadronic
PeVatrons. As the total source count is 12, the case of 10$\%$ indicates that
about one of them is hadronic. If one source is hadronic, parameters within
the “10$\%$” band are required. If zero sources are hadronic, the limit is
invoked; it is higher than the “10$\%$” band because we conservatively allow
statistical fluctuations. The “normal” SNR model is slightly outside the
“10$\%$” band. If the gas density surrounding typical SNRs is slightly larger
than 1 cm-3 due to, e.g., compression by the shock, there is a chance that the
LHAASO sources include one hadronic accelerator, though the uncertainties are
large. Overall, this figure suggests that the highest-energy gamma-ray sources
are likely dominated by leptonic accelerators, although they might contain one
(or even a few) hadronic accelerators.

SNRs have long been the leading candidate for hadronic PeVatrons. However, no
SNRs show evidence of gamma-ray emission beyond 100 TeV, increasing the
impression that they (or at least the majority of them) might not accelerate
PeV hadrons. But our results show that even in the scenario where all ordinary
SNRs are hadronic PeVatrons, LHAASO most likely does not expect to see 100 TeV
gamma rays from them. LHAASO is only starting to probe scenarios where SNRs
have high densities or nearby molecular clouds.

Figure 3 also shows the relative importance of unresolved hadronic sources and
truly diffuse emission. The diagonal dashed line represents $X=X_{\rm ISM}$,
where source column densities for CRs are comparable to what they will
encounter in the ISM before escaping the Galaxy. The grammage accumulated in
the source ($X_{\rm src}$) could be comparable to $X_{\rm ISM}$; in that case,
the LHAASO source counts must be dominated by hadronic sources (except for the
case of very high $n_{\rm src}^{\rm gas}$ and small $\tau_{\rm src}^{\rm
esc}$.) If only a fraction of LHAASO sources are hadronic, $X_{\rm src}$
should be smaller. This means that the total (source + diffuse) Milky Way
luminosity in the PeV range is likely dominated by the diffuse emission from
the ISM. This conclusion is affected by the uncertain source size (discussed
below) and escape time from the Galaxy (a shorter escape time would increase
the relative importance of sources).

Figure 4 shows the effects of taking into account the total gamma-ray and
neutrino fluxes from the Milky Way plane. The gamma-ray constraint from Tibet
AS$\gamma$ is compatible with the constraints by LHAASO source counts. If one
can isolate contributions from the truly diffuse ISM emission (which may be
dominant), the constraints will get even stronger. The constraint from the
IceCube neutrino total flux is comparable to the constraint from IceCube
source counts; both are in general somewhat weaker than gamma-ray constraints.
Figures 4 also displays the gamma-ray luminosity (the corresponding neutrino
luminosity is a factor of 2 smaller than this) and total source number in the
Milky Way as a reference; the same scale on Fig. 4 applies to all other
figures in the main text.

The source-count constraints get worse if the sources is larger or the
observational angular resolution is worse, while the total-emission
constraints do not. These indicate the potential power of the diffuse
measurements to probe the origin of PeVatrons, especially if they are extended
or have low gas densities. If IceCube confirms the emission from the Milky Way
plane, it would define an allowed band in the $n$–$\tau$ plane, similar to
that for the LHAASO sources. For the cases considered here, while source
crowding could be an issue for cases with a large escape time, the background
induced by faint sources in a given direction is never a problem.

Figure 5 extends our gamma-ray point-source considerations to somewhat lower
energies, where there are more detections. Comparing this to Fig. 3 shows that
there is no overlap between “100$\%$ Hadronic” areas for LHAASO and HAWC. If
the escape times for these two particle energies were the same, having all of
the HAWC sources be hadronic would be in conflict with the LHAASO
observations. However, HAWC data are based on gamma-ray measurements at 7 TeV,
probing $E_{p}\sim$ 70 TeV, while LHAASO data are $E_{p}\sim$ 1 PeV. It is
reasonable to expect that the escape time depends on energy, being longer for
lower-energy particles. Still, this analysis shows that the required $n_{\rm
src}^{\rm gas}$ and $\tau_{\rm src}^{\rm esc}$ values required to explain all
the HAWC sources with hadronic accelerators are high, suggesting a significant
contribution from leptonic sources or hadronic non-PeVatrons.

In the Appendix, we explore the following three model variations and their
corresponding $n$–$\tau$ planes.

  1. 1.

It may be that only a fraction of supernova explosions produce PeV hadrons. In
the scenario where $\Gamma_{\rm CR}=0.003$ yr-1 (10% of the SN rate), even the
“SNR (Dense)” point predicts that only one of the LHAASO sources are hadronic.
This scenario is also allowed by existing IceCube data.

This scenario should be distinguished from the case where all SNRs are
PeVatrons, but only a fraction reside in dense environments. Such a scenario
can be studied by Fig. 3. For example, it might be that 90% of SNRs are in a
“normal” environment, while 10% are in a dense environment. These cases
predict 10% and 100% gamma-ray source count (Fig. 3). In total, hadronic
PeVatrons may be $90\%\times 0.1+10\%\times 1\simeq 20\%$ of the highest-
energy gamma-ray sources, i.e., 2 or 3 of them.

  2. 2.

If the emitting regions are larger than 10 pc, more sources escape detection,
reducing the capability of source searches. In scenarios where $R_{\rm
src}=30$ pc, the standard “SNR” points are now further outside the “10$\%$”
band, indicating that they do not appear in the highest-energy gamma-ray
source count even with a factor of $\sim$10 increase in the gas densities. An
important new result is that the possibility that $X_{\rm src}>X_{\rm ISM}$ is
allowed over a wider parameter space than the case of 10-pc sources, which
suggests that sources could make a dominant contribution to the total Galactic
emission.

  3. 3.

Finally, hadronic accelerators might produce CRs over a long time, rather than
impulsively. A scenario that attracts increasing attention is stellar winds in
young star clusters over Myr timescales. We consider the scenario where they
are hadronic PeVatrons. The total number of OB stars in the Milky Way is
estimated to be ${\sim 10^{5}}$ Seo  _et al._ (2018). The number of OB stars
contained in a cluster, $N_{\rm OB}$, is distributed as ${dn/d\log{N_{\rm
OB}}\propto 1/N_{\rm OB}}$. This may suggest that clusters with smaller
$N_{\rm OB}$ contribute more to the total energetics. However, it is
reasonable to think that a larger clusters might have collective effects that
make particle acceleration efficient. Here, we consider clusters of $N_{\rm
OB}\gtrsim 10^{2}$ to be efficient PeVatrons Gupta  _et al._ (2018); the
number of such clusters are about $10^{2}$. The kinetic power is $\sim
10^{38}$ erg s-1 per cluster and $\sim 10^{40}$ erg s-1 in total for the
Galaxy. Comparison of this to the PeV CR energy budget indicates that 1% of
the kinetic energy must be converted to PeV CRs. This requires 10% of the
energy transferred to CRs with a hard spectrum, such as $E^{-2}$.

Our estimates for star clusters are encouraging. They may contribute half of
the LHAASO source count, and future Gen2 observations might find as many as
three of them. Indeed, one of the LHAASO sources is from the region that
contains the Cygnus OB association of massive stars. At the same time, if many
of LHAASO sources are associated with star clusters, we might expect more
associations than one. This might suggest that we are somewhat in the downward
fluctuations. Alternatively, it could be that instead of continuous injection,
impulsive events (e.g., SNR shocks embedded in compact clusters) could be the
PeV hadron sources. Solid conclusions require more work on both theory and
observation.

##  VI Future Prospects for Probing the Hadronic Source Population

Figure 6: Sensitivity of Gen2 to Galactic hadronic PeVatrons. Left: The case
of detections. Right: The case of non-detections. The other features follow
Fig. 3 (left panel). _Key takeaways:_ _Over a wide parameter space, it is
promising that Gen2 will find more than one PeVatron, though it may find zero.
Non-detection by Gen2 would rule our models where 100% of the LHAASO sources
are hadronic and also where $X_{\rm src}>X_{\rm ISM}$, but would still allow a
wide range of hadronic models._

In this section, we outline the power of upcoming observations to find
hadronic PeVatrons. In Sec. VI.1, we calculate the prospects for future
neutrino observations, which will be the most decisive. In Sec. VI.2, we
discuss improved gamma-ray and other observations, which will be available
sooner. In Sec. VI.3, we discuss how our population approach, combined with
improvements in multi-messenger studies of individual sources, will be
critical to solving the long-standing mystery of the Milky Way’s hadronic
PeVatrons.

###  VI.1 Future Neutrino Observations

IceCube will likely be upgraded to Gen2, which will be much more powerful.
Here we estimate its potential to discover Milky Way sources.

Figure 6 (left panel) shows a scenario where Gen2 does detect Milky Way
sources, with the examples of more than one or three sources. As shown, this
would strongly constrain the properties of PeVatrons, limiting the parameter
space to one of the green areas, as labeled. The joint parameter space where
the LHAASO sources can be 100% hadronic is even smaller. The upper right
region of the figure (white space) would be ruled out by LHAASO and IceCube
not detecting more sources.

Figure 6 (right panel) shows a scenario where Gen2 does not detect Milky Way
sources. This would rule out the parameter space where 100% of LHAASO sources
are hadronic and also where $X_{\rm src}$ is larger than $X_{\rm ISM}$. The
latter is particularly important, as it would indicate that the total Galactic
emission should be dominated by the truly diffuse ISM emission (or leptonic
sources). Nevertheless, Gen2 non-detection could still allow a wide range of
hadronic models.

The probability that Gen2 will find a source can be viewed optimistically or
pessimistically. If the source parameters are in the band of “$>$1 Source,”
there is a good chance that Gen2 finally detects the Milky Way’s PeVatrons.
However, even for models that predict one detection, the real source count can
easily fluctuate down to zero. On the other hand, if no sources are observed,
the region above the dotted line in Fig. 6 (right panel) would be robustly
excluded, even in the presence of statistical fluctuations. To put it
differently, the left and right panels of Fig. 6 displays prospects for Gen2
in the presence of favorable and unfavorable statistical fluctuations. Note
that sources Gen2 would discover are most likely those already found by
LHAASO, but not necessarily. Gen2, especially with a cascade analysis, can in
principle observe a much wider sky region and be sensitive to even extended
sources.

Even if Gen2 detects sources, it may not be many, and further improvements in
sensitivity will be needed to decisively probe Milky Way hadronic
accelerators. This conclusion is less optimistic than early work that
suggested that IceCube could detect several sources Kistler and Beacom (2006);
Beacom and Kistler (2007); Halzen  _et al._ (2008); Gonzalez-Garcia  _et al._
(2009). Those early studies assumed that most gamma-ray sources were hadronic
and the calculations were not constrained by detailed measurements of source
properties and calibrated to the CR energy budget, as here. To improve
sensitivity to Milky Way sources, ocean-based detectors may be especially
important because of their better angular resolution Adrián-Martínez  _et al._
(2016); Avrorin  _et al._ (2018); Agostini  _et al._ (2020); Ye  _et al._
(2022). Also, those detectors are planned to be in the Northern hemisphere,
and hence would have better sensitivity to the inner Milky Way.

Figure 7 shows how such future neutrino telescopes could probe PeVatron
models. We show the cases of non-detection with two sensitivities (“KM3” and
“KM3$\times$5”) as defined above. Non-detection by “KM3$\times$5” level
experiments would be very constraining. The limit would surpass even the “MC
near SNR” point by a factor of 4, ruling out scenarios where more than quarter
of SNRs coincide with MCs that produce gamma rays. Nevertheless, there would
still remain a wide range of models allowed, where hadronic PeVatrons are
thin, indicating the need for still-better sensitivity.

Figure 7: Projected constraints from future neutrino experiments for the case
of non-detection. Labels are same as the right panel of Fig. 6 (except for
Gen2) and omitted. Some arrows are removed to improve visibility. The neutrino
sensitivities are assumed to be uniform across the sky with values as expected
for KM3NeT (marked with “KM3”) and five times better than that (“KM3$\times
5$”). _Key takeaway:_ _To definitely probe scenario where hadronic PeVatrons
are thin, sensitivity better than that of even Gen2 and KM3NeT is needed._

One idea to attain better sensitivities is further efforts to improve the
angular reconstruction of the neutrino cascade events. In the current neutrino
source searches, the best sensitivities are obtained by the analysis of track-
like events due to muons. As pointed out by Ref. Beacom and Candia (2004),
electron neutrinos have much smaller backgrounds than muon neutrinos, making
this channel interesting if the angular resolution can be made sufficiently
good. We will report on this in a forthcoming paper.

Our analysis does not constrain CR-hidden sources, from which CR do not
escape. Though IceCube has not yet identified any such sources, Gen2 or other
future experiments might.

###  VI.2 Future Multi-Messenger Observations

New multi-messenger observations — covering the full electromagnetic spectrum
(not only gamma rays) plus cosmic rays — will also be important.

The most powerful input we have so far for constraining hadronic PeVatrons is
the count of LHAASO gamma-ray sources near the PeV range. In the next few
years, even better results are expected as the LHAASO construction is
completed, the observation time is increased, and more theoretical modeling is
done. The LHAASO source count provides firm upper limits on the properties of
hadronic PeVatrons, generally stronger than even those based on the non-
detection of IceCube neutrino sources. In the forthcoming years, as Tibet
AS$\gamma$, HAWC, and especially LHAASO increase the source counts in the 100
TeV range, this will lead to better constraints in the $n$–$\tau$ plane. Their
diffuse measurements will also help, primarily by constraining the properties
of the high-energy cosmic-ray spectrum.

In the further future, the Cherenkov Telescope Array (CTA) Acharya  _et al._
(2019), which will have outstanding flux sensitivity and angular resolution,
will discover many sources and will be able to conduct detailed morphological
studies, probing the spatial correlations between gamma-ray emission and gas
density. A possible future observatory is the Southern Wide-field Gamma-ray
Observatory (SWGO) Albert  _et al._ (2019). Like Tibet AS$\gamma$, HAWC, and
LHAASO, this would be a large ground array that simultaneously views a wide
field on the sky. Such detectors are more sensitive to extended sources than
pointed observatories like CTA. SWGO would be located in the Southern
hemisphere, providing a good view of the inner Galaxy, where a high density of
bright, interesting sources is found. Another experiment, the Andes Large-area
PArticle detector for Cosmic-ray physics and Astronomy (ALPACA), is also
planned to observe highest-energy gamma-ray sources in the Southern hemisphere
Kato  _et al._ (2021).

In principle, gamma-ray observations could separate hadronic and leptonic
sources without using neutrinos, though decisive answers have been challenging
to obtain. A lot of attention has been placed on specifically identifying
hadronic sources. This might be done through some combination of the spectrum
shape, spatial correlations between emission and gas density, and multi-
wavelength observations and modeling. Another approach would be specifically
identifying leptonic sources, using tests of the spectrum shape and multi-
wavelength studies. Crucially, the electrons that produce gamma rays through
inverse-Compton emission must also produce synchrotron x-rays. In the limit
that leptonic sources are dominant, this may be a more fruitful approach. If
more sources are identified as leptonic, the allowed parameter space in the
$n_{\rm src}^{\rm gas}$–$\tau_{\rm src}^{\rm esc}$ would get narrower and
lower, providing stronger constraints on the PeVatron models.

One of the key inputs in our population model is the escape time $\tau_{\rm
src}^{\rm esc}$. The recent discovery of “TeV halos” around pulsars, where the
escape time is quite long, has demonstrated the power of gamma rays to probe
the escape of very high energy particles from accelerators. Near-future
observations will certainly find more halos around CR accelerators, which
would be key input to theoretical efforts in understanding the confinement by
both leptonic and hadronic sources.

Multi-wavelength observations will help identify sources and characterize
their properties. In particular, radio and x-ray observations of gamma-ray
objects are important. First, they allow us to isolate components (e.g., SNRs,
pulsar wind nebulae, star clusters) due to having good angular resolution.
Second, these data probe shock structures and magnetic fields in the region,
which is a key input to understand particle acceleration. Third, they are
needed to constrain the leptonic component of the gamma-ray emission. In-depth
multi-wavelength observations toward sources will also be needed to
characterize the environment, constraining the values of $n_{\rm src}^{\rm
gas}$.

New CR data will help reduce uncertainties in the CR production rate. First,
an improvement in the measurement of the grammage, especially at higher
energies, would be crucial in determining the CR energy budget. Second, a
better understanding of the CR knee is needed. The existing measurements of CR
protons show discrepancies among observations; the amount is by a factor of
about 3 at $E_{p}\simeq 3$ PeV, which becomes larger for higher energies. This
prohibits us from pointing the exact energy of the proton knee and other
spectrum features. Although the all-particle CR spectrum shows the knee at
$\simeq 3$ PeV, the proton spectrum might have a spectrum break at $\simeq 1$
PeV.

Throughout this work, we use the proton energy budget near 1 PeV and also
gamma-ray data at $E_{\gamma}=100$ TeV. In fact, the LHAASO sources are
observed much above this energy, reaching more than 1 PeV in some cases. This
is curious, because hadronic sources that accelerate protons only to $\sim$1
PeV cannot generate bright emission at such energies, because typically
$E_{\gamma}/E_{p}\ll 1$. Such high-energy photons are also difficult to
produce by leptonic sources; in the Klein-Nishina regime of the inverse-
Compton scattering, the fractional energy from CRs to gamma rays are large,
$E_{\gamma}/E_{e}\sim 1$, but the cross section is suppressed. If we have more
detections of $\sim 1$ PeV photons, it might suggest that the knee is at
$E_{p}>1$ PeV or that there are super-knee sources. As the exact value of
$E_{\gamma}/E_{p}$ can significantly depend on the spectrum properties, better
theoretical and observational characterization of the knee spectrum will be
key.

###  VI.3 Importance of Our Population Approach

Our new theoretical framework will remain valuable even as studies of
individual sources advance. For example, dedicated IceCube studies of LHAASO
sources may lead to stronger constraints on hadronic emission, but those are
not connected to the Milky Way CR energy budget. Even once individual hadronic
sources are identified, our population approach will remain valuable. First,
it might be that one or a few hadronic sources are observed in both gamma rays
and neutrinos, but that this does not necessarily mean that this source class
contributes significantly to the Milky Way’s cosmic-ray budget. In our
approach, this can be tested by examining if the source parameters, $n_{\rm
src}^{\rm gas}$ and $\tau_{\rm src}^{\rm esc}$, are consistent with the
various constraints from gamma rays and neutrinos. Contradictions would
suggest that either the sources are not representative PeVatrons or that the
standard models of CR production and propagation as adopted above need to be
drastically modified.

Our approach can be improved in a number of ways. Besides the total source
count, other observables can be predicted by our population models. First, the
flux distributions, $dN/dF$, differ for different combinations of $n_{\rm
src}^{\rm gas}$ and $\tau_{\rm src}^{\rm esc}$, even if they predict the same
source count. Our calculations here can be regarded as only using the integral
of $dN/dF$ above $F_{\rm lim}$. In principle, the use of full distribution can
further constrain models. For example, our approach may be allowing too much
of the bottom-right region, where sources are luminous and rare. Second, the
latitude distribution could help to constrain models where gamma-ray sources
have low gas densities. If the sources are low-luminosity and nearby, then the
latitude distribution should have a larger scatter compared to the case of
high-luminosity and distant sources. This might further constrain models in
the upper left part the plane. Quantitatively, since the detection horizon for
gamma-ray sources is $\simeq 3\sqrt{n_{\rm src}^{\rm gas}/\rm cm^{-3}}$ kpc
and the CR sources are distributed within a height of $\simeq 30$ pc
(conservatively small), the scatter in the latitude could be $|b|\simeq
0.6^{\circ}\sqrt{n_{\rm src}^{\rm gas}/\rm cm^{-3}}$.

One difficulty in our approach is that we have to assume a CR source rate to
derive constraints in the $n$–$\tau$ plane, because we consider the _number_
of sources. Such an assumption would be eliminated if we instead consider the
total _flux_ from these sources, which we discussed only briefly. Such an
approach may give new and comparable constraints on the PeVatron models.
Alternatively, we could construct a plane of $n_{\rm src}^{\rm gas}$ and
$\Gamma_{\rm CR}$. Another difficulty in our approach is that we assume that
sources are the same as each other, which is a strong assumption, although we
briefly discussed how to consider source-to-source variations in an
approximate manner. At the present level of precision, this is reasonable, but
ultimately it will be necessary to take variations into account in a more
sophisticated manner.

We encourage the development of public codes to model point sources and
collections thereof. At the moment, there seem to be none. This is in contrast
to extensive work on modeling diffuse gamma-ray and neutrino emission in the
Milky Way, such as with GALPROP Moskalenko and Strong (1998); Strong and
Moskalenko (1998); Strong  _et al._ (2000); Moskalenko  _et al._ (2002);
Strong  _et al._ (2004); Ptuskin  _et al._ (2006); Porter  _et al._ (2017,
2019); Boschini  _et al._ (2020); Porter  _et al._ (2021), DRAGON Gaggero  _et
al._ (2014); Evoli  _et al._ (2017, 2018), USINE Maurin  _et al._ (2001);
Putze  _et al._ (2009, 2010); Maurin (2020), PICARD Kissmann (2014); Kissmann
_et al._ (2015), and CRIPTIC Krumholz  _et al._ (2022)), some of which are
public.

##  VII Conclusions

The Milky Way contains powerful but unidentified accelerators of PeV hadronic
CRs. And gamma-ray point sources have been observed at energies into the PeV
range. However, no neutrino sources have been detected. As illustrated in Fig.
1, these observations leave two extreme possibilities for the mysterious
hadronic PeVatrons: (1) Hadronic CR sources do produce both gamma rays and
neutrinos, making them exciting multi-messenger targets, but greater
sensitivity is needed to detect the neutrinos, versus (2) Hadronic CR sources
are so thin that they produce no detectable gamma rays or neutrinos, so that
the observed gamma-ray sources are only leptonic accelerators.

The main aim of this paper is to understand where the hadronic PeVatrons lie
between those extremes. We introduce a new population-based approach that
characterizes sources in the plane of source gas density and CR escape time
(the $n$–$\tau$ plane), explicitly calibrating point-source models to CR
observations. We quantify the counts of identifiable hadronic sources and the
contributions of hadronic sources to the total emission from the Milky Way
plane. We calculate the ranges of allowed models for hadronic accelerators.

We compare predictions for gamma-ray sources to data from LHAASO, which has
identified PeV-range sources. In the optimistic interpretation, these sources
could all be hadronic CR accelerators and hence multi-messenger sources. This
requires models to be within the “100$\%$” band in the $n$–$\tau$ plane shown
in Fig. 3. In this case, PeV hadrons accumulate a comparable grammage in the
sources and in the ISM. In a pessimistic interpretation, the LHAASO sources
could be all leptonic. This is allowed over a wide parameter space, quantified
in the “LHAASO (0$\%$ Hadronic)” limit in Fig. 3 (left). In this case, the
grammage accumulated in the ISM would dominate. Possibly the most realistic
scenario is that only one or a few of the LHAASO sources is hadronic, as in
Fig. 3 (right).

Neutrino observations are the key to decisively resolving the nature of the
hadronic PeVatrons. We show that IceCube non-detections rule out models with
high gas densities, but still leave a wide parameter space open, as per the
“IceCube (Source)” limit in Fig. 3. Improvements in the neutrino sensitivities
are needed. Gen2 is promising for finally identifying the hadronic PeVatrons,
as discovery can be expected for a wide range of the parameter space (Fig. 6,
left). Non-detections by Gen2 and KM3NeT would constrain significant and
important parameter space, but even then, a substantial discovery space for
multi-messenger source would remain. To quantify that, the detection range of
Gen2 for is only $\sim 1.3\sqrt{n_{\rm src}^{\rm gas}/\rm cm^{-3}}$ kpc (for
$\Gamma_{\rm CR}=0.03$ yr-1). Even larger detectors may be needed to fully
probe the origins of hadronic CRs in the Milky Way.

Identifying the origins of the cosmic rays is a century-old problem. While it
is possible that new observations may soon lead to breakthroughs, it is also
possible that this problem will remain challenging for decades more. A key
focus of work on Milky Way hadronic CR accelerators has been on trying to
determine if individual gamma-ray sources are hadronic or leptonic, with clear
progress but not meeting the goal of definitive answers for most sources. And
even when this goal is achieved, it remains a separate question to decide if
these sources are producing enough CRs to account for the Milky Way fluxes.
This paper is a first step in starting a new, population-based approach to
revealing the hadronic PeVatrons. Our hope is that it will be improved through
new observational constraints, new ideas for theoretical constraints, and the
development of public codes similar in sophistication to those used for
modeling the Galactic diffuse emission.

In the near term, results from IceCube on observations of the total neutrino
emission from the Milky Way plane will be quite important, especially if the
hint of a signal strengthens in significance. As shown in Fig. 4, an IceCube
signal near their present “Total” limit would suggest that source emission is
more important than diffuse emission and that the LHAASO gamma-ray sources
should be dominantly hadronic. Ironically, while that would confirm the
existence of hadronic PeVatrons, it would still not reveal where they are.
Still, it would indicate that the optimistic case in Fig. 1 is likely true,
indicating a bright future for multi-messenger astronomy as cosmic-ray, gamma-
ray, and neutrino observatories gain in sensitivity. In turn, that would imply
excellent prospects not only for resolving long-standing questions in
astrophysics, but also for developing new tests of physics beyond the standard
model using these sources.

## Acknowledgments

We are grateful for helpful discussions with Katsuaki Asano, Aya Bamba,
Stefano Gabici, Francis Halzen, William Luszczak, Yutaka Ohira, Hidetoshi
Sano, Hiromasa Suzuki, and especially Tim Linden and Kohta Murase. T.S. was
primarily supported by an Overseas Research Fellowship from the Japan Society
for the Promotion of Science (JSPS). T.S. was partially supported by and
J.F.B. was fully supported by National Science Foundation grant No.
PHY-2012955. This research made use of matplotlib Hunter (2007) and numpy van
der Walt  _et al._ (2011). Computing resources were provided by the Ohio
Supercomputer Center Ohio Supercomputer Center (1987).

##  Appendix A Population Model for Alternative Scenarios

In the main text, we focused on the case of impulsive CR injections at a rate
of $\Gamma_{\rm CR}=0.03~{}$yr-1 and assumed the size of the gamma-ray and
neutrino emitting regions to $R_{\rm src}=10$ pc. Here we present figures for
three alternative scenarios:

  1. 1.

Figures 8 and 11 show scenarios where the CR source rate is 0.003 yr-1, 10% of
the Galactic SN rate.

  2. 2.

Figures 9 and 12 show scenarios where $R_{\rm src}=30$ pc.

  3. 3.

Figures 10 and 13 show scenarios where the injections of CRs are continuous.
We use $\mathcal{N}_{\rm CR}=10^{2}$, as discussed in the main text, and we
assume $R_{\rm src}=30$ pc. Also, we show a star to represent the case of
“Star Cluster” sources, as discussed in the main text.

Figures 8–10 are the same as the right panel of Fig. 3, and Figs. 11–13 are
the same as the left panel of Fig. 6. In these appendix figures, arrows are
omitted for clarity.

Figure 8: Same as Fig. 3 (right panel), but for $\Gamma_{\rm CR}$ = 0.003
yr-1. Figure 9: Same as Fig. 3 (right panel), but for $R_{\rm src}=30$ pc.
Figure 10: Same as Fig. 3 (right panel), but for the case of continuous CR
injection with $\mathcal{N}_{\rm CR}=10^{2}$ and $R_{\rm src}=30$ pc. Figure
11: Same as Fig. 6 (right panel), but for $\Gamma_{\rm CR}$ = 0.003 yr-1.
Figure 12: Same Fig. 6 (right panel), but for $R_{\rm src}=30$ pc. Figure 13:
Same as Fig. 6 (right panel), but for the case of continuous CR injection with
$\mathcal{N}_{\rm CR}=10^{2}$ and $R_{\rm src}=30$ pc.

## References

  * Ginzburg and Syrovatskii (1964) V. L. Ginzburg and S. I. Syrovatskii,  _The Origin of Cosmic Rays_ (1964).
  * Berezinskii  _et al._ (1990) V. S. Berezinskii, S. V. Bulanov, V. A. Dogiel, and V. S. Ptuskin,  _Astrophysics of cosmic rays_ (1990).
  * Gaisser (1990) T. K. Gaisser,  _Cosmic rays and particle physics._ (1990).
  * Gaisser  _et al._ (2016) T. K. Gaisser, R. Engel, and E. Resconi,  _Cosmic Rays and Particle Physics_ (2016).
  * Lagage and Cesarsky (1983) P. O. Lagage and C. J. Cesarsky, A&A 125, 249 (1983).
  * Bell (2004) A. R. Bell, MNRAS 353, 550 (2004).
  * Antoni  _et al._ (2005) T. Antoni _et al._ , Astroparticle Physics 24, 1 (2005), arXiv:astro-ph/0505413 [astro-ph] .
  * Bell  _et al._ (2013) A. R. Bell, K. M. Schure, B. Reville, and G. Giacinti, MNRAS 431, 415 (2013), arXiv:1301.7264 [astro-ph.HE] .
  * Blasi (2013) P. Blasi, A&A Rev. 21, 70 (2013), arXiv:1311.7346 [astro-ph.HE] .
  * Anchordoqui  _et al._ (2014) L. A. Anchordoqui _et al._ , Journal of High Energy Astrophysics 1, 1 (2014), arXiv:1312.6587 [astro-ph.HE] .
  * Cristofari (2021) P. Cristofari, Universe 7, 324 (2021), arXiv:2110.07956 [astro-ph.HE] .
  * Stecker (1971) F. W. Stecker,  _Cosmic gamma rays_ , Vol. 249 (1971).
  * Aharonian (2004) F. A. Aharonian,  _Very high energy cosmic gamma radiation : a crucial window on the extreme Universe_ (2004).
  * Dermer and Menon (2009) C. D. Dermer and G. Menon,  _High Energy Radiation from Black Holes: Gamma Rays, Cosmic Rays, and Neutrinos_ (2009).
  * Abdo  _et al._ (2007) A. A. Abdo _et al._ , ApJ 664, L91 (2007), arXiv:0705.0707 [astro-ph] .
  * Abeysekara  _et al._ (2017a) A. U. Abeysekara _et al._ , ApJ 843, 40 (2017a), arXiv:1702.02992 [astro-ph.HE] .
  * Abdalla  _et al._ (2018a) H. Abdalla _et al._ (H. E. S. S. Collaboration), A&A 612, A1 (2018a), arXiv:1804.02432 [astro-ph.HE] .
  * Abdollahi  _et al._ (2020) S. Abdollahi _et al._ , ApJS 247, 33 (2020), arXiv:1902.10045 [astro-ph.HE] .
  * Cao  _et al._ (2021a) Z. Cao _et al._ , Nature 594, 33 (2021a).
  * Atoyan and Aharonian (1996) A. M. Atoyan and F. A. Aharonian, MNRAS 278, 525 (1996).
  * Ahlers and Halzen (2018) M. Ahlers and F. Halzen, Progress in Particle and Nuclear Physics 102, 73 (2018), arXiv:1805.11112 [astro-ph.HE] .
  * Aartsen  _et al._ (2019a) M. G. Aartsen _et al._ (IceCube Collaboration), ApJ 886, 12 (2019a), arXiv:1907.06714 [astro-ph.HE] .
  * Aartsen  _et al._ (2020) M. G. Aartsen _et al._ , Phys. Rev. Lett. 124, 051103 (2020), arXiv:1910.08488 [astro-ph.HE] .
  * Albert  _et al._ (2020a) A. Albert _et al._ (ANTARES Collaboration and Icecube Collaboration), ApJ 892, 92 (2020a), arXiv:2001.04412 [astro-ph.HE] .
  * Kheirandish (2020) A. Kheirandish, Ap&SS 365, 108 (2020), arXiv:2006.16087 [astro-ph.HE] .
  * Kistler and Beacom (2006) M. D. Kistler and J. F. Beacom, Phys. Rev. D 74, 063007 (2006), arXiv:astro-ph/0607082 [astro-ph] .
  * Beacom and Kistler (2007) J. F. Beacom and M. D. Kistler, Phys. Rev. D 75, 083001 (2007), arXiv:astro-ph/0701751 [astro-ph] .
  * Halzen  _et al._ (2008) F. Halzen, A. Kappes, and A. Ó Murchadha, Phys. Rev. D 78, 063004 (2008), arXiv:0803.0314 [astro-ph] .
  * Gonzalez-Garcia  _et al._ (2009) M. C. Gonzalez-Garcia, F. Halzen, and S. Mohapatra, Astroparticle Physics 31, 437 (2009), arXiv:0902.1176 [astro-ph.HE] .
  * Aartsen  _et al._ (2021) M. G. Aartsen _et al._ , Journal of Physics G Nuclear Physics 48, 060501 (2021), arXiv:2008.04323 [astro-ph.HE] .
  * Adrián-Martínez  _et al._ (2016) S. Adrián-Martínez _et al._ , Journal of Physics G Nuclear Physics 43, 084001 (2016), arXiv:1601.07459 [astro-ph.IM] .
  * Avrorin  _et al._ (2018) A. D. Avrorin _et al._ (Baikal-GVD Collaboration), arXiv e-prints , arXiv:1808.10353 (2018), arXiv:1808.10353 [astro-ph.IM] .
  * Agostini  _et al._ (2020) M. Agostini _et al._ , Nature Astronomy 4, 913 (2020), arXiv:2005.09493 [astro-ph.HE] .
  * Ye  _et al._ (2022) Z. P. Ye _et al._ , arXiv e-prints , arXiv:2207.04519 (2022), arXiv:2207.04519 [astro-ph.HE] .
  * Panov  _et al._ (2009) A. D. Panov _et al._ , Bulletin of the Russian Academy of Sciences, Physics 73, 564 (2009), arXiv:1101.3246 [astro-ph.HE] .
  * Adriani  _et al._ (2011) O. Adriani _et al._ , Science 332, 69 (2011), arXiv:1103.4055 [astro-ph.HE] .
  * Aguilar  _et al._ (2015) M. Aguilar _et al._ (AMS Collaboration), Phys. Rev. Lett. 114, 171103 (2015).
  * Yoon  _et al._ (2017) Y. S. Yoon _et al._ , ApJ 839, 5 (2017), arXiv:1704.02512 [astro-ph.HE] .
  * Atkin  _et al._ (2018) E. Atkin _et al._ , Soviet Journal of Experimental and Theoretical Physics Letters 108, 5 (2018), arXiv:1805.07119 [astro-ph.HE] .
  * An  _et al._ (2019) Q. An _et al._ , Science Advances 5, eaax3793 (2019), arXiv:1909.12860 [astro-ph.HE] .
  * Adriani  _et al._ (2019) O. Adriani _et al._ (Calet Collaboration), Phys. Rev. Lett. 122, 181102 (2019), arXiv:1905.04229 [astro-ph.HE] .
  * Höorandel (2006) J. R. Höorandel, in  _Journal of Physics Conference Series_, Journal of Physics Conference Series, Vol. 47 (2006) pp. 41–50, arXiv:astro-ph/0508014 [astro-ph] .
  * Albert  _et al._ (2022) A. Albert _et al._ , Phys. Rev. D 105, 063021 (2022), arXiv:2204.06662 [astro-ph.HE] .
  * Amenomori  _et al._ (2008) M. Amenomori _et al._ (Tibet AS$\gamma$ Collaboration), ApJ 678, 1165 (2008), arXiv:0801.1803 [hep-ex] .
  * Apel  _et al._ (2013) W. D. Apel _et al._ , Astroparticle Physics 47, 54 (2013).
  * Aartsen  _et al._ (2019b) M. G. Aartsen _et al._ (IceCube Collaboration), Phys. Rev. D 100, 082002 (2019b), arXiv:1906.04317 [astro-ph.HE] .
  * Lipari and Vernetto (2020) P. Lipari and S. Vernetto, Astroparticle Physics 120, 102441 (2020), arXiv:1911.01311 [astro-ph.HE] .
  * Ptuskin  _et al._ (2010) V. Ptuskin, V. Zirakashvili, and E.-S. Seo, ApJ 718, 31 (2010), arXiv:1006.0034 [astro-ph.CO] .
  * Ahnen  _et al._ (2017) M. L. Ahnen _et al._ , MNRAS 472, 2956 (2017), arXiv:1707.01583 [astro-ph.HE] .
  * Abeysekara  _et al._ (2020a) A. U. Abeysekara _et al._ , ApJ 894, 51 (2020a), arXiv:2003.13615 [astro-ph.HE] .
  * Cristofari  _et al._ (2013) P. Cristofari, S. Gabici, S. Casanova, R. Terrier, and E. Parizot, MNRAS 434, 2748 (2013), arXiv:1302.2150 [astro-ph.HE] .
  * Cristofari  _et al._ (2018) P. Cristofari, S. Gabici, R. Terrier, and T. B. Humensky, MNRAS 479, 3415 (2018), arXiv:1803.09728 [astro-ph.HE] .
  * Bykov (2014) A. M. Bykov, A&A Rev. 22, 77 (2014), arXiv:1511.04608 [astro-ph.HE] .
  * Fujita  _et al._ (2017) Y. Fujita, K. Murase, and S. S. Kimura, J. Cosmology Astropart. Phys 2017, 037 (2017), arXiv:1604.00003 [astro-ph.HE] .
  * Ioka  _et al._ (2017) K. Ioka, T. Matsumoto, Y. Teraki, K. Kashiyama, and K. Murase, MNRAS 470, 3332 (2017), arXiv:1612.03913 [astro-ph.HE] .
  * Guépin  _et al._ (2018) C. Guépin, L. Rinchiuso, K. Kotera, E. Moulin, T. Pierog, and J. Silk, J. Cosmology Astropart. Phys 2018, 042 (2018), arXiv:1806.03307 [astro-ph.HE] .
  * Ohira  _et al._ (2018) Y. Ohira, S. Kisaka, and R. Yamazaki, MNRAS 478, 926 (2018), arXiv:1702.05866 [astro-ph.HE] .
  * Kimura  _et al._ (2018) S. S. Kimura, K. Murase, and P. Mészáros, ApJ 866, 51 (2018), arXiv:1807.03290 [astro-ph.HE] .
  * Sudoh  _et al._ (2020) T. Sudoh, Y. Inoue, and D. Khangulyan, ApJ 889, 146 (2020), arXiv:1911.00013 [astro-ph.HE] .
  * Cooper  _et al._ (2020) A. J. Cooper, D. Gaggero, S. Markoff, and S. Zhang, MNRAS 493, 3212 (2020), arXiv:2002.01477 [astro-ph.HE] .
  * Kimura  _et al._ (2020) S. S. Kimura, K. Murase, and P. Mészáros, ApJ 904, 188 (2020), arXiv:2008.04515 [astro-ph.HE] .
  * Guépin  _et al._ (2020) C. Guépin, B. Cerutti, and K. Kotera, A&A 635, A138 (2020), arXiv:1910.11387 [astro-ph.HE] .
  * Morlino  _et al._ (2021) G. Morlino, P. Blasi, E. Peretti, and P. Cristofari, MNRAS 504, 6096 (2021), arXiv:2102.09217 [astro-ph.HE] .
  * Kimura  _et al._ (2021) S. S. Kimura, T. Sudoh, K. Kashiyama, and N. Kawanaka, ApJ 915, 31 (2021), arXiv:2103.15029 [astro-ph.HE] .
  * Engelmann  _et al._ (1990) J. J. Engelmann _et al._ , A&A 233, 96 (1990).
  * Blum  _et al._ (2013) K. Blum, B. Katz, and E. Waxman, Phys. Rev. Lett. 111, 211101 (2013), arXiv:1305.1324 [astro-ph.HE] .
  * Aguilar  _et al._ (2016) M. Aguilar _et al._ (AMS Collaboration), Phys. Rev. Lett. 117, 231102 (2016).
  * Blasi  _et al._ (2012) P. Blasi, E. Amato, and P. D. Serpico, Phys. Rev. Lett. 109, 061101 (2012), arXiv:1207.3706 [astro-ph.HE] .
  * Aloisio and Blasi (2013) R. Aloisio and P. Blasi, J. Cosmology Astropart. Phys 2013, 001 (2013), arXiv:1306.2018 [astro-ph.HE] .
  * Moskalenko and Strong (1998) I. V. Moskalenko and A. W. Strong, ApJ 493, 694 (1998), arXiv:astro-ph/9710124 [astro-ph] .
  * Strong  _et al._ (2007) A. W. Strong, I. V. Moskalenko, and V. S. Ptuskin, Annual Review of Nuclear and Particle Science 57, 285 (2007), arXiv:astro-ph/0701517 [astro-ph] .
  * Cowsik and Burch (2010) R. Cowsik and B. Burch, Phys. Rev. D 82, 023009 (2010).
  * Lipari (2017) P. Lipari, Phys. Rev. D 95, 063009 (2017), arXiv:1608.02018 [astro-ph.HE] .
  * Lipari (2019) P. Lipari, Phys. Rev. D 99, 043005 (2019), arXiv:1810.03195 [astro-ph.HE] .
  * Gabici  _et al._ (2019) S. Gabici, C. Evoli, D. Gaggero, P. Lipari, P. Mertsch, E. Orlando, A. Strong, and A. Vittino, International Journal of Modern Physics D 28, 1930022-339 (2019), arXiv:1903.11584 [astro-ph.HE] .
  * Murase and Fukugita (2019) K. Murase and M. Fukugita, Phys. Rev. D 99, 063012 (2019), arXiv:1806.04194 [astro-ph.HE] .
  * Kelner  _et al._ (2006) S. R. Kelner, F. A. Aharonian, and V. V. Bugayov, Phys. Rev. D 74, 034018 (2006), arXiv:astro-ph/0606058 [astro-ph] .
  * Kafexhiu  _et al._ (2014) E. Kafexhiu, F. Aharonian, A. M. Taylor, and G. S. Vila, Phys. Rev. D 90, 123014 (2014), arXiv:1406.7369 [astro-ph.HE] .
  * Celli  _et al._ (2020) S. Celli, F. Aharonian, and S. Gabici, ApJ 903, 61 (2020), arXiv:2009.05999 [astro-ph.HE] .
  * Amenomori  _et al._ (2019) M. Amenomori _et al._ (Tibet AS$\gamma$ Collaboration), Phys. Rev. Lett. 123, 051101 (2019), arXiv:1906.05521 [astro-ph.HE] .
  * Abeysekara  _et al._ (2020b) A. U. Abeysekara _et al._ (HAWC Collaboration), Phys. Rev. Lett. 124, 021102 (2020b).
  * Albert  _et al._ (2020b) A. Albert _et al._ (HAWC Collaboration), ApJ 896, L29 (2020b), arXiv:2005.13699 [astro-ph.HE] .
  * Albert  _et al._ (2021a) A. Albert _et al._ , ApJ 907, L30 (2021a), arXiv:2012.15275 [astro-ph.HE] .
  * Amenomori  _et al._ (2021a) M. Amenomori _et al._ (Tibet AS$\gamma$ Collaboration), Nature Astronomy 5, 460 (2021a), arXiv:2109.02898 [astro-ph.HE] .
  * Amenomori  _et al._ (2021b) M. Amenomori _et al._ (Tibet AS$\gamma$ Collaboration), Phys. Rev. Lett. 127, 031102 (2021b), arXiv:2107.01064 [astro-ph.HE] .
  * Abeysekara  _et al._ (2021) A. U. Abeysekara _et al._ , Nature Astronomy 5, 465 (2021), arXiv:2103.06820 [astro-ph.HE] .
  * Amenomori  _et al._ (2022) M. Amenomori _et al._ , ApJ 932, 120 (2022).
  * Cao  _et al._ (2021b) Z. Cao _et al._ , ApJ 917, L4 (2021b), arXiv:2107.02020 [astro-ph.HE] .
  * Abdalla  _et al._ (2021) H. Abdalla _et al._ , A&A 653, A152 (2021), arXiv:2106.06405 [astro-ph.HE] .
  * Xin  _et al._ (2019) Y. Xin, H. Zeng, S. Liu, Y. Fan, and D. Wei, ApJ 885, 162 (2019), arXiv:1907.04972 [astro-ph.HE] .
  * Fang  _et al._ (2022) K. Fang, M. Kerr, R. Blandford, H. Fleischhack, and E. Charles, Phys. Rev. Lett. 129, 071101 (2022), arXiv:2208.05457 [astro-ph.HE] .
  * Tanaka and Takahara (2011) S. J. Tanaka and F. Takahara, ApJ 741, 40 (2011), arXiv:1108.1690 [astro-ph.HE] .
  * Torres  _et al._ (2014) D. F. Torres, A. Cillis, J. Martín, and E. de Oña Wilhelmi, Journal of High Energy Astrophysics 1, 31 (2014), arXiv:1402.5485 [astro-ph.HE] .
  * Abeysekara  _et al._ (2017b) A. U. Abeysekara _et al._ , Science 358, 911 (2017b), arXiv:1711.06223 [astro-ph.HE] .
  * Linden  _et al._ (2017) T. Linden, K. Auchettl, J. Bramante, I. Cholis, K. Fang, D. Hooper, T. Karwal, and S. W. Li, Phys. Rev. D 96, 103016 (2017), arXiv:1703.09704 [astro-ph.HE] .
  * Abdalla  _et al._ (2018b) H. Abdalla _et al._ (H. E. S. S. Collaboration), A&A 612, A2 (2018b), arXiv:1702.08280 [astro-ph.HE] .
  * Sudoh  _et al._ (2019) T. Sudoh, T. Linden, and J. F. Beacom, Phys. Rev. D 100, 043016 (2019), arXiv:1902.08203 [astro-ph.HE] .
  * Di Mauro  _et al._ (2021) M. Di Mauro, S. Manconi, M. Negro, and F. Donato, Phys. Rev. D 104, 103002 (2021), arXiv:2012.05932 [astro-ph.HE] .
  * Albert  _et al._ (2021b) A. Albert _et al._ (HAWC Collaboration), ApJ 911, L27 (2021b), arXiv:2101.07895 [astro-ph.HE] .
  * Sudoh  _et al._ (2021) T. Sudoh, T. Linden, and D. Hooper, J. Cosmology Astropart. Phys 2021, 010 (2021), arXiv:2101.11026 [astro-ph.HE] .
  * Breuhaus  _et al._ (2022a) M. Breuhaus, B. Reville, and J. A. Hinton, A&A 660, A8 (2022a), arXiv:2109.05296 [astro-ph.HE] .
  * de Oña Wilhelmi  _et al._ (2022) E. de Oña Wilhelmi, R. López-Coto, E. Amato, and F. Aharonian, ApJ 930, L2 (2022), arXiv:2204.09440 [astro-ph.HE] .
  * Joshi  _et al._ (2022) J. C. Joshi, S. J. Tanaka, L. S. Miranda, and S. Razzaque, arXiv e-prints , arXiv:2205.00521 (2022), arXiv:2205.00521 [astro-ph.HE] .
  * Fiori  _et al._ (2022) M. Fiori, B. Olmi, E. Amato, R. Bandiera, N. Bucciantini, L. Zampieri, and A. Burtovoi, MNRAS 511, 1439 (2022), arXiv:2201.02221 [astro-ph.HE] .
  * Martin  _et al._ (2022) P. Martin, L. Tibaldo, A. Marcowith, and S. Abdollahi, arXiv e-prints , arXiv:2207.11178 (2022), arXiv:2207.11178 [astro-ph.HE] .
  * Amenomori  _et al._ (2021c) M. Amenomori _et al._ (Tibet AS$\gamma$ Collaboration), Phys. Rev. Lett. 126, 141101 (2021c), arXiv:2104.05181 [astro-ph.HE] .
  * Lipari and Vernetto (2018) P. Lipari and S. Vernetto, Phys. Rev. D 98, 043003 (2018), arXiv:1804.10116 [astro-ph.HE] .
  * De la Torre Luque  _et al._ (2022) P. De la Torre Luque, D. Gaggero, D. Grasso, O. Fornieri, K. Egberts, C. Steppa, and C. Evoli, arXiv e-prints (2022), arXiv:2203.15759 [astro-ph.HE] .
  * Berezhko and Völk (2000) E. G. Berezhko and H. J. Völk, ApJ 540, 923 (2000), arXiv:astro-ph/0004353 [astro-ph] .
  * Prodanović  _et al._ (2007) T. Prodanović, B. D. Fields, and J. F. Beacom, Astroparticle Physics 27, 10 (2007), arXiv:astro-ph/0603618 [astro-ph] .
  * Casanova and Dingus (2008) S. Casanova and B. L. Dingus, Astroparticle Physics 29, 63 (2008), arXiv:0711.2753 [astro-ph] .
  * Ahlers and Murase (2014) M. Ahlers and K. Murase, Phys. Rev. D 90, 023010 (2014), arXiv:1309.4077 [astro-ph.HE] .
  * Ahlers  _et al._ (2016) M. Ahlers, Y. Bai, V. Barger, and R. Lu, Phys. Rev. D 93, 013009 (2016), arXiv:1505.03156 [hep-ph] .
  * Linden and Buckman (2018) T. Linden and B. J. Buckman, Phys. Rev. Lett. 120, 121101 (2018), arXiv:1707.01905 [astro-ph.HE] .
  * Fang and Murase (2021) K. Fang and K. Murase, ApJ 919, 93 (2021), arXiv:2104.09491 [astro-ph.HE] .
  * Liu and Wang (2021) R.-Y. Liu and X.-Y. Wang, ApJ 914, L7 (2021), arXiv:2104.05609 [astro-ph.HE] .
  * Vecchiotti  _et al._ (2022) V. Vecchiotti, F. Zuccarini, F. L. Villante, and G. Pagliaroli, ApJ 928, 19 (2022), arXiv:2107.14584 [astro-ph.HE] .
  * Murase  _et al._ (2014) K. Murase, T. A. Thompson, and E. O. Ofek, MNRAS 440, 2528 (2014), arXiv:1311.6778 [astro-ph.HE] .
  * Cristofari  _et al._ (2020) P. Cristofari, P. Blasi, and E. Amato, Astroparticle Physics 123, 102492 (2020), arXiv:2007.04294 [astro-ph.HE] .
  * Stecker (1970) F. W. Stecker, Ap&SS 6, 377 (1970).
  * Gaisser and Schaefer (1992) T. K. Gaisser and R. K. Schaefer, ApJ 394, 174 (1992).
  * Mori (2009) M. Mori, Astroparticle Physics 31, 341 (2009), arXiv:0903.3260 [astro-ph.HE] .
  * Kachelriess  _et al._ (2014) M. Kachelriess, I. V. Moskalenko, and S. S. Ostapchenko, ApJ 789, 136 (2014), arXiv:1406.0035 [astro-ph.HE] .
  * Peron and Aharonian (2022) G. Peron and F. Aharonian, A&A 659, A57 (2022), arXiv:2110.08778 [astro-ph.HE] .
  * Breuhaus  _et al._ (2022b) M. Breuhaus, J. A. Hinton, V. Joshi, B. Reville, and H. Schoorlemmer, A&A 661, A72 (2022b), arXiv:2201.03984 [astro-ph.HE] .
  * Moskalenko  _et al._ (2006) I. V. Moskalenko, T. A. Porter, and A. W. Strong, ApJ 640, L155 (2006), arXiv:astro-ph/0511149 [astro-ph] .
  * Cao  _et al._ (2021c) Z. Cao _et al._ (Lhaaso Collaboration), Science 373, 425 (2021c), arXiv:2111.06545 [astro-ph.HE] .
  * Hinton and Hofmann (2009) J. A. Hinton and W. Hofmann, ARA&A 47, 523 (2009), arXiv:1006.5210 [astro-ph.HE] .
  * Aharonian  _et al._ (2002) F. A. Aharonian, A. A. Belyanin, E. V. Derishev, V. V. Kocharovsky, and V. V. Kocharovsky, Phys. Rev. D 66, 023005 (2002), arXiv:astro-ph/0202229 [astro-ph] .
  * Aharonian  _et al._ (1994) F. A. Aharonian, L. O. Drury, and H. J. Voelk, A&A 285, 645 (1994).
  * Bykov  _et al._ (2000) A. M. Bykov, R. A. Chevalier, D. C. Ellison, and Y. A. Uvarov, ApJ 538, 203 (2000), arXiv:astro-ph/0003235 [astro-ph] .
  * Uchiyama  _et al._ (2010) Y. Uchiyama, R. D. Blandford, S. Funk, H. Tajima, and T. Tanaka, ApJ 723, L122 (2010), arXiv:1008.1840 [astro-ph.HE] .
  * Zirakashvili and Aharonian (2010) V. N. Zirakashvili and F. A. Aharonian, ApJ 708, 965 (2010), arXiv:0909.2285 [astro-ph.HE] .
  * Caprioli (2011) D. Caprioli, J. Cosmology Astropart. Phys 2011, 026 (2011), arXiv:1103.2624 [astro-ph.HE] .
  * Inoue  _et al._ (2012) T. Inoue, R. Yamazaki, S.-i. Inutsuka, and Y. Fukui, ApJ 744, 71 (2012), arXiv:1106.3080 [astro-ph.HE] .
  * Aharonian (2013) F. A. Aharonian, Astroparticle Physics 43, 71 (2013).
  * Celli  _et al._ (2019a) S. Celli, G. Morlino, S. Gabici, and F. A. Aharonian, MNRAS 487, 3199 (2019a), arXiv:1804.10579 [astro-ph.HE] .
  * Baade and Zwicky (1934) W. Baade and F. Zwicky, Proceedings of the National Academy of Science 20, 259 (1934).
  * Sano  _et al._ (2020) H. Sano _et al._ , ApJ 904, L24 (2020), arXiv:2011.05018 [astro-ph.HE] .
  * Sano and Fukui (2021) H. Sano and Y. Fukui, Ap&SS 366, 58 (2021), arXiv:2106.00708 [astro-ph.HE] .
  * Yasuda  _et al._ (2021) H. Yasuda, S.-H. Lee, and K. Maeda, ApJ 919, L16 (2021), arXiv:2109.04032 [astro-ph.HE] .
  * Yasuda  _et al._ (2022) H. Yasuda, S.-H. Lee, and K. Maeda, ApJ 925, 193 (2022), arXiv:2111.09534 [astro-ph.HE] .
  * Aharonian and Atoyan (1996) F. A. Aharonian and A. M. Atoyan, A&A 309, 917 (1996).
  * Gabici and Aharonian (2007) S. Gabici and F. A. Aharonian, ApJ 665, L131 (2007), arXiv:0705.3011 [astro-ph] .
  * Aharonian  _et al._ (2008) F. Aharonian _et al._ , A&A 481, 401 (2008), arXiv:0801.3555 [astro-ph] .
  * Gabici  _et al._ (2009) S. Gabici, F. A. Aharonian, and S. Casanova, MNRAS 396, 1629 (2009), arXiv:0901.4549 [astro-ph.HE] .
  * Fujita  _et al._ (2009) Y. Fujita, Y. Ohira, S. J. Tanaka, and F. Takahara, ApJ 707, L179 (2009), arXiv:0911.4482 [astro-ph.HE] .
  * Ohira  _et al._ (2011) Y. Ohira, K. Murase, and R. Yamazaki, MNRAS 410, 1577 (2011), arXiv:1007.4869 [astro-ph.HE] .
  * Uchiyama  _et al._ (2012) Y. Uchiyama, S. Funk, H. Katagiri, J. Katsuta, M. Lemoine-Goumard, H. Tajima, T. Tanaka, and D. F. Torres, ApJ 749, L35 (2012), arXiv:1203.3234 [astro-ph.HE] .
  * Mitchell  _et al._ (2021) A. M. W. Mitchell, G. P. Rowell, S. Celli, and S. Einecke, MNRAS 503, 3522 (2021), arXiv:2103.01787 [astro-ph.HE] .
  * Morlino and Caprioli (2012) G. Morlino and D. Caprioli, A&A 538, A81 (2012), arXiv:1105.6342 [astro-ph.HE] .
  * Ackermann  _et al._ (2013) M. Ackermann _et al._ , Science 339, 807 (2013), arXiv:1302.3307 [astro-ph.HE] .
  * Jogler and Funk (2016) T. Jogler and S. Funk, ApJ 816, 100 (2016).
  * Abdalla  _et al._ (2018c) H. Abdalla _et al._ (H. E. S. S. Collaboration), A&A 612, A3 (2018c), arXiv:1802.05172 [astro-ph.HE] .
  * Drury (2011) L. O. Drury, MNRAS 415, 1807 (2011), arXiv:1009.4799 [astro-ph.GA] .
  * D’Angelo  _et al._ (2016) M. D’Angelo, P. Blasi, and E. Amato, Phys. Rev. D 94, 083003 (2016), arXiv:1512.05000 [astro-ph.HE] .
  * Nava  _et al._ (2016) L. Nava, S. Gabici, A. Marcowith, G. Morlino, and V. S. Ptuskin, MNRAS 461, 3552 (2016), arXiv:1606.06902 [astro-ph.HE] .
  * Gaggero  _et al._ (2018) D. Gaggero, F. Zandanel, P. Cristofari, and S. Gabici, MNRAS 475, 5237 (2018), arXiv:1710.05038 [astro-ph.HE] .
  * Nava  _et al._ (2019) L. Nava, S. Recchia, S. Gabici, A. Marcowith, L. Brahimi, and V. Ptuskin, MNRAS 484, 2684 (2019), arXiv:1903.03193 [astro-ph.HE] .
  * Celli  _et al._ (2019b) S. Celli, G. Morlino, S. Gabici, and F. A. Aharonian, MNRAS 490, 4317 (2019b), arXiv:1906.09454 [astro-ph.HE] .
  * Brose  _et al._ (2020) R. Brose, M. Pohl, I. Sushch, O. Petruk, and T. Kuzyo, A&A 634, A59 (2020), arXiv:1909.08484 [astro-ph.HE] .
  * Brose  _et al._ (2021) R. Brose, M. Pohl, and I. Sushch, A&A 654, A139 (2021), arXiv:2108.10773 [astro-ph.HE] .
  * Kamijima and Ohira (2021) S. F. Kamijima and Y. Ohira, Phys. Rev. D 104, 083028 (2021), arXiv:2107.05820 [astro-ph.HE] .
  * Suzuki  _et al._ (2022) H. Suzuki, A. Bamba, R. Yamazaki, and Y. Ohira, ApJ 924, 45 (2022), arXiv:2110.13304 [astro-ph.HE] .
  * Ptuskin and Zirakashvili (2003) V. S. Ptuskin and V. N. Zirakashvili, A&A 403, 1 (2003), arXiv:astro-ph/0302053 [astro-ph] .
  * Schure and Bell (2013) K. M. Schure and A. R. Bell, MNRAS 435, 1174 (2013), arXiv:1307.6575 [astro-ph.HE] .
  * Cardillo  _et al._ (2015) M. Cardillo, E. Amato, and P. Blasi, Astroparticle Physics 69, 1 (2015), arXiv:1503.03001 [astro-ph.HE] .
  * Marcowith  _et al._ (2018) A. Marcowith, V. V. Dwarkadas, M. Renaud, V. Tatischeff, and G. Giacinti, MNRAS 479, 4470 (2018), arXiv:1806.09700 [astro-ph.HE] .
  * Inoue  _et al._ (2021) T. Inoue, A. Marcowith, G. Giacinti, A. Jan van Marle, and S. Nishino, ApJ 922, 7 (2021), arXiv:2108.13433 [astro-ph.HE] .
  * Aharonian  _et al._ (2007) F. Aharonian _et al._ , ApJ 661, 236 (2007), arXiv:astro-ph/0612495 [astro-ph] .
  * Fukui  _et al._ (2012) Y. Fukui _et al._ , ApJ 746, 82 (2012), arXiv:1107.0508 [astro-ph.GA] .
  * Fukui  _et al._ (2021) Y. Fukui, H. Sano, Y. Yamane, T. Hayakawa, T. Inoue, K. Tachihara, G. Rowell, and S. Einecke, ApJ 915, 84 (2021), arXiv:2105.02734 [astro-ph.HE] .
  * Fukui  _et al._ (2017) Y. Fukui _et al._ , ApJ 850, 71 (2017), arXiv:1708.07911 [astro-ph.HE] .
  * Fukuda  _et al._ (2014) T. Fukuda, S. Yoshiike, H. Sano, K. Torii, H. Yamamoto, F. Acero, and Y. Fukui, ApJ 788, 94 (2014), arXiv:1405.2599 [astro-ph.HE] .
  * Sano  _et al._ (2019) H. Sano _et al._ , ApJ 876, 37 (2019), arXiv:1805.10647 [astro-ph.HE] .
  * Fujita  _et al._ (2022) Y. Fujita, R. Yamazaki, and Y. Ohira, arXiv e-prints , arXiv:2205.12276 (2022), arXiv:2205.12276 [astro-ph.HE] .
  * Dermer and Powale (2013) C. D. Dermer and G. Powale, A&A 553, A34 (2013), arXiv:1210.8071 [astro-ph.HE] .
  * Zeng  _et al._ (2019) H. Zeng, Y. Xin, and S. Liu, ApJ 874, 50 (2019), arXiv:1811.12644 [astro-ph.HE] .
  * Suzuki  _et al._ (2020) H. Suzuki, A. Bamba, R. Yamazaki, and Y. Ohira, PASJ 72, 72 (2020), arXiv:2006.03382 [astro-ph.HE] .
  * Bykov  _et al._ (2020) A. M. Bykov, A. Marcowith, E. Amato, M. E. Kalyashova, J. M. D. Kruijssen, and E. Waxman, Space Sci. Rev. 216, 42 (2020), arXiv:2003.11534 [astro-ph.HE] .
  * Aharonian  _et al._ (2019) F. Aharonian, R. Yang, and E. de Oña Wilhelmi, Nature Astronomy 3, 561 (2019), arXiv:1804.02331 [astro-ph.HE] .
  * Yang  _et al._ (2018) R.-z. Yang, E. de Oña Wilhelmi, and F. Aharonian, A&A 611, A77 (2018), arXiv:1710.02803 [astro-ph.HE] .
  * Sun  _et al._ (2020) X.-N. Sun, R.-Z. Yang, Y.-F. Liang, F.-K. Peng, H.-M. Zhang, X.-Y. Wang, and F. Aharonian, A&A 639, A80 (2020), arXiv:2006.00879 [astro-ph.HE] .
  * Higdon and Lingenfelter (2005) J. C. Higdon and R. E. Lingenfelter, ApJ 628, 738 (2005).
  * Gupta  _et al._ (2020) S. Gupta, B. B. Nath, P. Sharma, and D. Eichler, MNRAS 493, 3159 (2020), arXiv:1910.10168 [astro-ph.HE] .
  * Bhadra  _et al._ (2022) S. Bhadra, S. Gupta, B. B. Nath, and P. Sharma, MNRAS 510, 5579 (2022), arXiv:2201.00529 [astro-ph.HE] .
  * Vieu  _et al._ (2022a) T. Vieu, S. Gabici, V. Tatischeff, and S. Ravikularaman, MNRAS 512, 1275 (2022a), arXiv:2201.07488 [astro-ph.HE] .
  * Vieu  _et al._ (2022b) T. Vieu, B. Reville, and F. Aharonian, arXiv e-prints , arXiv:2207.01432 (2022b), arXiv:2207.01432 [astro-ph.HE] .
  * Sartore  _et al._ (2010) N. Sartore, E. Ripamonti, A. Treves, and R. Turolla, A&A 510, A23 (2010), arXiv:0908.3182 [astro-ph.GA] .
  * Faucher-Giguère and Kaspi (2006) C.-A. Faucher-Giguère and V. M. Kaspi, ApJ 643, 332 (2006), arXiv:astro-ph/0512585 [astro-ph] .
  * Aharonian  _et al._ (2021) F. Aharonian _et al._ (Lhaaso Collaboration), Chinese Physics C 45, 025002 (2021), arXiv:2010.06205 [astro-ph.HE] .
  * Aartsen  _et al._ (2017) M. G. Aartsen _et al._ (IceCube Collaboration), ApJ 849, 67 (2017), arXiv:1707.03416 [astro-ph.HE] .
  * Gaggero  _et al._ (2015) D. Gaggero, D. Grasso, A. Marinelli, A. Urbano, and M. Valli, ApJ 815, L25 (2015), arXiv:1504.00227 [astro-ph.HE] .
  * Gaggero  _et al._ (2017) D. Gaggero, D. Grasso, A. Marinelli, M. Taoso, and A. Urbano, Phys. Rev. Lett. 119, 031101 (2017), arXiv:1702.01124 [astro-ph.HE] .
  * Albert  _et al._ (2020c) A. Albert _et al._ (HAWC Collaboration), ApJ 905, 76 (2020c), arXiv:2007.08582 [astro-ph.HE] .
  * Aiello  _et al._ (2019) S. Aiello _et al._ , Astroparticle Physics 111, 100 (2019), arXiv:1810.08499 [astro-ph.HE] .
  * Seo  _et al._ (2018) J. Seo, H. Kang, and D. Ryu, Journal of Korean Astronomical Society 51, 37 (2018), arXiv:1804.07486 [astro-ph.HE] .
  * Gupta  _et al._ (2018) S. Gupta, B. B. Nath, and P. Sharma, MNRAS 479, 5220 (2018), arXiv:1804.05877 [astro-ph.GA] .
  * Beacom and Candia (2004) J. F. Beacom and J. Candia, J. Cosmology Astropart. Phys 2004, 009 (2004), arXiv:hep-ph/0409046 [hep-ph] .
  * Acharya  _et al._ (2019) B. S. Acharya _et al._ (Cherenkov Telescope Array Consortium),  _Science with the Cherenkov Telescope Array_ (2019).
  * Albert  _et al._ (2019) A. Albert _et al._ , arXiv e-prints , arXiv:1902.08429 (2019), arXiv:1902.08429 [astro-ph.HE] .
  * Kato  _et al._ (2021) S. Kato _et al._ , Experimental Astronomy 52, 85 (2021), arXiv:2109.13490 [astro-ph.HE] .
  * Strong and Moskalenko (1998) A. W. Strong and I. V. Moskalenko, ApJ 509, 212 (1998), arXiv:astro-ph/9807150 [astro-ph] .
  * Strong  _et al._ (2000) A. W. Strong, I. V. Moskalenko, and O. Reimer, ApJ 537, 763 (2000), arXiv:astro-ph/9811296 [astro-ph] .
  * Moskalenko  _et al._ (2002) I. V. Moskalenko, A. W. Strong, J. F. Ormes, and M. S. Potgieter, ApJ 565, 280 (2002), arXiv:astro-ph/0106567 [astro-ph] .
  * Strong  _et al._ (2004) A. W. Strong, I. V. Moskalenko, and O. Reimer, ApJ 613, 962 (2004), arXiv:astro-ph/0406254 [astro-ph] .
  * Ptuskin  _et al._ (2006) V. S. Ptuskin, I. V. Moskalenko, F. C. Jones, A. W. Strong, and V. N. Zirakashvili, ApJ 642, 902 (2006), arXiv:astro-ph/0510335 [astro-ph] .
  * Porter  _et al._ (2017) T. A. Porter, G. Jóhannesson, and I. V. Moskalenko, ApJ 846, 67 (2017), arXiv:1708.00816 [astro-ph.HE] .
  * Porter  _et al._ (2019) T. A. Porter, G. Jóhannesson, and I. V. Moskalenko, ApJ 887, 250 (2019), arXiv:1909.02223 [astro-ph.HE] .
  * Boschini  _et al._ (2020) M. J. Boschini _et al._ , ApJ 889, 167 (2020), arXiv:1911.03108 [astro-ph.HE] .
  * Porter  _et al._ (2021) T. A. Porter, G. Johannesson, and I. V. Moskalenko, arXiv e-prints , arXiv:2112.12745 (2021), arXiv:2112.12745 [astro-ph.HE] .
  * Gaggero  _et al._ (2014) D. Gaggero, L. Maccione, D. Grasso, G. Di Bernardo, and C. Evoli, Phys. Rev. D 89, 083007 (2014), arXiv:1311.5575 [astro-ph.HE] .
  * Evoli  _et al._ (2017) C. Evoli, D. Gaggero, A. Vittino, G. Di Bernardo, M. Di Mauro, A. Ligorini, P. Ullio, and D. Grasso, J. Cosmology Astropart. Phys 2017, 015 (2017), arXiv:1607.07886 [astro-ph.HE] .
  * Evoli  _et al._ (2018) C. Evoli, D. Gaggero, A. Vittino, M. Di Mauro, D. Grasso, and M. N. Mazziotta, J. Cosmology Astropart. Phys 2018, 006 (2018), arXiv:1711.09616 [astro-ph.HE] .
  * Maurin  _et al._ (2001) D. Maurin, F. Donato, R. Taillet, and P. Salati, ApJ 555, 585 (2001), arXiv:astro-ph/0101231 [astro-ph] .
  * Putze  _et al._ (2009) A. Putze, L. Derome, D. Maurin, L. Perotto, and R. Taillet, A&A 497, 991 (2009), arXiv:0808.2437 [astro-ph] .
  * Putze  _et al._ (2010) A. Putze, L. Derome, and D. Maurin, A&A 516, A66 (2010), arXiv:1001.0551 [astro-ph.HE] .
  * Maurin (2020) D. Maurin, Computer Physics Communications 247, 106942 (2020), arXiv:1807.02968 [astro-ph.IM] .
  * Kissmann (2014) R. Kissmann, Astroparticle Physics 55, 37 (2014), arXiv:1401.4035 [astro-ph.HE] .
  * Kissmann  _et al._ (2015) R. Kissmann, M. Werner, O. Reimer, and A. W. Strong, Astroparticle Physics 70, 39 (2015), arXiv:1504.08249 [astro-ph.HE] .
  * Krumholz  _et al._ (2022) M. R. Krumholz, R. M. Crocker, and M. L. Sampson, arXiv e-prints , arXiv:2207.13838 (2022), arXiv:2207.13838 [astro-ph.HE] .
  * Hunter (2007) J. D. Hunter, Computing in Science & Engineering 9, 90 (2007).
  * van der Walt  _et al._ (2011) S. van der Walt, S. C. Colbert, and G. Varoquaux, Computing in Science Engineering 13, 22 (2011).
  * Ohio Supercomputer Center (1987) Ohio Supercomputer Center, “Ohio supercomputer center,”  (1987).