# Neutrino Observations of LHAASO Sources: Present Constraints and Future
Prospects

Tian-Qi Huang,1,2 and Zhuo Li,1,2
1Department of Astronomy, School of Physics, Peking University, Beijing
100871, People’s Republic of China
2Kavli Institute for Astronomy and Astrophysics, Peking University, Beijing
100871, People’s Republic of China
E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The Large High Altitude Air Shower Observatory (LHAASO) observed a dozen of
gamma-ray sources with significant emission above 100 TeV, which may be strong
candidates of PeVatrons. Neutrino observations are crucial to diagnose whether
the gamma-ray radiative process is hadronic or leptonic. We use the Bayesian
method to analyze the ten-year (2008-2018) IceCube data, and hence constrain
the hadronic gamma-ray emission in the LHAASO sources. The present neutrino
data show that the hadronic gamma-ray flux from Crab Nebula is lower than the
observed gamma-ray flux at the 90% C.L. and contributes less than 86%, which
disfavors the hadronic origin of the gamma-rays below tens of TeV. For the
other LHAASO sources, the present neutrino observations cannot put useful
constraints on the gamma-ray radiative process. We consider the uncertainty of
the source extension: the upper limits on hadronic gamma-ray flux tend to
increase with the extension; and some sources, i.e., LHAASO J2032+4102, LHAASO
J1929+1745, and LHAASO J1908+0621, show relatively high statistical
significance of neutrino signals if the extension is $\lesssim 0.6^{\circ}$.
We finally estimate the future observational results of LHAASO sources by the
proposed neutrino telescopes. If the LHAASO-observed PeV gamma-rays are of
hadronic origin, Crab Nebula may be detected at $>100$ TeV at $3\sigma$ C.L.
within 20 years by a neutrino detector with the effective area 30 times that
of IceCube.

###### keywords:

neutrinos – gamma-rays:general – cosmic rays – methods: statistical

††pubyear: 2021††pagerange: Neutrino Observations of LHAASO Sources: Present
Constraints and Future Prospects–Neutrino Observations of LHAASO Sources:
Present Constraints and Future Prospects

##  1 Introduction

The origins of high-energy cosmic-rays are still unclear. Their propagation
trajectories are deflected from the direction of accelerators by the magnetic
field in the medium. Thus it is difficult to identify the sources by cosmic
ray observations. However, cosmic-rays interact with the background baryons or
radiation in/around the source and generate high energy gamma-rays (via
$\pi^{0}$ decay) and neutrinos (via $\pi^{\pm}$ decay), which are not affected
by the magnetic field in propagation and can be good indicators of cosmic-ray
sources.

The most energetic photon ever observed is with 1.4 PeV and from the direction
of Cygnus OB2. It was detected by the $\rm km^{2}$ array (KM2A) of the Large
High Altitude Air Shower Observatory (LHAASO) (Cao et al., 2021a). LHAASO-KM2A
had also detected significant gamma-ray emission above 100 TeV from 12 sources
(Cao et al., 2021a; The LHAASO Collaboration et al., 2021; Cao et al., 2021b),
which are probably the Galactic accelerators of PeV cosmic-rays, namely,
PeVatrons. According to gamma-ray observations, various PeVatron candidates
had been proposed, including, e.g., the Galactic center (e.g. HESS
Collaboration et al., 2016), supernova remnants (SNRs, e.g. Tibet AS$\gamma$
Collaboration et al., 2021), pulsar wind nebulae (PWNe, e.g. Arons, 2012) and
young massive star clusters (YMCs, e.g. Aharonian et al., 2019).

A detection of gamma-rays up to hundreds of TeV is not sufficient to identify
a source as PeVatron, because the gamma-rays can be produced by leptonic
processes as well, e.g. inverse Compton scattering off electrons. However, the
detection of neutrinos accompanying the $\pi^{0}$ decayed gamma-rays will
definitely help to identify PeVatrons.

Most of the 12 LHAASO sources with 100 TeV emission (LHAASO sources hereafter)
show TeV gamma-ray counterparts. The neutrino emission from these TeV gamma-
ray counterparts have been investigated in previous works, e.g., Kappes et al.
(2009); Halzen et al. (2017). In some cases, sophisticated models have been
proposed for the neutrino emission, for example, the Crab Nebula (e.g. Amato
et al., 2003), Cygnus Region (e.g. Anchordoqui et al., 2007; Bykov et al.,
2021) and SNR G106.3+2.7 (e.g. Ge et al., 2021).

IceCube Neutrino Observatory, a $\rm km^{3}$ scale detector at the South Pole,
has carried out searches for time-integrated neutrino signals from the TeV
gamma-ray sources, but results in non-detection so far (e.g. Abbasi et al.,
2011; Aartsen et al., 2013, 2017a, 2019, 2020). In the recent search using the
ten-year data of IceCube from 2008 to 2018, Aartsen et al. (2020) searched for
astrophysical neutrinos from the directions of a list of candidate sources,
where MGRO J1908+06 is the most significant Galactic source but still far from
the $3\sigma$ confidence level. The upper limits on the neutrino flux from
these candidate sources are provided, under the assumptions of point-like
neutrino sources and $E^{-2}$ or $E^{-3}$ neutrino spectra. Using these
neutrino upper limits Huang & Li (2022) put constraints on the hadronic gamma-
ray components from LHAASO sources.

Since the upper limits on neutrino flux depend on the assumptions of the
source extensions and neutrino spectral shapes, in this work we evaluate the
effect of these assumptions on the upper limit estimate. We will use the
Bayesian method and the ten-year IceCube muon-track data (IceCube
Collaboration et al., 2021) from 2008 to 2018 to search for the neutrino flux
from LHAASO sources and hence constrain the hadronic components in gamma-rays.
The current operating detectors (e.g. IceCube and Baikal-GVD Zaborov, 2020)
only shows non-detection, but several neutrino telescopes are under
construction (e.g. KM3NeT-ARCA Adrián-Martínez et al., 2016) or being planed
(e.g. P-ONE and IceCube-Gen2 Agostini et al., 2020; Aartsen et al., 2021),
thus in this work we further predict the future observation results of
neutrinos from the LHAASO sources.

This paper is organized as follows. In section 2, we introduce the IceCube
muon-track data that is used. In section 3, we describe the Bayesian method
used for deriving the neutrino upper limits. In section 4, we present the
connection between the spectra of neutrinos and the associated hadronic gamma-
rays. In section 5, we compare the gamma-ray observations with the upper
limits of hadronic gamma-ray flux and discuss the uncertainties in the upper
limit estimate. In section 6, we estimate the significance of neutrino signals
from LHAASO sources in the combined search by the future projects. Finally
section 7 is the summary and discussion.

Figure 1: An example of the background and signal PDFs. Left: the spatial
background PDF. Different colors represent different data samples. Right: the
energy background PDF and the energy signal PDF in the case of MGRO J1908+06
for data sample IC86-II. The thin (thick) purple line shows the energy signal
PDF for an $E^{-2}$ ($E^{-3}$) neutrino spectrum. The black line shows the
energy background PDF at $\delta_{s}=6.18^{\circ}$.

##  2 The muon-track data from IceCube

IceCube released the muon-track data from 2008 to 2018 which was used in the
ten-year search for point-like neutrino sources (IceCube Collaboration et al.,
2021). We will use these data to search for neutrino signals from the LHAASO
sources. The data is composed of three parts:

The experimental data events. The experimental data events are grouped into
five samples including IC40, IC59, IC79, IC86-I, and IC86-II, corresponding to
different construction levels of the detector. The number in the sample name
represents the number of strings in the detector. Digital Optical Modules on
the string record the Cherenkov light from the charged particles produced in
the neutrino-nucleus interactions ($\nu+N$). The interaction time, the
reconstructed direction and the reconstructed energy ($E_{\rm rec}$) are given
for each event. Note that $E_{\rm rec}$ is the reconstructed energy of the
muon passing through the detector, different from the energy of the incident
neutrino which generates the muon through the neutrino-nucleus charge current
(CC) interaction $\nu_{\mu}+N{\xrightarrow{\rm CC}}\mu+X$.

The instrument response functions. The instrument response functions include
the effective area and the smearing matrix. The effective area $A_{\rm
eff}(E_{\nu},\delta_{\nu})$ relies on the neutrino energy $E_{\nu}$ and the
declination angle $\delta_{\nu}$. The smearing matrix $\mathcal{M}(E_{\rm
rec}|E_{\nu},\delta_{\nu})$ gives the fractional count of simulated events in
the reconstructed energy bin relative to all events in the ($E_{\nu}$,
$\delta_{\nu}$) bin. The matrix tells us the probability to get the
reconstructed energy $E_{\rm rec}$ when a neutrino with the energy $E_{\nu}$
enters the detector from the declination $\delta_{\nu}$.

The detector uptime. The detector uptime records the time periods during which
the detector is running well. We can get the livetime of the detector for each
data sample.

##  3 Bayesian method for signal search

Since the present data only result in the non-detection of neutrino signals
from LHAASO sources (see below), we should find a way to give the upper limits
to the neutrino flux. We will use the Bayesian method to search for signal
events. The Bayesian method here is based on the form of unbinned likelihood
function which is widely used in the neutrino source searches (Braun et al.,
2008; Abbasi et al., 2011; Aartsen et al., 2014, 2017b, 2019, 2020; Kochocki
et al., 2021; Zhou et al., 2021). Given the likelihood function one can
calculate the probability density of neutrino event number, and then derive
the 90% C.L. upper limit of the neutrino flux
$\Phi^{90\%}_{\nu_{\mu}+\bar{\nu}_{\mu}}$.

###  3.1 Likelihood function

The first step is to build the likelihood function. The likelihood $L$ for
observing $n_{s}$ signal events from a source is given by the product of
probability density functions (PDFs) for each track event. Consider a source
at position $\vec{x}_{s}=(\alpha_{s},\delta_{s})$, with $\alpha_{s}$ the right
ascension and $\delta_{s}$ the declination angle, with angular extension
$\sigma_{s}$. Assuming a power law neutrino spectrum $dN_{\nu}/dE_{\nu}\propto
E_{\nu}^{-\gamma}$, a likelihood can be built for the source,

$\displaystyle L$
$\displaystyle=\prod_{k}\prod_{i}^{N_{k}}\left(\frac{n_{s}^{k}}{N_{k}}S_{i}+\frac{N_{k}-n^{k}_{s}}{N_{k}}B_{i}\right),$
(1)

where $S_{i}$ and $B_{i}$ are the PDFs for signal and background events,
$k\in\left\\{{\rm
IC40,\,IC59,\,IC79,\,\text{IC86-I},\,\text{IC86-II}}\right\\}$ represents five
data samples, $N_{k}$ is the total number of events in data sample $k$, and
$n^{k}_{s}$ is the number of signal events in data sample $k$.

If the source extension is small, the detection probability over the entire
source extension can be considered to be proportional to the effective area
$A_{\rm eff}^{k}$ for the source declination $\delta_{s}$. We write

$n^{k}_{s}=n_{s}\times\frac{\int t_{k}\int E_{\nu}^{-\gamma}A_{\rm
eff}^{k}(E_{\nu}\,,\delta_{s})dE_{\nu}}{\sum_{k^{{}^{\prime}}}\int
t_{k^{{}^{\prime}}}\int E_{\nu}^{-\gamma}A_{\rm
eff}^{k^{{}^{\prime}}}(E_{\nu},\delta_{s})dE_{\nu}},$ (2)

where $n_{s}$ is the total number of signal events and $t_{k}$ is the detector
livetime for data sample $k$.

The PDFs for signal and background events are given, respectively, by

$\displaystyle S_{i}$ $\displaystyle=S^{\rm
spat}(\vec{x}_{i}|\sigma_{i},\vec{x}_{s},\sigma_{s})\times S^{\rm
ener}(E_{i}|\vec{x}_{s},\gamma),$ (3) $\displaystyle B_{i}$
$\displaystyle=B^{\rm spat}(\delta_{i})\times B^{\rm ener}(E_{i}|\delta_{i}).$
(4)

Here $S^{\rm spat}$ is the spatial signal PDF that describes the distribution
of the reconstructed direction of signal events $\vec{x}_{i}$. If the source
follows a 2D Gaussian with extension $\sigma_{s}$, and the event angular
uncertainty is $\sigma_{i}$, we have

$S^{\rm
spat}(\vec{x}_{i}|\sigma_{i},\vec{x}_{s},\sigma_{s})=\frac{1}{2\pi(\sigma_{i}^{2}+\sigma_{s}^{2})}{\rm
exp}\left[-\frac{\phi^{2}}{2(\sigma_{i}^{2}+\sigma_{s}^{2})}\right]\frac{\phi}{{\rm
sin}\phi},$ (5)

where $\phi$ is the angular distance between $\vec{x}_{i}$ and $\vec{x}_{s}$.
The spatial signal PDF is normalized as $\int S^{\rm
spat}(\vec{x}|\sigma_{i},\vec{x}_{s},\sigma_{s})d\Omega=1$. Typically,
$\sigma_{i}\approx 0.64^{\circ}$ (Halzen et al., 2017).

Next, $S^{\rm ener}$ is the energy signal PDF that describes the distribution
of the reconstructed energy of signal events, $E_{i}$, from the source
direction $\vec{x}_{s}$. We use the effective area $A_{\rm eff}^{k}$ and the
smearing matrix $\mathcal{M}_{k}$ to build the energy signal PDF,

$S^{\rm ener}(E_{i}|\vec{x}_{s},\gamma)=\frac{\int E_{\nu}^{-\gamma}A_{\rm
eff}^{k}(E_{\nu},\delta_{s})dE_{\nu}\cdot\mathcal{M}_{k}(E_{i}|E_{\nu},\delta_{s})}{\int
E_{\nu}^{-\gamma}A_{\rm eff}^{k}(E_{\nu},\delta_{s})dE_{\nu}}.$ (6)

Finally, $B^{\rm spat}$ denotes the spatial background PDF, and $B^{\rm ener}$
the energy background PDF. Both are nearly uniform in right ascension. We
assume that the background events follow the PDF of experimental data
scrambled in right ascension. For data sample $k$, we denote $N_{ij}^{k}$ the
number of events with the declination ${\rm sin}\,\delta\in[{\rm
sin}\delta_{i},\,{\rm sin}\delta_{i}+0.02)$ and the reconstructed energy ${\rm
log_{10}}E_{\rm rec}\in[{\rm log_{10}}E_{j},{\rm log_{10}}E_{j}+0.1)$. Thus
the spatial and energy background PDFs are given, respectively, by

$\displaystyle B^{\rm
spat}(\delta)=\frac{\sum_{j}N_{ij}^{k}}{N_{k}\times\Delta\Omega},$ (7)
$\displaystyle B^{\rm ener}(E_{\rm
rec}|\delta)=\frac{N_{ij}^{k}}{\sum_{j}N_{ij}^{k}},$ (8)

where $\Delta\Omega=2\pi\Delta{\rm sin}\,\delta=0.04\pi$.

Figure 1 gives an example of the background and signal PDFs. The left panel
shows the spatial background PDF in five data samples. The right panel shows
the energy PDF for background and signals in the direction of MGRO J1908+06
for data sample IC86-II. MGRO J1908+06 shows the most significant neutrino
signals in Galactic sources in the ten-year search (Aartsen et al., 2020).

###  3.2 Probability density of event number and the neutrino flux upper limit

As there is no preference on the signal event number, the prior distribution
$\pi(n_{s})$ is assumed to be uniform in the range $[0,N]$, where
$N=\sum_{k}N_{k}$ is the total number of events. Given the spectral index
$\gamma$ and source extension $\sigma_{s}$, the probability density of event
number can be expressed as

$p(n_{s}|X,\gamma,\sigma_{s})=\frac{L(X|n_{s},\gamma,\sigma_{s})\pi(n_{s})}{\int
L(X|n_{s},\gamma,\sigma_{s})\pi(n_{s})dn_{s}},$ (9)

where $X$ represents the observed data. The cumulative probability
$\int^{n_{90}}p(n_{s}|X,\gamma,\sigma_{s})dn_{s}=0.9$ defines the 90% upper
limit on signal event number $n_{90}$. Given $n_{90}$, the 90% upper limit on
neutrino flux $\Phi^{90\%}_{\nu_{\mu}+\bar{\nu}_{\mu}}=dN_{\nu}/dE_{\nu}$ can
be derived by solving out

$n_{\rm 90}=\sum_{k}\int dt_{k}\int
dE_{\rm\nu}\Phi_{\rm\nu_{\mu}+\bar{\nu}_{\mu}}^{90\%}A_{\rm
eff}^{k}(E_{\rm\nu},\delta_{s}).$ (10)

For an example, Figure 2 shows the cumulative probability function (CDF) of
MGRO J1908+06 under the point source (PS) and extended source (ES) hypothesis,
respectively. The extension of MGRO J1908+06 is around $0.34^{\circ}$ by
Aharonian et al. (2009).

Figure 2: The cumulative probability of the 100 TeV neutrino flux in the case
of MGRO J1908+06. Different models correspond to different spectral index
$\gamma$ and source extension, $\sigma_{s}=0.0^{\circ}(\rm PS)$ and
$0.34^{\circ}(\rm ES)$. The horizontal line marks the 90% confidence level.
The vertical lines mark the 90% upper limits on astrophysical neutrino flux
$\Phi^{90\%}_{\nu_{\mu}+\bar{\nu}_{\mu}}$ for different models.

###  3.3 Caveats

Some notes should be made here. Firstly, the spectrum term $E_{\nu}^{-\gamma}$
in the likelihood function can be replaced by the other spectral shapes. So we
can measure the astrophysical neutrino flux under any spectral shape
assumptions.

Secondly, considering a constant effective area for the entire source
extension, i.e., taken to be that of the source center, is not accurate enough
if the source extension is large. In this case, Equation 2 should be changed
to

$\displaystyle n^{k}_{s}=n_{s}\times\frac{\int t_{k}\int
E_{\nu}^{-\gamma}\sum_{j}\omega_{j}A_{\rm
eff}^{k}(E_{\nu}\,,\delta_{j})dE_{\nu}}{\sum_{k^{{}^{\prime}}}\int
t_{k^{{}^{\prime}}}\int E_{\nu}^{-\gamma}\sum_{j}\omega_{j}A_{\rm
eff}^{k^{{}^{\prime}}}(E_{\nu},\delta_{j})dE_{\nu}},$ (11)

and Equation 6 to

$S^{\rm ener}=\frac{\int E_{\nu}^{-\gamma}\sum_{j}\omega_{j}A_{\rm
eff}^{k}(E_{\nu},\delta_{j})\mathcal{M}_{k}(E_{i}|E_{\nu},\delta_{j})dE_{\nu}}{\int
E_{\nu}^{-\gamma}\sum_{j}\omega_{j}A_{\rm
eff}^{k}(E_{\nu},\delta_{j})dE_{\nu}},$ (12)

where $\omega_{j}$ is the fraction of $S^{\rm spat}$ locating in the
declination range $[\delta_{j},\delta_{j+1})$.

##  4 The neutrino and gamma-ray connection

The spectral shape of neutrino emission is required to derive the upper limit
on the hadronic gamma-ray flux. In p-p interactions, gamma-rays are produced
accompanying the neutrino production. There is a clear connection of the
spectrum and flux between these neutrinos and gamma-rays. Here we provide the
way to derive the neutrino spectrum if the accompanying gamma-ray spectrum is
known. Similar to Huang et al. (2020), we use the parameterized energy spectra
of neutrinos and gamma-rays by Kelner et al. (2006) to convert a certain
hadronic gamma-ray spectrum into the neutrino spectrum in the hadronuclear
(p-p) scenario.

The differential spectrum of neutrinos or gamma-rays
$\Phi_{\nu,\gamma}(E)\equiv dN_{\nu,\gamma}/dE$ can be given by

$\displaystyle\begin{aligned}
\Phi_{\nu,\gamma}\left(E\right)=\int_{E}^{\infty}\sigma_{\text{pp}}\left(E_{p}\right)\mathcal{G}_{p}\left(E_{p}\right)F_{\nu,\gamma}\left(\frac{E}{E_{p}},E_{p}\right)\frac{dE_{p}}{E_{p}},\end{aligned}$
(13)

where

$\mathcal{G}_{p}(E_{p})\propto\int n_{\rm
H}\left(\vec{x}\right)J_{p}\left(E_{p},\vec{x}\right)dV,$ (14)

$\sigma_{\text{pp}}$ is the cross section of inelastic p-p interactions, $V$
is the volume of the emission region, $n_{\rm H}$ is the number density of the
background gas, $J_{p}$ is the differential cosmic ray flux density, and
$F_{\rm\nu}$ ($F_{\rm\gamma}$) is the energy distribution probability of
neutrinos (gamma-rays) decayed from the secondary particles (e.g.
$\rm\pi$-mesons and $\rm\eta$-mesons) in p-p interactions. We adopt the
analytical parameterization of $F_{\gamma}$ and $F_{\nu}$ by Kelner et al.
(2006) (see Equations 58, 62, and 66 therein). The accuracy is better than 10%
in the cosmic ray energy range $100~{}{\rm GeV}\leq E_{p}\leq 100~{}{\rm PeV}$
and for neutrino or gamma-ray energy of $E/E_{p}\geq 10^{-3}$.

Given a gamma-ray spectrum $\Phi_{\gamma}(E)$, we use Equation 13 to find out
function $\mathcal{G}_{p}(E_{p})$, and hence derive the neutrino spectrum
$\Phi_{\nu}(E)$ also with Equation 13. The p-p produced neutrino flavor ratio,
after mixing in propagation, is $\nu_{e}:\nu_{\mu}:\nu_{\tau}\approx 1:1:1$ at
the Earth. So we calculate the muon neutrino flux
$\Phi_{\nu_{\mu}+\bar{\nu}_{\mu}}$ as equal to one-third of the total all-
flavor neutrino flux.

The gamma-ray spectrum is obtained by the best fit of the observational data.
We assume a group of spectral models, as well as a pulsar wind model proposed
by Amato et al. (2003). The group of spectral models are:

$\Phi_{\gamma}\propto\left\\{\begin{aligned} &E_{10}^{-\gamma_{1}}&&({\rm
PL})\\\ &E_{10}^{-\gamma_{1}-\gamma_{2}\,{\rm log_{10}}E_{10}}&&({\rm
LOGP})\\\ &E_{10}^{-\gamma_{1}}{\rm exp}\left(-E/E_{c}\right)&&({\rm ECPL})\\\
&E_{10}^{-\gamma_{1}}{\rm exp}\left(-\sqrt{E/E_{c}}\right)&&({\rm ECPL2})\\\
\end{aligned}\right.$ (15)

where $E_{10}=E/10\,{\rm TeV}$. Table 1 shows the best-fit parameters for the
gamma-ray spectra of the LHAASO sources, and the parameter values in the
pulsar wind model111The dimensionless parameter $\mu$ describes the target
density in the nebula, $\Gamma$ is the wind Lorentz factor, and $f_{p}$ is the
fraction of spin down luminosity carried by protons.. The source extension in
Table 1 is taken according to either the measurements in TeV energy range or
the prior extension of LHAASO sources (see details in Appendix A).

Note, in the spectral fitting, the gamma-ray absorption
($\gamma\gamma\rightarrow e^{+}e^{-}$) due to the interstellar radiation field
(ISRF) and cosmic microwave background (CMB) is taken into account if the
distances to the TeV counterparts of LHAASO sources are available in the
TeVCat (Wakely & Horan, 2008). The ISRF energy density is taken from Popescu
et al. (2017). The absorption is usually not important (see Table 2).

All the neutrino spectra used for signal search are calculated through
Equation 13 except that the spectrum for the pulsar wind model is directly
taken from Amato et al. (2003). With the upper limit on neutrino flux, we can
also get the upper limit on hadronic gamma-ray flux with a similar process.

We have tested in our calculation the accuracy of the solved-out function
$\mathcal{G}_{p}$ by comparing the best-fit gamma-ray spectrum and the one
derived by the solution function of $\mathcal{G}_{p}$. The difference is less
than 2% over an energy range $[E_{\gamma,1},E_{\gamma,2}]$, which is
corresponding to the neutrino energy range $[0.5E_{\gamma,1},0.5E_{\gamma,2}]$
where the central 90% neutrino events are located.

LHAASO | $\sigma_{s}$ | $\gamma_{1}$ | $\gamma_{2}$ | $E_{\rm c}$ | Model
---|---|---|---|---|---
Source | [deg] |  |  | [TeV] |
J0534+2202 | 0.0 | 2.86 | 0.20 | — | LOGP
J1825-1326 | 0.30 | 2.40 | 0.45 | — | LOGP
| 0.0 | 2.13 | — | 286 | ECPL
J1839-0545 | 0.41 | 2.26 | — | 36 | ECPL
J1843-0338 | 0.24 | 2.03 | — | 48 | ECPL
J1849-0003 | 0.09 | 1.99 | — | — | SPL
| 0.09 | 1.99 | — | 300 | ECPL
J1908+0621 | 0.52 | 2.53 | 0.30 | — | LOGP
J2018+3651 | 0.20 | 1.57 | — | 26 | ECPL
J2032+4102 | 1.8 | 2.94 | — | — | SPL
J2226+6057 | 0.24 | 2.29 | 0.33 | — | LOGP
J1929+1745 | 0.30 | 2.40 | — | — | SPL
J1956+2845 | 0.30 | 2.09 | — | — | SPL
J2108+5157 | 0.0 | 1.95 | — | 20 | ECPL2

LHAASO | $\sigma_{s}$ | $\mu$ | $\Gamma$ | $f_{p}$ | Model
---|---|---|---|---|---
Source | [deg] |  |  |  |
J0534+2202 | 0.0 | 1 | $10^{7}$ | 0.15 | Amato et al.

Table 1: The spectral parameters of hadronic gamma-rays and the source
extensions ($\sigma_{s}$) used for signal search. The extensions refer to the
source extensions of TeV gamma-ray counterparts given in Table 2. Figure 3:
The 90% C.L. upper limits on hadronic gamma-ray flux (thick lines) from LHAASO
sources in comparison with the observed gamma-ray flux (thin lines). The dark
parts of the thick lines correspond to the central 90% energy ranges in the
neutrino detection. For each LHAASO source, the first model in Table 1 is
shown in red, while the second model in blue. The blue data points are from
LHAASO observations (Cao et al., 2021a; The LHAASO Collaboration et al., 2021;
Cao et al., 2021b). The other gamma-ray data are taken from various
observations. The references for each LHAASO source are: (a) Crab Nebula by
Tibet AS$\gamma$ (Amenomori et al., 2019), H.E.S.S. (Aharonian et al., 2006)
and MAGIC (Albert et al., 2008); (b) HESS J1826-130 (H. E. S. S. Collaboration
et al., 2020), HESS J1825-137 (H. E. S. S. Collaboration et al., 2019), eHWC
J1825-134 (Abeysekara et al., 2020); (c) HESS J1841-055 (Aharonian et al.,
2008); (d) HESS J1843-033 and HESS J1844-030 (H. E. S. S. Collaboration et
al., 2018), 3HWC J1843-034 (Albert et al., 2020); (e) HESS J1849-000 (H. E. S.
S. Collaboration et al., 2018); (f) HESS J1908+063 (Aharonian et al., 2009),
eHWC J1907+063 (Abeysekara et al., 2020); (g) VER J2019+368 (Abeysekara et
al., 2018), eHWC J2019+368 (Abeysekara et al., 2020), TASG J2019+368
(Amenomori et al., 2021); (h) TeV J2032+4130 (Aliu et al., 2014), MGRO
J2031+41 (Abdo et al., 2007, 2009), 3HWC J2031-415 (Albert et al., 2020), TASG
J2032+414 (Amenomori et al., 2021), ARGO J2031+4157 (Bartoli et al., 2014);
(i) G106+2.7 by VERITAS (Acciari et al., 2009) and Tibet $\rm{AS\gamma}$
(Tibet AS$\gamma$ Collaboration et al., 2021); (j) 3HWC J1928+178 (Albert et
al., 2020); (k) 3HWC J1954+286 (Albert et al., 2020); and (l) 4FGL
J2108.0+5155 (Cao et al., 2021b).

##  5 Present constraints from IceCube data

Figure 4: The flux ratio between the 90% C.L. upper limit on hadronic gamma-
ray flux $\Phi^{90\%}_{\gamma}$ and the observed gamma-ray flux
$\Phi^{\rm{obs}}_{\gamma}$ as function of source extension $\sigma_{s}$. The
red (blue) lines correspond to the spectral shapes of hadronic gamma-rays
following the red (blue) lines in Figure 3. These results come from the
analysis of the ten-year IceCube data.

Using the above-described method for the IceCube data, we derive the upper
limits of hadronic gamma-ray flux from the LHAASO sources. We compare the 90%
C.L. upper limits on the hadronic gamma-ray flux with the observed gamma-ray
flux, as shown in Figure 3. Only Crab Nebula is well constrained; its derived
upper limit of hadronic gamma-ray flux disagrees with the best-fit observed
gamma-ray flux within $5\sigma$ uncertainty, and is smaller than 86% of the
latter. This upper limit is somewhat higher than that by Huang & Li (2022),
becasue using Equation 13 the additional gamma-rays from $\eta$-meson decay
are considered. The IceCube data can only constrain the hadronic flux from the
Crab Nebula below 50 TeV. However, although the central energy range covers
the energy of hundreds of TeV in the cases of LHAASO J1825-1326, HAWC
J1825-134, LHAASO J1849-0003, LHAASO J1929+1745, and LHAASO J1956+2845, their
upper limits are higher than the gamma-ray fluxes observed, thus not
constraining the hadronic component.

Some factors leading to uncertainties in the analysis should be noted here.
The assumptions on the extension of neutrino sources affect the measurements
of neutrinos/hadronic gamma-rays from LHAASO sources. Figure 4 shows the
dependence of the upper limit on the source extension. The extended-source
hypothesis brings higher upper limits than those with the point-source
hypothesis. The hadronic component can only be constrained in two cases under
the point-like source assumption; besides Crab Nebula, the other is LHAASO
J1908+0621, where the 90% C.L. hadronic upper limit is lower than the gamma-
ray flux within $\sim 1\sigma$ uncertainty in the central 90% energy range.
Moreover, the change of considering large extension (i.e., using Equation 11
and Equation 12) is relatively large for the sources in the southern sky, for
example, the expected number of signal events change from 5.4 to 5.7 for
LHAASO J1839-0545 (the extension $\delta_{s}=0.9^{\circ}$), from 4.1 to 3.9
for LHAASO J1843-0338 ($\delta_{s}=0.9^{\circ}$), but for LHAASO J1908+0621
($\delta_{s}=0.9^{\circ}$) only increases by 0.4%.

The statistical method also affect the estimate of upper limits. We select 20
sources (including the seven sources associated with LHAASO sources) from
Table III in the ten-year search paper (Aartsen et al., 2020) and compare
their upper limits222The upper limits for $E^{-2}$ ($E^{-3}$) spectrum are
taken from Table III (Figure 3) of Aartsen et al. (2020). with the ones given
by our Bayesian method. Our upper limits tend to be $1.4^{+0.6}_{-0.2}$
($0.9^{+0.2}_{-0.3}$) times of the upper limits given by Aartsen et al. (2020)
for the $E^{-2}$ ($E^{-3}$) neutrino spectrum.

The difference between the upper limits given by the two methods may be due to
several reasons. The first is the difference in statistical approaches. The
Frequentist approach is used in the ten-year search, while we use the Bayesian
approach. The results given by Bayesian approach relies on the uniform prior
distribution $\pi(n_{s})$. The second is the difference in the likelihood
functions. The likelihood function used in the ten-year search includes two
free parameters ($n_{s}$ and $\gamma$), while our likelihood function only one
($n_{s}$). The third is the difference in the background PDFs. The spatial and
energy background PDFs ($B^{\rm spat}$ and $B^{\rm ener}$) are precisely
parameterized with good data/Monte Carlo agreement (Aartsen et al., 2016,
2019) while ours are obtained from the data.

Furthermore, we also estimate the significance of signal events as function of
the unknown source extension (see Appendix B for detailed calculation). We
find that three LHAASO sources show relatively strong significance of signal
neutrino events for some extension values: LHAASO J2032+4102 if
$\delta_{s}=0.3^{\circ}$; LHAASO J1929+1745 if $\delta_{s}=0.6^{\circ}$; and
LHAASO J1908+0621 under PS assumption.

Figure 5: The statistical significance (p-value) of the future track-event
observation as function of observation time. The starting point (zero time) is
the present status, already with the 10 year (2008-2018) IceCube data. The red
and blue lines correspond to the spectral models in red and blue in Figure 3,
respectively. The dark and light lines represent the observation with low and
high threshold energy, $E_{\rm th}=1$ and $10\,{\rm TeV}$, respectively. In
the first panel for Crab Nebula, the gray lines present the pulsar wind model
of Amato et al. (2003) but with $E_{\rm th}=100\,{\rm TeV}$: the gray solid
line is for the case of 20 year observation by a detector 30 times larger than
IceCube at the South Pole; and the gray dashed-dotted line is for 20 year
observation by a detecor 30 times larger than IceCube at the position of Lake
Baikal. As for the red and blue lines, the solid lines present the cases of
the 10 year IceCube observation followed by 10 year $\rm{PLE\nu M\hbox{-}1}$
and then 10 year $\rm{PLE\nu M\hbox{-}2}$ observation; the dashed lines
present the cases for the sources at $\delta_{s}>-5^{\circ}$, the 10 year
IceCube observation followed by 10 year more IceCube and then 10 year IceCube-
Gen2 observation; and the dashed-dotted lines are the cases of the sources at
$\delta_{s}<-5^{\circ}$ observed by a IceCube-like detector at the position of
Lake Baikal for 20 years.

##  6 Future Prospects

In this section we evaluate the detection of neutrinos from LHAASO sources by
current and future neutrino telescopes, assuming that all the observed gamma-
ray emission is hadronic origin. In particular, we estimate the observational
time it takes before obtaining a $5\sigma$ detection of the source. We
consider the Planetary Neutrino Monitoring System ($\rm{PLE\nu M}$), proposed
by Schumacher et al. (2021), a concept of global repository of high-energy
neutrino observations by current and future neutrino telescopes. Two proposed
systems are considered: $\rm{PLE\nu M\hbox{-}1}$ consisting of IceCube and
three IceCube-like detectors, located at the latitudes of Baikal-GVD, KM3NeT-
ARCA and P-ONE; and $\rm{PLE\nu M\hbox{-}2}$ consisting of IceCube-Gen2 and
the three IceCube-like detectors as those in $\rm{PLE\nu M\hbox{-}1}$. Here by
IceCube-like we mean the effective area is similar to IceCube, and we assume
IceCube-Gen2 has an effective area 7.5 times that of IceCube.

To quantitatively describe the future observation results, we estimate the
statistical significance of observation with a p-value analytically expressed
as (The ATLAS Collaboration, 2011; Halzen et al., 2017)

$\displaystyle p_{\text{value}}$
$\displaystyle=\frac{1}{2}\left[1-\operatorname{erf}\left(\sqrt{q_{0}^{\rm
obs}/2}\right)\right],$ (16)

where

$\displaystyle q_{0}^{\rm obs}$
$\displaystyle=2\left[Y_{b}-N_{D}+N_{D}\ln\left(\frac{N_{D}}{Y_{b}}\right)\right],$
(17)

$Y_{b}$ is the expected number of background events, and $N_{D}$ is the median
of Poisson-distributed events containing both signal and background. We count
events within a solid angle $\Omega=\pi\sigma_{\rm eff}^{2}$ around the
source, where $\sigma_{\rm eff}=\sqrt{\sigma_{s}^{2}+(1.6\Delta\xi_{\rm
det})^{2}}$. The angular resolution $\Delta\xi_{\rm det}$ is around
$0.4^{\circ}$ for IceCube and $0.2^{\circ}$ or smaller for the detectors in
water, e.g., Baikal-GVD, KM3NeT and P-ONE. The angle $\Omega$ is the one that
contains roughly 72% of signal events from the source. This method has been
applied by Halzen et al. (2017) to evaluating the detection of the neutrino
emission from Galactic sources with IceCube.

When we calculate the number of detected neutrino events, we take the
following approach. First, only neutrinos with energy above some detection
threshold $E_{\rm th}$ are counted.Second, as the background events in up-
going track events are mainly induced by atmospheric muon neutrinos, the
software MCEq (Fedynitch et al., 2015) is employed to model the background
from atmospheric neutrinos, assuming the Gaisser-H4a model (Gaisser, 2012) for
primary cosmic rays and the Sibyll2.3c model (Fedynitch et al., 2019) for
hadronic interactions. Third, only the up-going track events with the zenith
angle $\theta_{z}>85^{\circ}$ are considered, because the discovery potential
for IceCube is much larger in the northern sky than in the south, by a few or
a few tens (see Figure 3 of Aartsen et al., 2020). Fourth, we consider that a
certain source is only visible if $\theta_{z}>85^{\circ}$, and the visible
time is calculated taking into account the Earth’s rotation.

The results about the prediction of the significance of future observations on
the 12 LHAASO sources are presented in Figure 5. In the calculation, we assume
that the hadronic gamma-rays follow the energy spectra in Figure 3. We also
assume the situation of future observations are: following the 10-year
(2008-2018) observation by IceCube is the operation of $\rm{PLE\nu
M\hbox{-}1}$ for 10 years, and then $\rm{PLE\nu M\hbox{-}2}$ for the next 10
years.

We see that within 20 years, five LHAASO sources can be discovered at the
level of $5\sigma$ if the threshold is $E_{\rm th}=1\,{\rm TeV}$: Crab Nebula,
LHAASO J1825-1326, LHAASO J1839-0545, LHAASO J1908+0621, and LHAASO
J2018+3651; one source can also be detected at the level of $5\sigma$ but for
$E_{\rm th}=10\,{\rm TeV}$: LHAASO J2226+6057; and two sources are possible to
be detected with $3\sigma$ C. L. within 20 years: HAWC J1825-134 and LHAASO
J1849-0003 (SPL model).

Consider the neutrino search only by IceCube and IceCube-Gen2 (shown as the
dashed lines in Figure 5), 4 sources can be discovered within 20 years, i.e.,
Crab Nebula, LHAASO J1908+0621, LHAASO J2018+3651, and LHAASO J2226+6057. The
first three are consistent with the predictions in Figure 18 of Aartsen et al.
(2021).

For the combined search by $\rm{PLE\nu M\hbox{-}1}$ and $\rm{PLE\nu
M\hbox{-}2}$, the IceCube-like detectors located at the northern hemisphere is
very helpful in observing sources around the equatorial plane or in the
southern sky ($\delta_{s}<-5^{\circ}$). For an example of LHAASO J1825-1326,
if all the gamma-rays are hadronic in origin then even one IceCube-like
detector at the latitude of Baikal-GVD can provide a $5\sigma$ detection
within about 6 years.

As for Crabe Nebula, if consider only the neutrino emission in the pulsar wind
model proposed by Amato et al. (2003), the signals are hard to be discovered.
In order to observe at the level of $3\sigma$ within 20 years, a large
neutrino detector with the effective area 30 times larger than IceCube is
required (see the gray lines in Figure 5).

##  7 Summary and discussion

We analyze the ten-year IceCube data with Bayesian approach to constrain the
hadronic gamma-ray flux from the directions of LHAASO sources, and give the
90% C.L. upper limits for both extended- and point-source hypotheses. The
results given by the Bayesian approach rely on the uniform prior distribution.
We further evaluate the combined search for high energy neutrinos from LHAASO
sources by IceCube, IceCube-Gen2 and three IceCube-like detectors located at
the latitudes of Baikal-GVD, KM3NeT-ARCA, and P-ONE, respectively.

To summarize, the main conclusions are: i) The 90% C.L. upper limit on
hadronic gamma-ray flux from Crab Nebula is 86% of total gamma-ray flux
observed, disfavoring hadronic origin of its gamma-rays; However, no
constraint can be made for gamma-ray radiation mechanism in the other LHAASO
sources by the present IceCube data; ii) Six LHAASO sources are expected to be
discovered at the level of $5\sigma$ by $\rm{PLE\nu M\hbox{-}1}$ and
$\rm{PLE\nu M\hbox{-}2}$ within 20 years, if their gamma-rays are totally
hadronic in origin; iii) An IceCube-like detector at the location of Baikal-
GVD is expected to give a $5\sigma$ discovery around 6 years if all the gamma-
rays from LHAASO J1825-1326 are hadronic in origin.

In our analysis, the three most significant sources in neutrino searches are
LHAASO J2032+4102, LHAASO J1929+1745 and LHAASO J1908+0621. Under the point-
source hypothesis, the pre-trial p-value is 0.051 for LHAASO J1908+0621 and
0.091 for LHAAS0 J2032+4102 which are similar with the p-values in the
directions of MGRO J1908+06 and 2HWC J2031+415 in the ten-year search by
IceCube (Aartsen et al., 2020). The TeV counterparts of LHAASO J1929+1745 are
not in the source list of ten-year search because they are not bright in the
TeVCat (Wakely & Horan, 2008). LHAASO J1929+1745 is even one of the faintest
LHAASO sources at 100 TeV (Cao et al., 2021a). Mori et al. (2020) tried to
explain the gamma-rays from 2HWC J1928+177 ($0.10^{\circ}$ from LHAASO
J1929+1745) with a hadronic accelerator model.

The neutrino constraint on the Crab Nebula disfavors that hadronic process
dominates the bulk of the observed gamma-ray emission. However, the gamma-ray
spectrum by LHAASO shows a possible hardening around PeV energies (The LHAASO
Collaboration et al., 2021). If the PeV gamma-rays are hadronic in origin, a
neutrino detector with the effective area $\sim 30$ times larger than IceCube
is required to observe the neutrinos associated with these PeV gamma-rays at
the level of $3\sigma$ within 20 years (see the first panel in Figure 5).

Baikal-GVD is with 8 clusters deployed and two more clusters per year are
planned from 2022 to 2024. It will occupy a water volume around $0.7\,{\rm
km^{3}}$ in total (Baikal-GVD Collaboration et al., 2021). A KM3NeT-ARCA block
comprises 115 strings and eleven detection units have been deployed (KM3NeT
Collaboration et al., 2021). P-ONE has two pathfinders deployed and the design
is underway (Bailly et al., 2021). The P-ONE Explorer (10 strings) is planned
to be deployed in 2023-2024 and remainder of the array (70 strings) is planned
for deployment between 2028-2030 (Agostini et al., 2020). IceCube-Gen2 is
planned for deployment between 2027-2033. Thus $\rm{PLE\nu M\hbox{-}1}$ will
start to operate after 2030, while $\rm{PLE\nu M\hbox{-}2}$ after 2033. The
neutrino source candidates are expected to be identified or excluded at the
high confidence level from 2030 to 2050.

## Acknowlegdments

We thank Hao Zhou and Xiao-Yuan Huang for useful discussions. This work is
supported by the Natural Science Foundation of China (No. 11773003, U1931201)
and the China Manned Space Project (CMS-CSST-2021-B11).

## Data Availability

The ten-year (2008-2018) muon-track data by IceCube are available at
https://icecube.wisc.edu/data-releases/2021/01/all-sky-point-source-icecube-
data-years-2008-2018/. The TeVCat online source catalog is available at
http://tevcat.uchicago.edu. The ISRF energy density is available at
https://cdsarc.cds.unistra.fr/viz-bin/cat/J/MNRAS/470/2539. The code MCEq is
available at https://github.com/afedynitch/MCEq. Other data are available
through the references in the captions of Figure 3 and Table 2.

## References

  * Aartsen et al. (2013) Aartsen M. G., et al., 2013, ApJ, 779, 132
  * Aartsen et al. (2014) Aartsen M. G., et al., 2014, ApJ, 796, 109
  * Aartsen et al. (2016) Aartsen M. G., et al., 2016, ApJ, 833, 3
  * Aartsen et al. (2017a) Aartsen M. G., et al., 2017a, ApJ, 835, 151
  * Aartsen et al. (2017b) Aartsen M. G., et al., 2017b, ApJ, 849, 67
  * Aartsen et al. (2019) Aartsen M. G., et al., 2019, European Physical Journal C, 79, 234
  * Aartsen et al. (2020) Aartsen M. G., et al., 2020, Phys. Rev. Lett., 124, 051103
  * Aartsen et al. (2021) Aartsen M. G., et al., 2021, Journal of Physics G Nuclear Physics, 48, 060501
  * Abbasi et al. (2011) Abbasi R., et al., 2011, ApJ, 732, 18
  * Abdo et al. (2007) Abdo A. A., et al., 2007, ApJ, 664, L91
  * Abdo et al. (2009) Abdo A. A., et al., 2009, ApJ, 700, L127
  * Abeysekara et al. (2018) Abeysekara A. U., et al., 2018, ApJ, 861, 134
  * Abeysekara et al. (2020) Abeysekara A. U., et al., 2020, Phys. Rev. Lett., 124, 021102
  * Acciari et al. (2009) Acciari V. A., et al., 2009, ApJ, 703, L6
  * Adrián-Martínez et al. (2016) Adrián-Martínez S., et al., 2016, Journal of Physics G Nuclear Physics, 43, 084001
  * Agostini et al. (2020) Agostini M., et al., 2020, Nature Astronomy, 4, 913
  * Aharonian et al. (2006) Aharonian F., et al., 2006, A&A, 457, 899
  * Aharonian et al. (2008) Aharonian F., et al., 2008, A&A, 477, 353
  * Aharonian et al. (2009) Aharonian F., et al., 2009, A&A, 499, 723
  * Aharonian et al. (2019) Aharonian F., Yang R., de Oña Wilhelmi E., 2019, Nature Astronomy, 3, 561
  * Albert et al. (2008) Albert J., et al., 2008, ApJ, 674, 1037
  * Albert et al. (2020) Albert A., et al., 2020, ApJ, 905, 76
  * Albert et al. (2021) Albert A., et al., 2021, ApJ, 907, L30
  * Aliu et al. (2014) Aliu E., et al., 2014, ApJ, 783, 16
  * Amato et al. (2003) Amato E., Guetta D., Blasi P., 2003, A&A, 402, 827
  * Amenomori et al. (2019) Amenomori M., et al., 2019, Phys. Rev. Lett., 123, 051101
  * Amenomori et al. (2021) Amenomori M., et al., 2021, Phys. Rev. Lett., 127, 031102
  * Anchordoqui et al. (2007) Anchordoqui L. A., Beacom J. F., Goldberg H., Palomares-Ruiz S., Weiler T. J., 2007, Phys. Rev. D, 75, 063001
  * Arons (2012) Arons J., 2012, Space Sci. Rev., 173, 341
  * Baikal-GVD Collaboration et al. (2021) Baikal-GVD Collaboration et al., 2021, arXiv e-prints, p. arXiv:2106.06288
  * Bailly et al. (2021) Bailly N., et al., 2021, European Physical Journal C, 81, 1071
  * Bartoli et al. (2014) Bartoli B., et al., 2014, ApJ, 790, 152
  * Braun et al. (2008) Braun J., Dumm J., De Palma F., Finley C., Karle A., Montaruli T., 2008, Astroparticle Physics, 29, 299
  * Bykov et al. (2021) Bykov A. M., Petrov A. E., Kalyashova M. E., Troitsky S. V., 2021, ApJ, 921, L10
  * Cao et al. (2021a) Cao Z., et al., 2021a, Nature, 594, 33
  * Cao et al. (2021b) Cao Z., et al., 2021b, ApJ, 919, L22
  * Cowan et al. (2011) Cowan G., Cranmer K., Gross E., Vitells O., 2011, European Physical Journal C, 71, 1554
  * Fedynitch et al. (2015) Fedynitch A., Engel R., Gaisser T. K., Riehn F., Stanev T., 2015, in European Physical Journal Web of Conferences. p. 08001 (arXiv:1503.00544), doi:10.1051/epjconf/20159908001
  * Fedynitch et al. (2019) Fedynitch A., Riehn F., Engel R., Gaisser T. K., Stanev T., 2019, Phys. Rev. D, 100, 103018
  * Fesen et al. (1997) Fesen R. A., Shull J. M., Hurford A. P., 1997, AJ, 113, 354
  * Gaisser (2012) Gaisser T. K., 2012, Astroparticle Physics, 35, 801
  * Ge et al. (2021) Ge C., Liu R.-Y., Niu S., Chen Y., Wang X.-Y., 2021, The Innovation, 2, 100118
  * H. E. S. S. Collaboration (2020) H. E. S. S. Collaboration 2020, Nature Astronomy, 4, 167
  * H. E. S. S. Collaboration et al. (2018) H. E. S. S. Collaboration et al., 2018, A&A, 612, A1
  * H. E. S. S. Collaboration et al. (2019) H. E. S. S. Collaboration et al., 2019, A&A, 621, A116
  * H. E. S. S. Collaboration et al. (2020) H. E. S. S. Collaboration et al., 2020, A&A, 644, A112
  * HESS Collaboration et al. (2016) HESS Collaboration et al., 2016, Nature, 531, 476
  * Halzen et al. (2017) Halzen F., Kheirandish A., Niro V., 2017, Astroparticle Physics, 86, 46
  * Huang & Li (2022) Huang T.-Q., Li Z., 2022, ApJ, 925, 85
  * Huang et al. (2020) Huang Y., Li Z., Wang W., Zhao X., 2020, MNRAS, 492, 4246
  * IceCube Collaboration et al. (2021) IceCube Collaboration et al., 2021, arXiv e-prints, p. arXiv:2101.09836
  * KM3NeT Collaboration et al. (2021) KM3NeT Collaboration et al. 2021, Liveblog about the ongoing sea operation for ARCA, https://www.km3net.org/blog-about-the-ongoing-sea-operation-for-arca/
  * Kappes et al. (2009) Kappes A., Halzen F., Murchadha A. Ó., 2009, Nuclear Instruments and Methods in Physics Research A, 602, 117
  * Kelner et al. (2006) Kelner S. R., Aharonian F. A., Bugayov V. V., 2006, Phys. Rev. D, 74, 034018
  * Kochocki et al. (2021) Kochocki A., Takhistov V., Kusenko A., Whitehorn N., 2021, ApJ, 914, 91
  * Mori et al. (2020) Mori K., et al., 2020, ApJ, 897, 129
  * Popescu et al. (2017) Popescu C. C., Yang R., Tuffs R. J., Natale G., Rushton M., Aharonian F., 2017, MNRAS, 470, 2539
  * Schumacher et al. (2021) Schumacher L. J., Huber M., Agostini M., Bustamante M., Oikonomou F., Resconi E., 2021, arXiv e-prints, p. arXiv:2107.13534
  * The ATLAS Collaboration (2011) The ATLAS Collaboration The CMS Collaboration T. L. H. C. G., 2011, Technical report, Procedure for the LHC Higgs boson search combination in Summer 2011\. CERN, Geneva
  * The LHAASO Collaboration et al. (2021) The LHAASO Collaboration et al., 2021, Science, 373, 425
  * Tibet AS$\gamma$ Collaboration et al. (2021) Tibet AS$\gamma$ Collaboration et al., 2021, Nature Astronomy, 5, 460
  * Wakely & Horan (2008) Wakely S. P., Horan D., 2008, International Cosmic Ray Conference, 3, 1341
  * Wald (1943) Wald A., 1943, Transactions of the American Mathematical Society, 54, 426
  * Zaborov (2020) Zaborov D., 2020, arXiv e-prints, p. arXiv:2011.09209
  * Zhou et al. (2021) Zhou B., Kamionkowski M., Liang Y.-f., 2021, Phys. Rev. D, 103, 123018

##  Appendix A TeV gamma-ray Counterparts

In this appendix, we describe in some details about the TeV counterparts of
LHAASO sources. Table 2 gives the distances and the extensions of these TeV
counterparts.

LHAASO source | Extension | TeV Counterpart | Distance | $\rm exp(-\tau)$ | Angular Distance | Extension | Reference
---|---|---|---|---|---|---|---
| [deg] |  | [kpc] |  | [deg] | [deg] |
J0534+2202 | PS | Crab Nebula | 2.0 | 0.99 | 0.08 | 0.01 | [1]
J1825-1326 | 0.30 | HESS J1826-130 | 4.0 | 0.92 | 0.39 | 0.21 | [2]
|  | HESS J1825-137 | 3.9 | 0.92 | 0.39 | 0.55 (0.51) | [3]
|  | eHWC J1825-134 | — | — | 0.09 | 0.36 | [4]
|  | HAWC J1825-134 | — | — | 0.03 | PS | [5]
J1839-0545 | ES | HESS J1841-055 | — | — | 0.29 | 0.41 | [6]
J1843-0338 | ES | HESS J1843-033 | — | — | 0.22 | 0.24 | [6]
|  | HESS J1844-030 | — | — | 0.69 | PS | [6]
|  | 3HWC J1843-034 | — | — | 0.30 | PS | [7]
J1849-0003 | ES | HESS J1849-000 | 7 | 0.84 | 0.11 | 0.09 | [6]
J1908+0621 | 0.58 | HESS J1908+063 | — | — | 0.11 | 0.34 | [8]
|  | eHWC J1907+063 | — | — | 0.14 | 0.52 | [4]
J2018+3651 | ES | VER J2019+368 | — | — | 0.10 | 0.34 (0.14) | [9]
|  | eHWC J2019+368 | — | — | 0.17 | 0.20 | [4]
|  | TASG J2019+368 | — | — | 0.19 | 0.28 | [10]
J2032+4102 | ES | 3HWC J2031+415 | — | — | 0.47 | PS | [7]
|  | TeV J2032+4130 | 1.8 | 0.98 | 0.52 | 0.16 (0.07) | [11]
|  | ARGO J2031+4157 | 1.4 | 0.99 | 1.5 | 1.8 | [12]
|  | TASG J2032+414 | — | — | 0.41 | PS | [13]
|  | MGRO J2031+41 | — |  | 0.51 | 3.0 | [14]
J2226+6057 | 0.36 | G106.3+2.7 (VERITAS) | 0.8 | 0.99 | 0.14 | 0.27 (0.18) | [15]
|  | G106.3+2.7 (${\rm Tibet\,AS_{\gamma}}$) | 0.8 | 0.99 | 0.11 | 0.24 | [16]
J1929+1745 | ES | 3HWC J1928+178 | — | — | 0.16 | PS | [7]
J1956+2845 | ES | 3HWC J1954+286 | — | — | 0.33 | PS | [7]
J2108+5157 | PS | — | — | — | — | — | —
Table 2: The TeV gamma-ray counterparts associated with the LHAASO sources.
The second column is the intrinsic extension from LHAASO measurements. PS
represents point-like source whereas ES extended source. The prior intrinsic
extension for extended LHAASO source is $0.3^{\circ}$ (Cao et al., 2021a). The
seventh column is the intrinsic extension of TeV counterparts. The number
inside (outside) the brackets is the extension along the minor (major) axis.
The extension in boldface indicates the source extension adopted in Table 1.
The distance from the TeV counterpart to us is taken from TeVCat (Wakely &
Horan, 2008). The flux attenuation $\rm exp(-\tau)$ due to the ISRF and CMB is
given at $E_{\gamma}=100\,{\rm TeV}$. The extensions of TeV counterparts and
their angular distances to the LHAASO source refer to [1] H. E. S. S.
Collaboration (2020), [2] H. E. S. S. Collaboration et al. (2020), [3] H. E.
S. S. Collaboration et al. (2019), [4] Abeysekara et al. (2020), [5] Albert et
al. (2021), [6] H. E. S. S. Collaboration et al. (2018), [7] Albert et al.
(2020), [8] Aharonian et al. (2009), [9] Abeysekara et al. (2018), [10]
Amenomori et al. (2021), [11] Aliu et al. (2014), [12] Bartoli et al. (2014),
[13] Amenomori et al. (2021), [14] Abdo et al. (2007), [15] Acciari et al.
(2009) and [16] Tibet AS$\gamma$ Collaboration et al. (2021).

Crab Nebula. LHAASO-KM2A has detected an 1.1 PeV photon from the direction of
Crab Nebula and its energy spectrum has a possible hardening around PeV
energies which indicates a hadronic component (The LHAASO Collaboration et
al., 2021). So we consider the model by (Amato et al., 2003), in which protons
take $f_{p}=60\%$ of pulsar wind energy and produce gamma-rays and neutrinos
through p-p interactions. The target density $n_{t}$ is expressed as

$n_{t}=10\mu(M_{N\odot}/R_{\rm pc}^{3})\,{\rm cm}^{-3}$ (18)

with the parameter $\mu$ defined in the Equation 9 of Amato et al. (2003). The
p-p interactions are dominant in comparison with the photo-meson ($p\gamma$)
productions if $\mu$ is not very much smaller than unity. As for the target
density $\mu=1$, the gamma-ray flux will be higher than the observed flux if
$f_{p}>15\%$ and the peak energy will be much lower than PeV if the wind
Lorentz factor $\Gamma<10^{7}$. The density parameter $\mu$ equals to 5 if the
mass estimated by Fesen et al. (1997) uniformly distributes in the nebula
(Amato et al., 2003). The higher target density requires lower $f_{p}$ to
explain the observation.

LHAASO J1825-1326. HAWC has resolved the region around eHWC J1825-134 into
three sources (HAWC J1825-138, HAWC J1826-128 and HAWC J1825-134) and
discovered the point-like source HAWC J1825-134 whose energy spectrum extends
well beyond 200 TeV without a cutoff (Albert et al., 2021). As the angular
distance between HAWC J1825-134 and LHAASO J1825-1326 is only $0.03^{\circ}$
(see Table 2), we fit the energy spectrum of HAWC J1825-134 in combination of
the measurements by HAWC and LHAASO. The extrapolated flux ($>100\,{\rm TeV}$)
of HAWC J1825-138 and HAWC J1826-128 are removed from the measurements by
LHAASO.

LHAASO J2032+4102. The angular distances between LHAASO J2032+4102 and its TeV
counterparts are larger than the angular resolution of LHAASO-KM2A around
15-20 arcmin at 100 TeV (Cao et al., 2021a). Such angular distances indicate
that hundred TeV gamma-rays and TeV gamma-rays are probably generated in
different emission regions or astrophysical sources. Further studies on the
morphology and spectrum of this region are required to reveal the origin of
PeV gamma-rays from Cygnus region.

LHAASO J2108+5157. LHAASO J2108+5157 is the only LHAASO source without TeV
counterparts. As the constraint on TeV gamma-ray flux is not that strong, we
refer the hadronic model shown in Figure 5 of Cao et al. (2021b) and use the
ECPL2 in Equation 15 to model the profile of hadronic gamma-ray spectrum.

Rescaling due to extension. In Figure 3, the flux measurements of HESS
J1826-130, HESS J1841-055, HESS J1908+063, VER J2019+368, TeV J2032+4130 and
G106.3+2.7 (VERITAS) are rescaled according to their intrinsic extensions and
integration regions. See details in Appendix B of Huang & Li (2022).

Figure 6: The statistical significance (p-value) of neutrinos from LHAASO
sources as function of source extension $\sigma_{s}$. The red (blue) lines
correspond to the spectral models in red (blue) lines in Figure 3. These
results are from the ten-year IceCube data.

##  Appendix B Signal Significance for Different Source Extensions

In addition to the upper limits shown in Figure 3 and Figure 4, we give the
significance of signal neutrino events from the direction of LHAASO sources
here. A test statistic ($TS$) is built as

${TS}=2{\rm ln}\left[\frac{L\left(X\big{|}\hat{n}_{s},\sigma_{s},\Phi^{\rm
model}_{s}\right)}{L\left(X\big{|}n_{s}=0,\sigma_{s},\Phi^{\rm
model}_{s}\right)}\right],$ (19)

where $L$ is the likelihood function in Equation 1, $X$ represents the
observed data and $\hat{n}_{s}$ is the number of signal events ($n_{s}\geq 0$)
maximizing the likelihood with given source extension $\sigma_{s}$ and
spectral shape $\Phi^{\rm model}_{s}$ (e.g. $E^{-\gamma}$) of incident
astrophysical neutrinos. Assuming the validity of the Wald approximation
(Wald, 1943), $TS$ will follow the PDF as

$f({TS})=\frac{1}{2}\delta\left({TS}\right)+\frac{1}{2}\frac{1}{\sqrt{2\pi}}\frac{1}{\sqrt{TS}}e^{-{TS}/2}$
(20)

for pure background hypothesis (Cowan et al., 2011). So the statistical
significance of signal neutrino events can be described with the pre-trial
p-value as $1-F(TS)$, where $F(TS)$ is the CDF of $f(TS)$:

$F(TS)=\frac{1}{2}\left[1+\operatorname{erf}\left(\sqrt{TS/2}\right)\right].$
(21)

The p-values for different source extensions are shown in Figure 6.