Determining the Origin of Very-high-energy Gamma Rays from Galactic Sources by Future Neutrino Observations

Recently, the Large High Altitude Air Shower Observatory (LHAASO) identified 12 gamma-ray sources emitting gamma rays with energies above 100 TeV, making them potential PeV cosmic-ray accelerators (PeVatrons). Neutrino observations are crucial in determining whether the gamma-ray radiation process is of hadronic or leptonic origin. In this paper, we study three detected sources, LHAASO J1908+0621, LHAASO J2018+3651, and LHAASO J2032+4102, which are also the most promising Galactic high-energy neutrino candidate sources with the lowest pretrial p-value based on the stacking searches testing for excess neutrino emission by the IceCube Neutrino Observatory. We study the lepto-hadronic scenario for the observed multiband spectra of these LHAASO sources considering the possible counterpart source of the LHAASO sources. The very-high-energy gamma rays are entirely attributed to the hadronic contribution; therefore, the most optimistic neutrino flux can be derived. Then, we evaluate the statistical significance (p-value) as a function of the observation time of IceCube and the next-generation IceCube-Gen2 neutrino observatory, respectively. Our results tend to disfavor that all gamma rays above 100 GeV from LHAASO J1908+0621 are of purely hadronic origin based on current IceCube observations, but the purely hadronic origin of gamma rays above 100 TeV is still possible. By IceCube-Gen2, the origin of gamma rays above 100 TeV from LHAASO J1908+0621 can be further determined at a 5σ significance level within a running time of ∼3 yr. For LHAASO J2018+3651 and LHAASO J2032+4102, the required running time of IceCube-Gen2 is ∼10 yr (3σ) and ∼10 yr (5σ), respectively. Future observations by the next-generation neutrino telescope will be crucial to understanding the particle acceleration and radiation processes inside the sources.


INTRODUCTION
The origin of high-energy cosmic rays (CRs) has been a long-standing question in particle astrophysics.The cosmic ray spectrum is typically described by a power law with an index of ∼ 2.7 up to the so-called "knee" at around 3 PeV, beyond which the spectrum softens (Abbasi et al. 2018).CR composition changes from light (protons and helium) to heavier as a function of energy near the knee and is known to be dominated by protons (Gaisser et al. 2013).This suggests the existence of powerful astrophysical proton accelerators in Corresponding author: Kai Wang kaiwang@hust.edu.cnour galaxy, which can accelerate protons to energies up to a PeV, commonly referred to as "PeVatrons".The potential galactic PeVatrons can be identified by the very-high-energy (VHE, > 100 GeV) gamma-ray detection and have been explored by ground-based telescopes, such as H.E.S.S. (High Energy Stereoscopic System) (HESS Collaboration et al. 2016;Abdalla et al. 2018), MAGIC (Major Atmospheric Gamma Imaging Cherenkov) (Acciari et al. 2020), HAWC (High-Altitude Water Cherenkov) (Albert et al. 2020a), and LHAASO (Large High Altitude Air Shower Observatory) (Cao et al. 2021a(Cao et al. ,b, 2023)).Especially, LHAASO can capture gamma rays with energies from hundreds of GeV to exceeding PeV (LHAASO collaboration et al. 2010), the sources of which have strong possibilities of being galactic PeVatrons.However, the origin of VHE gammarays is still in debate.The VHE gamma rays can be produced through the decay of pions which are generated by hadronic processes between the accelerated high-energy cosmic rays and the surrounding medium.An alternative scenario is leptonic processes, such as inverse Compton scattering and bremsstrahlung of the accelerated high-energy electrons, which can also produce high-energy gamma rays.Therefore, confirming the gamma-ray origin is crucial to identifying the composition of accelerated particles and the radiation processes for PeVatrons.
Neutrinos are a significant probe in determining the origin of VHE gamma rays, which are concomitantly produced with gamma rays of hadronic origin.Therefore, the detection or non-detection of high-energy neutrinos can be a diagnosis of the hadronic or leptonic origin of VHE gamma rays.Very recently, neutrino emission from the galactic plane has been identified at the 4.5σ level of significance by IceCube Neutrino Observatory (Abbasi et al. 2023), which implies that galactic sources can generate high-energy neutrinos.In addition, for the LHAASO sources, Abbasi et al. (2023) conducted stacking searches testing for excess neutrino emission from 12 LHAASO sources which are identified with emissions above 100 TeV (Cao et al. 2021a) and thought to be PeVatron candidates.Although no significant neutrino emissions were found, three LHAASO sources, i.e., LHAASO J1908+0621, LHAASO J2018+3651, and LHAASO J2032+4102, present the lowest p-values, making them as the promising sources to identify the possible neutrino emission in the future and then judge the origin of VHE gamma rays.
Due to the complex spatial morphology of three LHAASO sources, the gamma-ray counterparts for these sources are still uncertain, which can be supernova remnants (SNRs), pulsar wind nebulae (PWNe), or young massive star clusters (YMCs).SNRs are widely thought of as the primary galactic cosmic-ray sources, and the particles can be accelerated by diffusive shock acceleration in their forward shocks generated by the interaction of supernova ejecta with the interstellar medium (ISM).The production of gamma rays and high-energy neutrinos from interactions of accelerated CR protons and nuclei with ambient medium (Gabici & Aharonian 2007).CR protons can also be accelerated and trapped in PWNe and YMCs and then produce gamma rays and neutrinos (Di Palma et al. 2017).
In this paper, we collect the multiband spectra observed from the direction of these sources and consider the possible counterpart sources of the LHAASO sources.With the proposed theoretical scenario, the multiband spectral modeling is implemented and VHE gamma rays are mainly attributed to the hadronic process.With the most optimistic neutrino production in the sources, we evaluate the statistical significance as a function of observation time for three sources, i.e., LHAASO J1908+0621, LHAASO J2018+3651, and LHAASO J2032+4102, using the IceCube and nextgeneration IceCube-Gen2 neutrino observatory respectively.
The remaining part of this paper is organized as follows.In section 2, we obtain the SED of sources through the lepto-hadronic scenario.In section 3, we calculate the corresponding neutrino SED from the hadronic interaction and compare the calculation with the sensitivity of the IceCube-Gen2 observatory.In section 4, we estimate the statistical significance of neutrino signals from LHAASO sources using both the IceCube and the proposed IceCube-Gen2.Finally, section 5 is the discussion and summary.
The origin of the gamma-ray emission from LHAASO J1908+0621 is still under debate due to its complex spatial morphology.In principle, PSR J1907+0631 can potentially power the entire TeV emissions, however, considering that it is not a gamma-ray pulsar, we neglect its contributions to gamma-ray emissions (Duvidovich et al. 2020;Lyne et al. 2017).Other possible scenarios can be concluded as follows.A leptonic component from the PWN powered by PSR J1907+0602 was initially proposed as the origin of VHE gamma rays (Abdo et al. 2010).The combination of leptonic and hadronic scenarios has been proposed as well for the origin of gamma-ray emission in this region (Duvidovich et al. 2020;Crestan et al. 2021;Li et al. 2021;Albert et al. 2022;De Sarkar & Gupta 2022).The hadronic component is usually related to SNR G40.5-0.5.As suggested in De Sarkar & Gupta (2022), we consider the pp interaction of escaped protons accelerated by the shock of SNR G40.5-0.5 with the materials of MCs in the MC region, and a leptonic component from SNR G40.5-0.5, which is located at a distance of 2.37 kpc, similar to Albert et al. (2021) and Cao et al. (2023).
The acceleration mechanisms of cosmic rays have been studied for a long time.One of the most popular acceleration mechanisms is the diffusive shock acceleration related to the supernova remnants (Fermi (1949); Drury (1983); Schure et al. (2012)).Charged particles would be accelerated by the shock waves produced by the supernova explosion.The shock wave from the supernova explosion expands and eventually reaches the MCs (Fujita et al. 2009).Subsequently, CRs escape the confinement region and begin to seep into the MCs when the escaping boundary contacts the surface of the MCs.Consequently, the hadronic component comprises the γ-ray produced from the interaction between escaped protons from SNR G40.5-0.5 and cold protons inside the associated MCs (Makino et al. 2019;De Sarkar & Gupta 2022).
We use the semi-analytical formulation by Kelner et al. (2006) to calculate the γ-ray spectra from interactions between injected CR protons and the MC materials that surround the SNR G40.5-0.5.CR protons with an exponential cutoff power law distribution are adopted, i.e., N p ∝ γ −αp p e −γp/γp,cut , where N p is the number of protons in a unit volume and in the energy interval (γ p , γ p + dγ p ), γ p is the proton Lorentz factor, α p is the spectral index and γ p,cut is the cutoff Lorentz factor.In this work, we invoke two scenarios to be responsible for the origin of VHE gamma-rays.The first scenario (Case 1) is that gamma rays above 100 GeV are totally attributed to the hadronic origin, and the second (Case 2) is that only gamma rays above 100 TeV are contributed by the hadronic process as a possible spectral hardening at higher energies may originate from the hadronic process (Albert et al. 2022) and Case 1 presents a slight spectral deviation from observations above 10 TeV (see next spectral fittings).
In addition to the contributions of hadrons discussed above, we also consider the contribution from the leptonic emission of relativistic electrons in the SNR+MC system.We have considered different leptonic radiation mechanisms, such as synchrotron, bremsstrahlung and IC (Blumenthal & Gould 1970;Baring et al. 1999).To calculate the IC contribution from MCs, we have considered the contribution from the interstellar radiation field (ISRF) model (Popescu et al. 2017) and Cosmic Microwave Background (CMB) (see the left panel of Fig. 1).The spectrum of the electron distribution is assumed to be a single power law with an exponential cutoff as for protons, i.e., N e ∝ γ −αe e e −γe/γe,cut , where N e is the number of electrons in a unit volume and in the energy interval (γ e , γ e + dγ e ), γ e is the electron Lorentz factor, α e is the index and γ e,cut is the electron cutoff Lorentz factor.
The multi-wavelength spectral energy distribution (SED) of the source LHAASO J1908+0621 is shown in Fig. 2 including two different scenarios.We implement a spectral fitting and the adopted model parameters are summarized in Table 1.It is worth noting that in the energy range of 1-10 TeV, we give priority to the HAWC data rather than H.E.S.S. data since HAWC is better at measuring the flux of extended sources.Imaging Atmospheric Cherenkov Telescopes (IACT, e.g., H.E.S.S.) will underestimate the flux of extended sources as they usually use blank sky regions near gamma-ray sources as background emissions.However, for a large extended source, there might be dim emission in the rim taken as background, causing a lower flux measurement.In terms of the additional component in the MeV-GeV range shown in Fig. 2, a consistent emission from the bremsstrahlung process is derived for Case 1, while for Case 2, the bremsstrahlung radiations move to higher energy band and therefore a direct fitting result from Li et al. ( 2021) is adopted to avoid involving the second population of electrons or protons.For two cases, the total energies of the injected protons required during spectral fittings are 1.1 × 10 47 or 5 × 10 46 erg, respectively (see Table 1).Both values are lower than the usual 1−10% of the kinetic energy released in SNRs (typically, E SN = 10 51 erg) (Aharonian et al. 2004).This could be attributed to the SNR age or the strength of the magnetic field in the SNR+MCs system.The adopted magnetic field strength and the MC number density are the same as in Albert et al. (2022).1.

LHAASO J2018+3651
MGRO J2019+37, the MILAGRO counterpart of LHAASO J2018+3651, is one of the brightest sources in the sky at TeV energies.This source is suspected to be associated with the GeV pulsar J2021+3651 (Abdo et al. 2009).The estimated distance of PSR J2021+3651 ranges from 2 to 12 kpc (Hou et al. 2014).Here we adopt a distance of 1.8 kpc as the same as in Fang et al. (2020) and the radius of its PWN is 24.6 pc.The PWN G75.2+0.1 of the pulsar PSR J2021+3651 is treated as a source of radiations of LHAASO J2018+3651.
In the PWN scenario, the leptonic origin for the multiwavelength emissions of LHAASO J2018+3651 can be naturally expected.For the leptonic processes including synchrotron radiations and IC scatterings, the same electron distribution as used for LHAASO J1908+0621 is adopted, i.e., a single power law with an exponential cutoff.For the IC scatterings, as shown in the middle panel of Fig. 1, the contributions from the ISRF model and CMB are taken into account as well.In fact, apart from gamma rays produced by leptons, there may also be hadronic components involved.A portion of the pulsar's spin-down power can be converted into a stream of nuclei and nuclei can be accelerated in pulsar magnetospheres (Hoshino et al. 1992;Arons & Tavani 1994;Gallant & Arons 1994).Accelerated nuclei can undergo photodisintegration when they collide with low-energy photons generated in the nonthermal radiation fields of the pulsar's outer magnetosphere.This process results in the release of energetic neutrons that subsequently decay either within or outside the Nebula.When the protons resulting from neutron decay collide with the matter within the nebula, they generate gamma-rays and neutrinos (Bednarek & Protheroe 1997;Liu & Wang 2021).Therefore, in addition to the IC process, the hadronic components can contribute γrays in the PWN scenario.An exponential cutoff power law distribution of protons is employed as the same as used for LHAASO J1908+0621.The multi-wavelength SED of LHAASO J2018+3651 and the spectral modeling are shown in Fig. 3(a) and the adopted parameters are summarized in Table 2.
Several studies have associated LHAASO J2018+3651, the counterpart of HAWC J2019+368 or MGRO J2019+37, with the PWN G75.2+0.1 powered by PSR J2021+3651 (Albert et al. 2021;Fang et al. 2020).They explain the multiband non-thermal emission via synchrotron radiation and inverse Compton scattering in the leptonic scenario (also see Woo et al. ( 2023)).Hou et al. (2014) respectively employed the separated leptonic model or hadronic model responsible for the VHE gamma rays.Here, we adopt a leptohadronic model to interpret multiband emission.The ambient proton density of the PWN is still unknown (Beacom & Kistler 2007) and we adopt 1 cm 3 , which is the same as the average density of the interstellar medium.The magnetic field is consistent with Albert et al. (2021).

LHAASO J2032+4102
The extended TeV gamma-ray source ARGO J2031+4157 (or MGRO J2031+41), the counterpart of LHAASO J2032+4102, is positionally coincident with the Cygnus Cocoon.No significant changes in morphology or spectrum have been observed for this extensive region (Ackermann et al. 2011).However, the energy spectrum from 1 GeV to 10 TeV suggests that the Cygnus Cocoon might either be an unidentified SNR or that the particle acceleration within a superbubble is similar to that within an SNR (Bartoli et al. 2014).
Although PSR J2032+4127 may contribute to the multiband emission of the source LHAASO J2032+4102, we choose an unknown SNR+MC system as the potential single source to power multiband observations from the LHAASO J2032+4102 region since the observed extended X-ray emission may be the counterpart of the TeV emission (Bartoli et al. 2014).The magnetic field B ≈ 3 µG is involved in order to be consistent with previous works (Horns et al. 2007).For the IC scatterings, the seed photons from the ISRF model and CMB are considered (see right panel of Fig. 1).Additionally, we attribute the observed γ-rays to the decay of π 0 mesons generated through inelastic collisions between accelerated protons and target gas in an unidentified SNR+MC system as for LHAASO 1908+0621.Both the electron and proton distributions are assumed to be a power-law The radio flux of a possible nonthermal extended radio source at λ = 20 cm (Paredes et al. 2006) is assumed to be an upper limit of the actual radio emission.The upper limits between 0.5-5 keV (CX, Chandra) and 20-40 keV (ISGRI, INTEGRAL) are taken from (Butt et al. 2006).The model parameters of the two panels are summarized in Table 2.
function with an exponential cutoff as above.The multiwavelength SED and spectral modelings of LHAASO J2032+4102 are shown in Figure 3(b) and the adopted parameters are concluded in Table 2.The MC number density is adopted as a value of 30 cm −3 .The total energy (≃ W p ) is 3 × 10 50 erg, which can be reasonably provided by one supernova, which typically releases ∼ 10 51 erg, and about 10% of which can be transferred to the accelerated particles.

NEUTRINO FLUX
Neutrinos are produced alongside γ-rays in hadronic pp interactions.Therefore, if there are γ-ray sources powered by hadronic interactions, it is also expected that neutrinos will be emitted from the same source region.To calculate the flux of muon neutrinos reaching Earth from the three sources, we use the semi-analytical formulation developed in Kelner et al. (2006).The derived muon neutrino flux is shown in Fig. 4 after considering neutrino oscillations.Neutrino oscillations induce an equal flux among different neutrino flavors when they reach Earth, namely, the (anti-)muon neutrino flux reaching Earth is one-third of all-flavor neutrino flux, i.e., (Φ νµ + Φ νµ + Φ νe + Φ νe )/3 with produced neutrino fluxes (Φ νµ , Φ νµ , Φ νe , Φ νe ) by pp interactions at the source.
As shown in Fig. 4, the model predicts a neutrino flux that exceeds the sensitivity limit of next-generation IceCube-Gen2.This indicates that if the total observed TeV/PeV gamma-ray emission from these LHAASO sources originates from the hadronic processes, the accompanying neutrino flux can be identified by IceCube-Gen2.As no accompanying neutrinos will be expected if VHE gamma rays are produced by leptonic processes, the detection or non-detection by IceCube-Gen2 in the future can provide further insights and help to determine the origin of VHE gamma rays from these LHAASO sources.

FUTURE PROSPECTS
In this section, we evaluate the detection of neutrinos from three sources by IceCube and IceCube-Gen2.The number of signal events from the point-like source at the declination δ s is expressed as  where A eff is the effective area for through-going track events.The background events from the three sources are mainly induced by atmospheric neutrinos.Thus, the number of background events is expressed as where I ν,atm is the atmospheric muon neutrino flux calculated by MCEq (Fedynitch et al. 2015), and I ν,astro is the diffuse astrophysical muon neutrino flux (Abbasi et al. 2022).As for IceCube, the effective area and the smearing matrix are released with the ten-year (2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018) muon track data (IceCube Collaboration 2021).
The smearing matrix gives the reconstructed energy distribution for neutrinos with different energies and incident directions.As for IceCube-Gen2, the effective area is assumed to be 7.5 times that of IceCube (Schumacher et al. 2022), and the smearing matrix is assumed to be the same as that of IceCube.Low-energy cut is applied to both IceCube and IceCube-Gen2 observations, aiming for the neutrino-induced muons with reconstructed energy E rec < 50 TeV (Omeliukh et al. 2021).
We estimate the statistical significance of observation with a p-value analytically expressed as (The AT-LAS Collaboration 2011; Halzen et al. 2017) where The results for the statistical significance are reported in Fig. 5.With the observation by IceCube, LHAASO J1908+0621 is expected to reach 3σ in ∼2020 (∼2 years after 2018) for Case 1.In the latest observation by Ice-Cube (2008-2020) (Abbasi et al. 2023), the pre-trial pvalue for LHAASO J1908+0621 is only 0.046 under the assumption of a power-law spectrum.Therefore, our results tend to disfavor Case 1 for LHAASO J1908+0621, which corresponds to that gamma rays above 100 GeV are of purely hadronic origin.LHAASO J1908+0621 will be detected at 5σ level in less than 10 months for Case 1, while for Case 2 will be detected at 5σ level in less than 3 years with IceCube-Gen2 detec-tor, considering the relevant parameters as reported in Table 1.Therefore, future observations by the next generation neutrino telescope can further help distinguish the origin of VHE gamma-rays from LHAASO J1908+0621.The source LHAASO J2018+3651 and LHAASO J2032+4102 could be detected at 3σ level in ∼10 years and 5σ level in ∼4 years by IceCube-Gen2, respectively, considering the relevant parameters as reported in Table 2.However, no significant statistical significance will be expected by IceCube for LHAASO J2018+3651 and LHAASO J2032+4102 in the next 10 years.
5. DISCUSSION AND SUMMARY Galactic high-energy neutrinos have been long expected.Galactic sources detected by LHAASO are thought of as the potential PeVatrons since gamma rays with energies larger than 100 TeV have been found from them.However, it is still under debate whether VHE gamma rays are produced by the leptonic process or the hadronic process.Neutrino observations can be an important probe to distinguish the origin of VHE gamma rays since the hadronic process (pp collisions or photomeson production process) will produce the accompanying high-energy neutrinos.In this paper, we investigate multiband spectra of three LHAASO sources, i.e., LHAASO J1908+0621, LHAASO J2018+3651, and LHAASO J2032+4102, which are the most promising galactic neutrino sources.We propose reasonable leptohadronic scenarios to implement the multiband spectral modeling.
Assuming gamma rays are entirely hadronic, we calculate the most optimistic flux of muon neutrinos generated from the hadronic process.Furthermore, we estimate the statistical significance (p-value) as a function of time for three sources using both the IceCube neutrino observatory and the proposed second-generation IceCube-Gen2.Our results tend to disfavor that all gamma rays above 100 GeV (Case 1) from LHAASO J1908+0621 are of purely hadronic origin based on current IceCube observations.However, the purely hadronic origin of gamma rays above 100 TeV (Case 2) from LHAASO J1908+0621 is still possible.No significant statistical significance will be expected by IceCube for LHAASO J2018+3651 and LHAASO J2032+4102 in the next 10 years.Besides, for IceCube-Gen2, LHAASO J1908+0621 can be detected at a 5σ level within 10 months for Case 1, and within 3 years for Case 2. Similarly, high-energy neutrinos from LHAASO J2018+3651 and LHAASO J2032+4102 can be respectively detected by IceCube-Gen2 at a 3σ level in ∼10 years and a 5σ level in ∼4 years if the VHE gamma rays are entirely hadronic.Future observations by IceCube-Gen2 or other more advanced next-generation neutrino telescopes at the positions of three sources will be important to untangle the exact nature of these enigmatic sources.

Figure 1 .Figure 2 .
Figure 1.Radiation field of each source.The red dashed line is the cosmic microwave background (CMB) energy density.The grey curve is the sum of CMB and interstellar radiation field (ISRF) at diverse galactocentric positions.

Figure 4 .
Figure 4.The estimated total muon neutrino flux reaching the Earth of LHAASO J1908+0621, J2018+3651, and J2032+4102, respectively.The red, brown, purple, and pink lines represent the total muon neutrino fluxes above 1 TeV produced by pp interactions.The light blue solid line indicates the sensitivity of IceCube-Gen2 for a point source at the celestial equator with an average significance of 5σ after 10 years of observations (Grant et al. 2019).

Figure 5 .
Figure 5.The statistical significance (p-value) as a function of observation time for three sources.The blue lines show the expected p-values obtained using IceCube in addition to the ten years from 2008 to 2018.The orange lines show the expected p-values obtained using IceCube-Gen2.

Y
b is the expected number of background events, and N D is the median of Poisson-distributed events containing both signal and background.The event numbers are counted within the solid angle Ω = 1.6σ,where σ is the angular resolution of the detector.This solid angle contains 72% of the signal from a point-like source(Alexandreas et al. 1993).The angular resolution σ of Ice-Cube is assumed to be the median angular uncertainty of the muon-track events (E rec < 50 TeV) observed within the declination band δ s ± 1 • .Omeliukh et al. (2021) reports the angular resolutions for low-energy samples simulated by IceCube-Gen2 (see Table2therein), but only three zenith bins are offered.Thus, we extrapolate the angular resolution to the source declination δ s as σ(δ s ) = σ 0 (δ s ) + 0.05 • , where σ 0 represents the angular resolution of IceCube-Gen2 for 10 TeV muons (see Figure24inAartsen et al. (2021)).