# Multi-messenger Study of Galactic Diffuse Emission with LHAASO and IceCube
Observations

Chengyu Shao  School of Physics and Astronomy, Sun Yat-Sen University, No.2
Daxue Rd, 519082, Zhuhai China  Sujie Lin<EMAIL_ADDRESS>School of
Physics and Astronomy, Sun Yat-Sen University, No.2 Daxue Rd, 519082, Zhuhai
China  Lili Yang<EMAIL_ADDRESS>School of Physics and Astronomy,
Sun Yat-Sen University, No.2 Daxue Rd, 519082, Zhuhai China  Centre for Astro-
Particle Physics, University of Johannesburg, PO Box 524, Auckland Park 2006,
South Africa

###### Abstract

With the breakthrough in PeV gamma-ray astronomy brought by the LHAASO
experiment, the high-energy sky is getting richer than before. Lately, LHAASO
Collaboration reported the observation of a gamma-ray diffuse emission with
energy up to the PeV level from both the inner and outer Galactic plane. In
these spectra, there is one bump that is hard to explain by the conventional
cosmic-ray transport scenarios. Therefore, we introduce two extra components
corresponding to unresolved sources with exponential-cutoff-power-law (ECPL)
spectral shape, one with an index of 2.4, and 20 TeV cutoff energy, and
another with index of 2.3 and 2 PeV cutoff energy. With our constructed model,
we simulate the Galactic diffuse neutrino flux and find our results are in
full agreement with the latest IceCube Galactic plane search. We estimate the
Galactic neutrino contributes of $\sim 9\%$ of astrophysical neutrinos at 20
TeV. In the high-energy regime, as expected most of the neutrinos observed by
IceCube should be from extragalactic environments.

††preprint: APS/123-QED

##  I Introduction

The origin of cosmic rays (CRs) is one of the key questions in astrophysics.
The CR energy spectrum shows the knee and ankle features. It is generally
believed that CRs with energies below the spectral knee at $\sim 10^{15}$ eV,
mainly come from our Galaxy, so-called Galactic cosmic rays (GCRs). While
those with energies above the spectral ankle at $\sim 10^{18}$ eV are mostly
from extra-galactic energetic sources. Most CR particles may lose their
directional information due to the deflection and interaction with extra-
galactic and Galactic magnetic fields and medium during their propagation.
This additional uncertainty means that the origin of CRs near the knee remains
mysterious.

To resolve the puzzles, alternative methods have been adopted. Collisions
between energetic CRs and ambient and interstellar medium generate neutral
($\pi^{0}$) and charged pions ($\pi^{\pm}$), which decay to gamma rays and
neutrinos. These secondary products detected on Earth will encode details of
both the CR and target populations. The accurate interpretation of such
measurements can provide direct information on the propagation and sources of
CRs.

In the last few decades, progress has been made in detecting high-energy
gamma-ray and neutrino emissions. The continuum diffuse gamma-ray emission has
been well measured by Fermi Large Area Telescope (LAT) up to a few hundred GeV
Ackermann et al. (2012, 2015). Later on, in the TeV energy regime, Milagro
Abdo et al. (2008), ARGO-YBJ Bartoli et al. (2015), H.E.S.S. Aharonian et al.
(2008); Abdalla et al. (2018a); Abdalla et al. (2018b) and HAWC Aab et al.
(2019); Albert et al. (2020) have been contributing data in the Galactic
plane. These measurements have only recently reached to PeV range thanks to
the Tibet AS$\gamma$ and LHAASO Amenomori et al. (2021); Cao et al. (2021,
2023). This discovery suggests the existence of PeVatrons Sudoh and Beacom
(2023), which are sources capable of accelerating particles up to PeV
energies. It is absolutely a big step towards understanding cosmic-ray physics
by exploring the knee region of the CR spectrum.

On the other hand, since the first detection of astrophysical neutrinos in
2012, IceCube has been accumulating neutrino data for more than 10 years
Aartsen et al. (2017, 2019a); Albert et al. (2018). With the development of
machine learning techniques and more statistics, the neutrino emission from
the Galactic plane has been identified lately Abbasi et al. (2023).

All these achievements can provide a hint to the injection, distribution and
propagation of CRs in our Galaxy. However, the analysis of the Galactic
diffuse emission (GDE) can be seriously contaminated by unresolved Galactic
point sources which may have a distribution similar to the interstellar gas.
This brings a challenge to recognize the accelerator of CRs. Previously, a few
groups have performed the studies on the diffuse emission from TeV to PeV,
where they discussed the possibility of the Galactic diffuse gamma-ray and
neutrino emission coming from cosmic-ray interaction, known sources, and
unresolved sources Luque et al. (2023); Egberts et al. (2023); De La Torre
Luque et al. (2023); Schwefer et al. (2023); Gabici et al. (2019).

In this work, based on the current cosmic-ray and Fermi-LAT data, and the most
recent LHAASO and IceCube Galactic plane observation, we apply the popular
GALPROP code Strong and Moskalenko (1998) to model CR transport and generate
simulated spectra and maps of the diffuse gamma-ray and neutrino emissions.
Specifically, we adopt a Diffusion plus Reacceleration (DR) model, and employ
DR-high and DR-low models to take into account the uncertainties of the
measurements obtained from the ground-based air-shower experiments, IceTop and
KASCADE respectively. However, we find a tension between the predicted gamma-
ray flux with our constructed models and the observations. To illustrate the
characteristic of the LHAASO Galactic plane spectrum and explain the excess,
we invent two populations of Galactic sources (EXTRA1 and EXTRA2) with
exponential-cutoff-power-law (ECPL) spectra shape. In the energy up to
$10^{5}~{}\mathrm{GeV}$, the spectrum of EXTRA1 has an index of 2.40 and
$\sim$ 20 TeV cutoff energy. In the higher energy end up to
$10^{6}~{}\mathrm{GeV}$, another component EXTRA2 with an index of 2.3 and a 2
PeV cutoff is introduced. This can be naturally explained by the two types of
unresolved sources in our Galaxy with different maximum cosmic-ray energy.

Based on the constructed models, we also estimate the diffuse Galactic
neutrino flux that is consistent with the latest IceCube Galactic plane
search. We found the Galactic neutrinos can contribute around 9 percent to the
all-sky neutrino events at 20 TeV. At PeV energy, most of the neutrinos are
coming from outside of our Galaxy. However, due to a few uncertain factors
like the mechanisms, numbers, and distribution of these unresolved sources in
our Galaxy and limited observation capability, there is still some space
allowed for modeling. Therefore, to further reveal the puzzles, the next-
generation Imaging Air Cherenkov Telescopes (IACTs) and neutrino detectors
with a larger effective area and better angular and energy resolution, which
can dramatically provide a precise location and morphology of sources, are in
high demand.

The paper is organized as follows. Section II provides the description of the
multi-messenger data, including cosmic-ray, gamma-ray, and neutrino
observations used in this work. In Section III, we present the injection and
propagation models of cosmic rays, together with the addition of extra source
components for fitting the gamma-ray data. Based on the constructed models, we
show the calculated Galactic diffuse gamma-ray and neutrino emission in
Section IV. In Section V, we give a discussion about the obtained results and
the origins of the two extra components. In Section VI, we give a summary and
future outlook for the multi-messenger observations.

##  II Multi-messenger observation

Thanks to the development of both satellites and ground-based observatories,
diffuse high-energy neutrinos with energies between 10 TeV to PeV Aartsen et
al. (2013), ultra-high-energy cosmic rays (UHECRs, $>10^{18}$ eV) Pie (2015),
and high-energy gamma rays from MeV to PeV have been measured, or upper limits
have been provided Atwood et al. (2009); Aab et al. (2019); Cao et al. (2023).
As there is a natural connection between these three messengers, neutrinos and
gamma rays are produced during the CR propagation and can directly point back
to the origin of CRs. Their joint detection and analysis should be a very
efficient way to explore the Universe. Aartsen et al. (2018); Abbott et al.
(2017). Moreover, the energy budgets of UHECRs, PeV neutrinos, and isotropic
sub-TeV gamma rays are comparable Fang and Murase (2018), which supports the
unification of high-energy cosmic particles.

Before exploring, understanding, and identifying the mechanisms and physical
processes of the astrophysical sources of CRs, the diffuse backgrounds
originating from our Galaxy should be seriously studied. One accurate diffuse
template can provide great help in analyzing the upcoming data. For this
purpose, we attempt to constrain the diffuse emission with current
observation. The measurements used in this work are presented below.

###  II.1 High-energy cosmic rays

The high-energy CR particles are accelerated by energetic astrophysical
sources like supernova remnants (SNRs) and propagate inside the Galactic
magnetic field around the Galactic disk after escaping. Although only the CR
fluxes around the sun could be measured, their distribution throughout the
Galaxy can be predicted by the propagation model. Generally, the propagation
model is constrained by the secondary-to-primary flux ratio observation, such
as B/C Collaboration (2022) and $\mathrm{{}^{10}Be/^{9}Be}$ Aguilar et al.
(2016). More details regarding the propagation model can be found in Section
III.

Their fluxes around the earth have been directly measured by space-born
experiments like AMS-02 Aguilar et al. (2015a, 2016, 2017, 2019a, 2019b) and
DArk Matter Particle Explorer (DAMPE) Alemanno et al. (2021); An et al.
(2019); Ambrosi et al. (2017); Collaboration (2022), and also indirectly
measured by the ground-based experiments like IceTop Aartsen et al. (2019b)
and KASCADE Antoni et al. (2005).

One has to notice that the measurements of the energies of the knees disagree
between IceTop and KASCADE, as shown in Figure 1. As the KASCADE experiment
uses the QGSJET-II-02 model while IceTop uses the Sibyll 2.1 model instead,
the discrepancies are caused by the large systematic uncertainty of the
hadronic model. In our study, we refer to the models derived from KASCADE and
IceTop data as DR-low and DR-high respectively.

(a) proton

(b) helium

Figure 1: Best-fitting spectra of protons (top), and helium nuclei (bottom),
along with the observation data from IceTop (blue circles), KASCADE (yellow
crosses), AMS-02 (red triangles), and DAMPE (purple stars). The solid line
represents the DR-low model, while the dashed line represents the DR-high
model.

In this work, to estimate the diffuse gamma-ray and neutrino emission, the
proton and electron plus positron spectra observed by Voyager, AMS-02, IceTop,
and KASCADE are adopted to constrain the Galactic CR distribution as seen in
Figure 1 and Figure 2.

Figure 2:  The black solid line shows the fitted electrons plus positrons
spectrum, and the measurements from AMS-02 (red dots) and DAMPE (yellow
crosses) are also marked.

###  II.2 Gamma-ray sky

The diffuse gamma-ray emission has been well measured by a few satellites
below TeV energies, such as EGRET, and followed by Fermi-LAT Hunter et al.
(1997); Atwood et al. (2009). Recently, the Galactic plane has been observed
up to 1 PeV, thanks to the Tibet AS$\gamma$ and LHAASO Amenomori et al.
(2021); Cao et al. (2023) experiments. These discoveries show the evidence of
hadronic origin of sub-PeV diffuse gamma rays, which are generated during the
propagation of tens of PeV CRs.

The LHAASO experiment announced the source-subtracted galactic diffuse gamma-
ray fluxes from the inner galactic plane (15${}^{\circ}<$l$<125{{}^{\circ}}$,
$|b|<5{{}^{\circ}}$) and outer plane (125${{}^{\circ}}<$l$<235{{}^{\circ}}$,
$|b|<5{{}^{\circ}}$) for the first time lately. Where a simple power law is
adopted to describe the spectra for both regions with similar spectral indices
of $-2.99$, which is consistent with the CR spectral index of the knee region.

In Figure 3, the data for the window of 15${}^{\circ}<$l$<125{{}^{\circ}}$,
$|b|<5{{}^{\circ}}$ from LHAASO and Fermi-LAT experiments Zhang et al. (2023),
and for 25${}^{\circ}<$l$<100{{}^{\circ}}$, $|b|<5{{}^{\circ}}$ from Tibet
AS$\gamma$ are presented. As can be seen, both the LHAASO and Tibet AS$\gamma$
data are in agreement with the Fermi-LAT data. However, the result from LHAASO
is a few times lower than that from AS$\gamma$. This is due to the different
analysis methods of these two experiments. Where LHAASO analyzes the data by
masking sources included in the TeVCAT with a radius of five times the
Gaussian extension widths. Therefore, this cut procedure may lose a large part
data of the Galactic plane, where the diffuse CR and unresolved sources are
located.

Figure 3: The gamma-ray data from Fermi-LAT (black crosses) and LHAASO
experiments (blue triangles) in the region of
15${}^{\circ}<$l$<125{{}^{\circ}}$, $|b|<5{{}^{\circ}}$, gamma-ray data from
Tibet AS$\gamma$ (red plus), IceCube total $\nu$ (blue shaded region) and
their results with $\pi^{0}$ model in the Galactic plane (red shaded region)
are shown.

###  II.3 Neutrino sky

Since the first observation of the astrophysical neutrino signal in the TeV -
PeV energy range in 2012 Aartsen et al. (2013), IceCube has kept updating the
neutrino sky for more than 10 years. The event distribution is consistent with
being isotropic, and the origin of these neutrino signals is still unresolved.
With larger statistics, IceCube has recently shown that there are more events
at lower Galactic latitudes and a deficit of neutrino events at high Galactic
latitudes Aartsen et al. (2019a). The IceCube Neutrino Observatory has
provided 6-year all-sky total and 10-year Galactic plane data Aartsen et al.
(2019a); Abbasi et al. (2023). In this recently updated data sample for
Galactic neutrino search, they performed the analysis for the cascade events
with lower energy thresholds. The neutrino emission from the Galactic plane is
reported at the 4.5$\sigma$ level of significance Abbasi et al. (2023) with a
total of 59,592 events selected over the entire sky in the energy range of 500
GeV to several PeV. As shown in Figure 3, the best-fitting Galactic plane
neutrino flux is comparable with the gamma-ray flux.

The total neutrino observation as shown with blue shaded region in Figure 3
includes events from Galactic and extra-galactic diffuse backgrounds and
astrophysical sources, whose spectrum follows a simple power-law distribution
as can be found below,

$\Phi_{\nu}=\Phi_{0}(\frac{E}{100\mathrm{TeV}})^{-\gamma}.$ (1)

here the normalization factor $\Phi_{0}$ is 1.66$\times
10^{-18}\mathrm{GeV^{-1}cm^{-2}s^{-1}sr^{-1}}$, and the common spectral index
$\gamma$ is 2.53. The observed neutrino spectrum is softer than $E^{-2}$ which
is comparable with the observed diffuse extragalactic gamma-ray background
Bechtol et al. (2017).

##  III Models

Both the diffuse gamma-ray and neutrino flux are generated by CR particles
when they propagate in the Galaxy. The hadronic component of CRs can induce
gamma-ray and neutrino emission through proton-proton interaction, as well as
bremsstrahlung radiation. On the other hand, the leptonic component of CRs
contributes to gamma-ray emissions through the inverse Compton (IC) effect. To
calculate these processes, along with the propagation effect of CRs, we
utilize the well-established GALPROP code Strong and Moskalenko (1998) in this
study. In this section, we provide a detailed description of the CR
propagation and emission model that are employed.

###  III.1 Cosmic-ray propagation

In the propagation model, CRs are assumed to undergo diffusion within the
Galactic magnetic field, taking into account the possible effects such as
reacceleration, energy loss, fragmentation, and decay. The diffusion
coefficient is parameterized as
$D(R)=\beta^{\eta}D_{0}(R/4\mathrm{GV})^{\delta}$, where $D_{0}$ is the
normalization factor at a reference rigidity of 4 GV, R is the particle’s
rigidity, $\beta$ is the velocity of the particles in natural units, $\delta$
is the slope of rigidity dependence, and $\eta$ is a phenomenological
parameter introduced to fit the low energy secondary-to-primary ratios.
Besides the diffusion, the convection or reacceleration effect is also
required by the observed $B/C$ data.

In some recent studies, with more secondary CR species like Li, Be, and B
precisely measured by AMS-02, it was found that the reacceleration effect is
favored Yuan et al. (2020). In this work, we adopt a Diffusion plus
Reacceleration (DR) model as a benchmark model. The model parameters that
describe the propagation processes are adopted following the work of other
groups Yuan et al. (2020), corresponding to the best-fit values obtained by
fitting the Li, Be, B, C, and O measurements from AMS-02. As listed in Table
1, the half height of diffuse zone $z_{h}$ is 6.3 kpc, and the Alfven speed
$v_{A}$ that describes the strength of reacceleration is
$33.76~{}\mathrm{km~{}s^{-1}}$.

Table 1: Propagation parameters. $D_{0}$ | $\delta$ | $z_{h}$ | $v_{A}$ | $\eta$ |
---|---|---|---|---|---
| ($10^{28}$ cm2 s-1) | (kpc) | (km s-1) |  |
$7.69$ | $0.362$ | $6.3$ | $33.76$ | $-0.05$ |

###  III.2 Cosmic-Ray injection

It is widely accepted that SNRs are the most promising galactic high-energy CR
sources, whose shock provides the ideal environment for first-order Fermi
accelerations of relativistic particles. In the discovery of 12 Galactic PeV
accelerators by LHAASO, eight of them are somehow linked to SNRs Cao et al.
(2021). Therefore we make a simple assumption that CRs are injected into the
Galaxy by the SNRs. As it is not possible to gather information on all
historical SNRs, for an estimation, we employ a continuous source distribution
for SNRs as follows

$f(r,z)=\left(\frac{r}{r_{\odot}}\right)^{1.25}\exp\left[-\frac{3.56(r-r_{\odot})}{r_{\odot}}\right]\exp\left(-\frac{|z|}{z_{s}}\right),$
(2)

where $r_{\odot}=8.3\mathrm{kpc}$ is the distance of the sun,
$z_{s}=0.2\mathrm{kpc}$ is a scale factor that indicates the thickness of the
Galactic disk.

Given the many spectral structures revealed by recent direct detection
experiments, the injection spectra of CRs may be quite complicated. A
multiple-broken-power-law spectrum is employed to describe these features as
seen below,

$f(x)=\left\\{\begin{aligned} &R^{-\nu_{0}}e^{-\frac{R}{Rc}},&R<R_{1}\\\
&(\prod_{i=1}^{n}R_{i}^{\nu_{i}-\nu_{i-1}})R^{-\nu_{n}}e^{-\frac{R}{Rc}},&R_{n}\leq
R<R_{n+1}\\\
&(\prod_{i=1}^{4}R_{i}^{\nu_{i}-\nu_{i-1}})R^{-\nu_{4}}e^{-\frac{R}{Rc}},&R_{4}\leq
R\end{aligned}\right.$ (3)

where n in the second row is from 1 to 3. The corresponding observation could
constrain the injection parameters in Equation 3 for different CR species. As
there exist discrepancies between the IceTop and KASCADE measurements, we
construct two models, differing in their injections, to indicate the upper and
lower boundary of theoretical estimation, so-called DR-high and DR-low models
respectively.

Following Ref. Zhang et al. (2023), the spectral structures are assumed to be
mainly due to the source injection as in Equation 3, without taking into
account the change of propagation parameters. The constraints for proton,
helium, and electron plus positron were performed. These three kinds of CR
particles play a major role in contributing to Galactic $\gamma$ and neutrino
emission. We develop our neutrino and gamma-ray sky map based on these best-
fitting parameters, listed in Table 2 and 3. The comparisons between
observation and model are shown in Figure 1 and 2.

Table 2: Source injection and solar modulation parameters as in Equation 3 for proton and Helium nuclei. | Proton | Helium
---|---|---
| DR-high | DR-low | DR-high | DR-low
| (IceTop) | (KASCADE) | (IceTop) | (KASCADE)
$\nu_{0}$ | $2.06$ | $2.06$ | $1.46$ | $1.46$
$\nu_{1}$ | $2.43$ | $2.43$ | $2.36$ | $2.36$
$\nu_{2}$ | $2.22$ | $2.22$ | $2.12$ | $2.12$
$\nu_{3}$ | $2.52$ | $2.52$ | $2.42$ | $2.42$
$\nu_{4}$ | $2.18$ | $2.32$ | $2.08$ | $2.28$
$R_{1}/{\rm GV}$ | $13.9$ | $13.9$ | $1.99$ | $1.99$
$R_{2}/{\rm TV}$ | $0.50$ | $0.50$ | $0.65$ | $0.65$
$R_{3}/{\rm TV}$ | $15.0$ | $15.0$ | $15.0$ | $15.0$
$R_{4}/{\rm TV}$ | $100.0$ | $100.0$ | $100.0$ | $100.0$
$R_{c}/{\rm PV}$ | $12.0$ | $4.0$ | $6.0$ | $4.0$
$\Phi/{\rm GV}$ | $0.700$ | $0.700$ | $0.700$ | $0.700$
Table 3: Source injection and solar modulation parameters of electron plus positron. $\nu_{0}^{-}$ | $\nu_{1}^{-}$ | $\nu_{2}^{-}$ | $\nu_{3}^{-}$ |  $R_{1}^{-}$/GV |  $R_{2}^{-}$/GV |  $R_{3}^{-}$/GV |  $R_{c}^{-}$/TV |  $\Phi^{-}$/GV
---|---|---|---|---|---|---|---|---
$2.33$ | $0.01$ | $2.88$ | $2.45$ | $0.950$ | $4.19$ | $55.7$ | $6.27$ | 1.1
$c_{e^{+}}$ | $\nu_{1}^{+}$ | $\nu_{2}^{+}$ |  $R_{1}^{+}$/GV |  $R_{c}^{+}$/TV |  $\Phi^{+}$/GV |  |  |
$1.00$ | $3.04$ | $2.08$ | $31.2$ | $3.42$ | 1.1 |  |  |

As seen in Table 2, for both DR-high and -low model, most of the best-fitting
parameters $\nu_{i}$ and $R_{i}$ are identical to each other, except $R_{c}$
and $\nu_{4}$. Here $R_{c}$ represents the characteristic cutoff rigidity of
the exponential cutoff spectral, describing the knee energy of those
particles. Apparently, the $R_{c}$ of DR-high is higher than that of DR-low,
because of the different knee energies from IceTop and KASCADE. On the other
hand, the $\nu_{4}$ of DR-high is smaller than that of DR-low, which is due to
a harder spectrum from the IceTop measurement.

###  III.3 Gamma-Ray expectation

With the propagation and injection of CRs fixed, we analyze the gamma-ray sky
map. We apply the GALPROP code to calculate the diffuse emissions from a few
processes, including natural pion decay, bremsstrahlung, and inverse Compton
scattering (ICS). The AAfrag package Kachelrieß et al. (2019) is adopted to
estimate the secondary gamma-ray and neutrino production from inelastic
hadronic interactions.

We show the diffuse gamma-ray spectra measured by the LHAASO and Fermi-LAT
experiments, along with our model predictions for both the inner region
(Figure 4(a) and 4(c)) and outer region (Figure 4(b) and 4(d)). To ensure a
self-consistent comparison, we apply the same masks as in the LHAASO analysis
Cao et al. (2023) for all calculated results and data.

Compared with the gamma-ray data from Fermi-LAT and LHAASO, the predicted flux
with the DR-high-only model is consistent with the data both at energies less
than a few GeV and above 60 TeV. However, between a few GeV and 60 TeV, DR-
high-only model can not explain the LHAASO data, as can be seen in Figure 4(c)
and 4(d).

This excess below 60 TeV was initially identified through the analysis of GeV
Fermi-LAT observations Ackermann et al. (2012). To account for this excess,
some studies have proposed a spatially dependent diffusion model Guo and Yuan
(2018). However, this modification of the propagation model is insufficient to
explain the data obtained by LHAASO.

In this work, we attribute this TeV excess to unresolved sources along the
Galactic plane, which are expected to be numerous and faint within the field
of view of LHAASO and Fermi-LAT. Various physical interpretations have been
previously discussed in the literature Kantzas et al. (2023); Linden and
Buckman (2018); Vecchiotti et al. (2022). Among these interpretations, the
pulsar TeV halo and pulsar wind nebulae (PWNe) have emerged as potential
candidates Linden and Buckman (2018); Vecchiotti et al. (2022). Therefore to
fit the bump $\sim\mathcal{O}(1)\mathrm{TeV}$ in the LHAASO spectrum, we
employ an ECPL component (named EXTRA1) with an index of 2.4 and a cutoff of
20 TeV to describe these unresolved sources that follow the spatial
distribution of pulsars.

The cutoff energy of the introduced EXTRA1 in this work is lower than that of
the extra component in Reference Zhang et al. (2023) (30 TeV), as we have
introduced the EXTRA2 to account for the high-energy data.

However, this component is insufficient for the DR-low case, where an
additional component is required at PeV energy. Therefore, another ECPL
component (EXTRA2) with an index of 2.3 and a cutoff at 2 PeV is introduced
for the DR-low case. These EXTRA1 and EXTRA2 components have close spectral
indices, around the average spectral index of sources in the H.E.S.S. Galactic
Plane Survey, but different cutoff energies. The cutoff energies of gamma rays
for different scenarios indicate different maximum energies of CR particles.
For instance, the leptonic (hadronic) origin, the 20 TeV gamma-ray cutoff
energy corresponds to 700 TeV (100 TeV) electron/positron (proton) cutoff
energy. This suggests that EXTRA1 and EXTRA2 likely represent at least two
distinct types of unresolved sources in the Galaxy.

Even though recent studies indicate that this excess, so as our EXTRA1
components, prefers to be leptonic origin, as strongly constrained by the
hardening of the local cosmic-ray proton spectrum observed by AMS. No source
class has been uncovered.

For EXTRA2 sources contributing at higher energies, there is no constrain from
current cosmic-ray observation. If they are of leptonic origin, such as PWNe,
the Klein-Nishina regime is dominant. Therefore a very high acceleration rate
is required and should also be higher than the electron radiative losses. This
is quite stringent. If EXTRA2 sources are TeV halos, some studies argued a
slower diffusion of the electrons in the interstellar medium is needed
Aharonian et al. (2021), which is still not understood. Therefore, the
hadronic model cannot be excluded. To confirm and deep explore the source
mechanisms of both EXTRA1 and EXTRA2, neutrino signals as the smoking gun
would provide the direct evidence for this mystery.

(a) Inner region $\&$ DR-low

(b) Outer region $\&$ DR-low

(c) Inner region $\&$ DR-high

(d) Outer region $\&$ DR-high

Figure 4:  The diffuse gamma-ray emission calculated from the DR model. The
physical radiation of ICS (green dot-dashed line), bremsstrahlung (pink dotted
line), and pion decay (blue dashed line) are shown. Two extra source
components, EXTRA1 (red dotted line) and EXTRA2 (red dot-dashed line) with
ECPL spectra are presented. Panel (a) and (b) are the spectra obtained from
the DR-low model and panel (c) and (d) for the DR-high model. The (a) and (c)
panels show the results for the inner Galaxy region of
15${}^{\circ}<$l$<125{{}^{\circ}}$, $|b|<5{{}^{\circ}}$, while the (b) and (d)
panels display the results for the outer Galaxy region of
125${{}^{\circ}}<$l$<235{{}^{\circ}}$, $|b|<5{{}^{\circ}}$.

##  IV Results

###  IV.1 Galactic diffuse gamma-ray emission

Based on the constructed model, we generate a diffuse gamma-ray emission map
which can be used as a template for future studies. This map consists of four
components: ICS, bremsstrahlung, natural pion decay, and extra source
contributions. Except for bremsstrahlung and neutral pion decay, the spatial
distributions of all these components are different from each other. In Figure
5, we show the gamma-ray energy spectrum for the region of
25${}^{\circ}<$l$<100{{}^{\circ}}$, $|b|<5{{}^{\circ}}$ without masking as an
example. In general, this spectrum is higher than that of the region
15${}^{\circ}<$l$<125{{}^{\circ}}$, which might be due to the masking effect
from the LHAASO analysis. For any other region of interest, the predicted
gamma-ray emission can be selected in the same manner, to serve as a
background template for point source analysis.

Figure 5:  The diffuse gamma-ray emission calculated from the DR-low model.
The physical radiation of ICS (green dot-dashed line), bremsstrahlung (pink
dotted line), and pion decay (blue dashed line) are shown. Two extra source
components, EXTRA1 (red dotted line) and EXTRA2 (red dot-dashed line) with
ECPL spectra are presented. This figure shows the result for the inner Galaxy
region of 25${}^{\circ}<$l$<100{{}^{\circ}}$, $|b|<5{{}^{\circ}}$

.

###  IV.2 Galactic diffuse neutrino

We show the neutrino sky map from 100 TeV to 10 PeV resulting from Section III
in Figure 6. As one can see in Figure 7(a) and Figure 7(b), our prediction for
Galactic diffuse neutrino emission for both all-sky and Galactic plane with
DR-low model are in agreement with IceCube best-fitting flux normalizations
from the data Abbasi et al. (2023). However, for the $\pi^{0}$ template of
IceCube, an extra source contribution with a hadronic origin is needed. This
appears to contradict the fact that these sources only contribute to gamma-ray
emissions and not cosmic rays.

For comparison, IceCube’s total neutrino is also shown here. Our calculated
Galactic diffuse neutrino flux shows that the contribution of Galactic
neutrinos to the total neutrino observation is around $9\%$ at 20 TeV, as seen
in Figure 7(a).

In Figure 7(b), we present a comparison of the surface brightness of one
flavor neutrino between the Galactic contribution in the disk region
($|b|<5^{\circ}$, $25^{\circ}<l<100^{\circ}$) and the total contribution
averaged over the all-sky region. This shows the distinctiveness of the
neutrino Galactic disk compared to the isotropic neutrino background. The
neutrino flux of the Galactic disk is prominent in the energy range from 10
TeV to 100 TeV and decreases significantly at higher energies. This is
constrained by the gamma-ray and cosmic-ray measurements. Our results are in
agreement with other groups’ study Abbasi et al. (2023); Kovalev et al.
(2022). The Milky Way is a source of high-energy neutrinos consistent with the
gamma-ray observation, as seen in Figure 7(a) and 5.

Figure 6:  Calculated galactic diffuse neutrino map with energies from 100 TeV
to 10 PeV. The morphology follows the gas distribution in our Galaxy.

(a) all-sky

(b) 25${}^{\circ}<$l$<100{{}^{\circ}}$, $|b|<5{{}^{\circ}}$

Figure 7:  The predicted neutrino flux per flavor from the DR-low model
compared with the IceCube total data (blue shaded region), their $\pi^{0}$
model (red shaded region), $KRA_{\gamma}^{5}$ model (region with brown solid
edge) and $KRA_{\gamma}^{50}$ model (region with brown dashed edge). Other
components including EXTRA1 (red dotted line), EXTRA2 (red dot-dashed line),
neutrino flux with DR-low model (green dashed line), and total $\nu$ flux
(black solid line) are shown. Panel (a) is for the all-sky region, and panel
(b) is in the region of 25${}^{\circ}<$l$<100{{}^{\circ}}$,
$|b|<5{{}^{\circ}}$.

In the case of DR-high-only model, as seen in Figure 8, calculated neutrino
flux is consistent with two best-fitting results with KRAγ model. At a few
PeV, the Glashow resonance is shown in the spectrum Glashow (1960). However,
to explain the results for the $\pi^{0}$ model, EXTRA1 would be necessary.

##  V Discussion

In this work, based on the most recent PeV Galactic diffuse gamma-ray
observation from LHAASO, with two sets of CR data from IceTop and KASCADE, we
construct our DR-high and DR-low models separately. For both models, we find
it is hard to explain the LHAASO Galactic plane search with conventional CR
propagation. After adding extra source contribution, the diffuse gamma-ray
emission can be well explained both by DR-high with EXTRA1 (Model 1) and DR-
low with both EXTRA1 and EXTRA2 (Model 2).

For Model 1, one extra source spectrum EXTRA1 is introduced, with a spectral
index of 2.4 and 20 TeV cutoff. For Model 2, two extra source contributions
are introduced, where one has an index of 2.4 and 20 TeV cutoff energy and
another with an index of 2.3 and 2 PeV cutoff energy. That means there could
be two populations of sources in our Galaxy with faint gamma-ray emission
which is lower than the sensitivity of our current instruments, so they have
not been identified. They follow similar CR accelerated mechanisms with close
spectral index, but various maximum CR energy.

Based on the obtained model, we simulated the Galactic diffuse neutrino flux,
obtaining the sky map as shown in Figure 6. For example with Model 2, we
estimate the Galactic contribution of the astrophysical flux is around 9 $\%$
at 20 TeV. It is uncertain if these Galactic neutrinos are from the CR
propagation or point sources because of insufficient statistical power.
Therefore, we believe the future Imaging Air Cherenkov Telescope Marinos et
al. (2022) and upgraded neutrino observatory will resolve the point sources
and precisely provide the diffuse map and reveal the origin and propagation of
cosmic rays.

The best-fitting Galactic neutrino flux from IceCube is model dependent. Where
the $\pi^{0}$ model is constrained by the Fermi MeV to GeV gamma-ray emission
and is extrapolated to TeV, where the same spatial emission profile is
assumed. While KRAγ models take into account the spatial distribution of
spectra, and cutoff energies of 5 and 50 PeV respectively. Therefore, the
$\pi^{0}$ model provides an even event distribution along the Galactic plane,
and KRAγ models give a higher neutrino flux at the Galactic center region. So
that for the interested region of 25${}^{\circ}<$l$<100{{}^{\circ}}$,
$|b|<5{{}^{\circ}}$, the $\pi^{0}$ model gives higher flux than that from KRAγ
models. On the other hand, the cosmic-ray diffuse modeling with Galprop for
the DR-low and DR-high models in this work, is not consistent with the KRA
models from the Dragon analysis. The discrepancy between all these models is
due to the low statistics and the uncertainty of the current templates.
Further accurate measurements and studies are quite essential. We summarize
the differences in Table 4.

With only the gamma-ray and cosmic-ray observation, EXTRA1 sources prefer to
be the leptonic origin, which has been discussed by a few groups. However, in
the case of IceCube best-fitting flux for the $\pi^{0}$ model which is the
only one consistent with the recent observations of 100 TeV gamma rays by the
Tibet AS$\gamma$ Amenomori et al. (2021), a population of EXTRA1 sources with
a hadronic scenario would be necessary no matter of DR-high or DR-low model,
as seen in Figure 7 and 8. It would require this kind of source to inject
fewer high-energy protons. So the identification of neutrinos can reveal the
origin of CRs, modify the CR propagation and distribution models drastically
and explore the history of our Galaxy. If EXTRA1 sources are the leptonic
origin, where no neutrino is produced. A tension exists between the predicted
diffuse Galactic neutrino flux and the IceCube results for the $\pi^{0}$
model.

Figure 8:  The predicted neutrino flux per flavor for the DR-high model
compared with the IceCube total data (blue shaded region), their $\pi^{0}$
model (red shaded region), $KRA_{\gamma}^{5}$ model (region with brown solid
edge) and $KRA_{\gamma}^{50}$ model (region with brown dashed edge). Neutrino
flux from EXTRA1 (red dotted line), DR-high model (green dashed line), and
total $\nu$ flux (black solid line) are shown for the region of
25${}^{\circ}<$l$<100{{}^{\circ}}$, $|b|<5{{}^{\circ}}$.

For the case of IceCube best-fitting fluxes for the KRAγ models, it provides
the lower limit for the neutrino emission from the Galactic plane. No other
extra hadronic scenario sources are needed. In other words, the EXTRA1 sources
would be of the leptonic origin. Which needs no extra proton injection and
releases the tension between the data and models.

No matter which results for the different model templates obtained by IceCube,
both leptonic and hadronic origins of EXTRA2 source are allowed by data. High-
energy neutrino emission is a unique diagnostic of hadronic content. With the
future PeV neutrino detection with improved sensitivity, the EXTRA2 sources
could be identified. If EXTRA2 of Model 2 is the hadronic origin, the Galactic
neutrino will contribute around $1\%$ to the total IceCube neutrino at PeV.
Otherwise, the contribution is $\sim$ $0.4\%$.

Table 4: The scenario of extra source meets the measurements of cosmic rays, gamma rays, and neutrinos. Origin | $p$ | $\nu$ | $\gamma$
---|---|---|---
Leptonic | $\checkmark$ | KRAγ | $\checkmark$
Hadronic | hadron-less | $\pi^{0}$ | $\checkmark$

##  VI Summary

In summary, thanks to the recent observations from LHAASO Cao et al. (2023)
and IceCube Abbasi et al. (2023), the Galactic diffuse sky has become richer,
especially in the high-energy regime. The LHAASO measurements show a bump in
the gamma-ray spectrum, where contributions from extra unresolved sources are
needed. The IceCube Collaboration confirms the high-energy neutrinos from the
Galactic plane. Our calculated flux with models obtained from gamma-ray
observation is consistent with the neutrino data. However for the best-fitting
results for the $\pi^{0}$ model from IceCube data, the EXTRA1 sources with
hadronic scenario is a must to fill the gap between calculated flux and data.
Even though it would be disfavored by CR measurements.

The joint analysis of cosmic rays, gamma rays, and neutrinos has shown strong
power in understanding the high-energy sky. For example, the diffuse gamma-ray
detection by LHAASO can probe the CR density in our Galaxy and solve the
problem of the disagreement between IceTop and KASCADE. Secondly, the neutrino
detection can reveal the hidden sources which are not transparent for gamma-
ray emission.

The current results from all these three messengers are in agreement with each
other. More evidence shows the existence of PeVatrons in our Galaxy. The next
step forward should be identifying the mysterious astronomical origin of high-
energy cosmic rays with upgraded neutrino and gamma-ray detectors.

##  VII Acknowledgements

We thank the referee for the useful and helpful comments and suggestions. This
work is supported by the National Natural Science Foundation of China (NSFC)
grants 12005313, 12205388, and 12261141691.

## References

  * Ackermann et al. (2012) M. Ackermann, M. Ajello, W. Atwood, L. Baldini, J. Ballet, G. Barbiellini, D. Bastieri, K. Bechtol, R. Bellazzini, B. Berenji, et al., The Astrophysical Journal 750, 3 (2012).
  * Ackermann et al. (2015) M. Ackermann, M. Ajello, A. Albert, W. Atwood, L. Baldini, J. Ballet, G. Barbiellini, D. Bastieri, K. Bechtol, R. Bellazzini, et al., The Astrophysical Journal 799, 86 (2015).
  * Abdo et al. (2008) A. A. Abdo et al., Astrophys. J. 688, 1078 (2008), eprint 0805.0417.
  * Bartoli et al. (2015) B. Bartoli et al. (ARGO-YBJ), Astrophys. J. 806, 20 (2015), eprint 1507.06758.
  * Aharonian et al. (2008) F. Aharonian, A. Akhperjanian, U. B. De Almeida, A. Bazer-Bachi, Y. Becherini, B. Behera, W. Benbow, K. Bernlöhr, C. Boisson, A. Bochow, et al., Physical Review Letters 101, 261104 (2008).
  * Abdalla et al. (2018a) H. Abdalla et al. (HESS), Astron. Astrophys. 612, A9 (2018a), eprint 1706.04535.
  * Abdalla et al. (2018b) H. Abdalla, A. Abramowski, F. Aharonian, F. A. Benkhali, E. Angüner, M. Arakawa, M. Arrieta, P. Aubert, M. Backes, A. Balzer, et al., Astronomy & Astrophysics 612, A1 (2018b).
  * Aab et al. (2019) A. Aab, P. Abreu, M. Aglietta, I. Albuquerque, J. M. Albury, I. Allekotte, A. Almela, J. A. Castillo, J. Alvarez-Muñiz, G. A. Anastasi, et al., Journal of Cosmology and Astroparticle Physics 2019, 022 (2019).
  * Albert et al. (2020) A. Albert, R. Alfaro, C. Alvarez, J. A. Camacho, J. Arteaga-Velázquez, K. Arunbabu, D. A. Rojas, H. A. Solares, V. Baghmanyan, E. Belmont-Moreno, et al., The Astrophysical Journal 905, 76 (2020).
  * Amenomori et al. (2021) M. Amenomori et al. (Tibet ASgamma), Phys. Rev. Lett. 126, 141101 (2021), eprint 2104.05181.
  * Cao et al. (2021) Z. Cao et al. (LHAASO), Nature 594, 33 (2021).
  * Cao et al. (2023) Z. Cao et al. (LHAASO) (2023), eprint 2305.05372.
  * Sudoh and Beacom (2023) T. Sudoh and J. F. Beacom, Phys. Rev. D 107, 043002 (2023), eprint 2209.03970.
  * Aartsen et al. (2017) M. G. Aartsen et al. (IceCube), Astrophys. J. 849, 67 (2017), eprint 1707.03416.
  * Aartsen et al. (2019a) M. G. Aartsen et al. (IceCube), Astrophys. J. 886, 12 (2019a), eprint 1907.06714.
  * Albert et al. (2018) A. Albert et al. (ANTARES, IceCube), Astrophys. J. Lett. 868, L20 (2018), eprint 1808.03531.
  * Abbasi et al. (2023) R. Abbasi et al. (IceCube), Science 380, adc9818 (2023), eprint 2307.04427.
  * Luque et al. (2023) P. D. l. T. Luque, D. Gaggero, D. Grasso, O. Fornieri, K. Egberts, C. Steppa, and C. Evoli, Astron. Astrophys. 672, A58 (2023), eprint 2203.15759.
  * Egberts et al. (2023) K. Egberts, C. Steppa, and K. P. Peters (2023), eprint 2303.11850.
  * De La Torre Luque et al. (2023) P. De La Torre Luque, D. Gaggero, and D. Grasso, PoS ECRS, 103 (2023).
  * Schwefer et al. (2023) G. Schwefer, P. Mertsch, and C. Wiebusch, Astrophys. J. 949, 16 (2023), eprint 2211.15607.
  * Gabici et al. (2019) S. Gabici, C. Evoli, D. Gaggero, P. Lipari, P. Mertsch, E. Orlando, A. Strong, and A. Vittino, Int. J. Mod. Phys. D 28, 1930022 (2019), eprint 1903.11584.
  * Strong and Moskalenko (1998) A. W. Strong and I. V. Moskalenko, Astrophys. J. 509, 212 (1998), eprint astro-ph/9807150.
  * Aartsen et al. (2013) M. G. Aartsen et al. (IceCube), Science 342, 1242856 (2013), eprint 1311.5238.
  * Pie (2015) _The Pierre Auger Observatory: Contributions to the 34th International Cosmic Ray Conference (ICRC 2015)_ (2015), eprint 1509.03732.
  * Atwood et al. (2009) W. B. Atwood et al. (Fermi-LAT), Astrophys. J. 697, 1071 (2009), eprint 0902.1089.
  * Aartsen et al. (2018) M. G. Aartsen et al. (IceCube, Fermi-LAT, MAGIC, AGILE, ASAS-SN, HAWC, H.E.S.S., INTEGRAL, Kanata, Kiso, Kapteyn, Liverpool Telescope, Subaru, Swift NuSTAR, VERITAS, VLA/17B-403), Science 361, eaat1378 (2018), eprint 1807.08816.
  * Abbott et al. (2017) B. P. Abbott et al. (LIGO Scientific, Virgo, Fermi-GBM, INTEGRAL), Astrophys. J. Lett. 848, L13 (2017), eprint 1710.05834.
  * Fang and Murase (2018) K. Fang and K. Murase, Nature Phys. 14, 396 (2018), eprint 1704.00015.
  * Collaboration (2022) D. Collaboration, Science Bulletin 67, 2162 (2022), ISSN 20959273, eprint 2210.08833.
  * Aguilar et al. (2016) M. Aguilar et al. (AMS), Phys. Rev. Lett. 117, 231102 (2016).
  * Aguilar et al. (2015a) M. Aguilar et al. (AMS), Phys. Rev. Lett. 114, 171103 (2015a).
  * Aguilar et al. (2017) M. Aguilar et al. (AMS), Phys. Rev. Lett. 119, 251101 (2017).
  * Aguilar et al. (2019a) M. Aguilar et al. (AMS), Phys. Rev. Lett. 122, 101101 (2019a).
  * Aguilar et al. (2019b) M. Aguilar et al. (AMS), Phys. Rev. Lett. 122, 041102 (2019b).
  * Alemanno et al. (2021) F. Alemanno et al., Phys. Rev. Lett. 126, 201102 (2021), eprint 2105.09073.
  * An et al. (2019) Q. An, R. Asfandiyarov, P. Azzarello, P. Bernardini, X. J. Bi, M. S. Cai, J. Chang, D. Y. Chen, H. F. Chen, J. L. Chen, et al., Science Advances 5, eaax3793 (2019), ISSN 2375-2548, arXiv:1909.12860 [astro-ph], URL http://arxiv.org/abs/1909.12860.
  * Ambrosi et al. (2017) G. Ambrosi, Q. An, R. Asfandiyarov, P. Azzarello, P. Bernardini, B. Bertucci, M. S. Cai, J. Chang, D. Y. Chen, H. F. Chen, et al., Nature 552, 63 (2017), ISSN 0028-0836, 1476-4687, arXiv:1711.10981 [astro-ph, physics:hep-ex, physics:hep-ph], URL http://arxiv.org/abs/1711.10981.
  * Aartsen et al. (2019b) M. G. Aartsen et al. (IceCube), Phys. Rev. D 100, 082002 (2019b), eprint 1906.04317.
  * Antoni et al. (2005) T. Antoni et al. (KASCADE), Astropart. Phys. 24, 1 (2005), eprint astro-ph/0505413.
  * Hunter et al. (1997) S. D. Hunter et al., Astrophys. J. 481, 205 (1997).
  * Zhang et al. (2023) R. Zhang, X. Huang, Z.-H. Xu, S. Zhao, and Q. Yuan, _Galactic diffuse gamma-ray emission from GeV to PeV energies in light of up-to-date cosmic ray measurements_ (2023), eprint 2305.06948.
  * Bechtol et al. (2017) K. Bechtol, M. Ahlers, M. Di Mauro, M. Ajello, and J. Vandenbroucke, The Astrophysical Journal 836, 47 (2017).
  * Yuan et al. (2020) Q. Yuan, C.-R. Zhu, X.-J. Bi, and D.-M. Wei, JCAP 11, 027 (2020).
  * Kachelrieß et al. (2019) M. Kachelrieß, I. V. Moskalenko, and S. Ostapchenko, Comput. Phys. Commun. 245, 106846 (2019), eprint 1904.05129.
  * Guo and Yuan (2018) Y.-Q. Guo and Q. Yuan, Physical Review D 97, 063008 (2018), ISSN 2470-0010, 2470-0029, eprint 1801.05904.
  * Kantzas et al. (2023) D. Kantzas, S. Markoff, A. J. Cooper, D. Gaggero, M. Petropoulou, and P. D. L. T. Luque, _Possible contribution of X-ray binary jets to the Galactic cosmic ray and neutrino flux_ (2023), eprint 2306.12715.
  * Linden and Buckman (2018) T. Linden and B. J. Buckman, Physical Review Letters 120, 121101 (2018), ISSN 0031-9007, 1079-7114, eprint 1707.01905.
  * Vecchiotti et al. (2022) V. Vecchiotti, F. Zuccarini, F. L. Villante, and G. Pagliaroli, The Astrophysical Journal 928, 19 (2022), ISSN 0004-637X, 1538-4357, eprint 2107.14584.
  * Aharonian et al. (2021) F. Aharonian, Q. An, Axikegu, L. X. Bai, Y. X. Bai, Y. W. Bao, D. Bastieri, X. J. Bi, Y. J. Bi, H. Cai, et al. (LHAASO Collaboration), Phys. Rev. Lett. 126, 241103 (2021), URL https://link.aps.org/doi/10.1103/PhysRevLett.126.241103.
  * Kovalev et al. (2022) Y. Y. Kovalev, A. V. Plavin, and S. V. Troitsky, Astrophys. J. Lett. 940, L41 (2022), eprint 2208.08423.
  * Glashow (1960) S. L. Glashow, Phys. Rev. 118, 316 (1960).
  * Marinos et al. (2022) P. D. Marinos, G. P. Rowell, T. A. Porter, and G. Jóhannesson, Mon. Not. Roy. Astron. Soc. 518, 5036 (2022), eprint 2211.01619.
  * Aguilar et al. (2015b) M. Aguilar et al. (AMS), Phys. Rev. Lett. 115, 211101 (2015b).
  * Abbasi et al. (2021) R. Abbasi et al. (IceCube), Phys. Rev. D 104, 022002 (2021), eprint 2011.03545.
  * Ackermann et al. (2016) M. Ackermann et al. (Fermi-LAT), Phys. Rev. Lett. 116, 151105 (2016), eprint 1511.00693.
  * Abbasi et al. (2022) R. Abbasi et al. (IceCube), Phys. Rev. D 106, 022005 (2022), eprint 2111.10169.
  * Evoli et al. (2012) C. Evoli, I. Cholis, D. Grasso, L. Maccione, and P. Ullio, Phys. Rev. D 85, 123511 (2012), eprint 1108.0664.