Origin of Galactic Sub-PeV Diffuse Gamma-Ray Emission: Constraints from High-energy Neutrino Observations

The gamma-ray sky is dominated by the diffuse emission from the Galactic plane. The survey by gamma-ray satellites such as SAS-2 (Fichtel et al. 1975), EGRET on the Compton Gamma-Ray Observatory (Hunter et al. 1997), and Fermi Large Area Telescope (Fermi-LAT; Ackermann et al. 2012; Neronov & Semikoz 2020) have measured diffuse gamma-ray emission from the Galactic plane from several tens of MeV up to TeV energies. Ground-based gamma-ray instruments can also measure Galactic diffuse gamma rays from a fraction of the Galactic plane due to the limited observable sky. Milagro (Atkins et al. 2005) and ARGO-YBJ (Bartoli et al. 2015) have extended the Galactic diffuse gamma-ray spectrum up to several TeV. The leading radiation mechanism for the diffuse Galactic gamma-ray emission (DGE) is believed to be the decay of neutral pions generated from hadronic interactions of cosmic-ray (CR) protons with the interstellar medium (e.g., Dermer 1986; Mori 1997; Strong et al. 2010), whereas energetic CR electron/positron pairs escaping from pulsar wind nebulae around middle-aged pulsars may also have an important contribution to the DGE at the TeV band (Linden & Buckman 2018). Therefore, the DGE can serve as a probe of the CR distribution in the Galactic disk, and be used to study the CR propagation and their origin. Very recently, Amenomori et al. (2021) reported detection of 0.1–1 PeV DGE by the Tibet AS+MD array. All gamma rays above 398 TeV are observed apart from 05 of any known TeV sources and hence unlikely originate from leptonic sources since high-energy electrons cannot propagate far from sources before being cooled. The detection of diffuse sub-PeV gamma-ray emission of the hadronic origin provides a good opportunity to study the origin of PeV CRs.

Since high-energy neutrinos always accompany with the production of pionic gamma-rays, a diffuse high-energy neutrino background is expected from the Galactic plane. Recently, the IceCube neutrino telescope has obtained an upper limit on the neutrino flux of the Galactic plane in the energy range of 1–500 TeV (Aartsen et al. 2017). This motivates us to compare the diffuse gamma-ray emission measured by the Tibet AS+MD array and the neutrino flux constraint from IceCube. As will be shown in this Letter, we find that if we require the model neutrino intensity to be consistent with the 90% C.L. upper limit of Galaxy, the model gamma-ray intensity is lower than the measured flux by the Tibet AS+MD array especially at the highest-energy bin in 398–1000 TeV. It may indicate that some additional sources contribute to the diffuse sub-PeV gamma-rays measured by the Tibet AS+MD array. For example, Ahlers & Murase (2014) suggested that remnants of historical yet unresolved hypernovae, which are believed to be capable of energizing protons up to 1 EeV (Wang et al. 2007; Budnik et al. 2008), may contribute to sub-PeV DGE while not overproduce Galactic neutrinos.

Interestingly, as pointed out by Amenomori et al. (2021), 4 out of 10 measured events above 398 TeV are detected within 4° from the center of the so-called Cygnus Cocoon. Therefore, a nonnegligible fraction of the diffuse gamma-ray flux measured by the Tibet AS+MD array may originate from the Cygnus region. Indeed, various instruments have detected an excess of gamma-ray emission from GeV up to 100 TeV from Cygnus Cocoon (Ackermann et al. 2011; Bartoli et al. 2012; Abeysekara et al. 2018, 2021), which is probably a superbubble related to the massive star cluster Cygnus OB2 and is also in spatial coincidence with a pulsar wind nebula, PWN TeV J2032+4130 and a supernova remnant, γ Cygni. All of these objects are possible sources of energetic protons and electrons so that they are potential sources of these gamma rays. We speculate that a significant fraction of sub-PeV gamma rays detected by the Tibet AS+MD array may originate from the Cygnus region and/or some other extended sources in the Galactic disk, rather than the true diffuse gamma rays. We will study whether this scenario can alleviate the tension.

To calculate the diffuse gamma-ray and neutrino emission produced by the proton–proton (pp) collisions between CRs and the interstellar medium, we need to know their spatial distributions, i.e., n_CR(E_p , r, z) and n_ISM(r, z) with E_p being the CR energy, R the radius from the Galactic center projected in the Galactic plane, and z the height from the Galactic plane. We employ here the model developed by Lipari & Vernetto (2018), which can successfully reproduce both the flux (with an error of order 10%–20% depending on the direction) and the main features of the angular distribution of the diffuse gamma-ray flux measured by Fermi-LAT . The model parameterizes n_CR(E_p , R, z) and n_ISM(R, z) as well as the CR slope as a function of R and z analytically so it can save a lot of computation time for modeling the CR transport and spatial distribution. Following the "factorized model" in that paper, we can then obtain the diffuse gamma ray or neutrino intensity (in units of GeVcm⁻² s⁻¹ sr⁻¹) in a certain direction with Galactic coordinates (l, b) by performing the line-of-sight integration

$$\begin{eqnarray}&&\begin{array}{l}{I}_{\gamma /\nu }({E}_{\gamma /\nu },l,b)\\ \,=\,\displaystyle \int \displaystyle \int {{ \mathcal F }}_{\gamma /\nu }\left\{{E}_{\gamma /\nu },{n}_{\mathrm{CR}}({E}_{p},r,l,b),{n}_{\mathrm{ISM}}(r,l,b)\right\}\\ \quad \times \exp \left[-\tau ({E}_{\gamma },r)\right]{{dE}}_{p}{dr},\end{array}\end{eqnarray} \tag{ 1 }$$

where ${{ \mathcal F }}_{\gamma /\nu }$ denotes the operator calculating the gamma-ray or neutrino emissivity of pp collisions following the semianalytical method developed by Kelner et al. (2006), r is the distance of a certain place from Earth in the direction (l, b), and the corresponding R and z can be found by $R={\left({R}_{E}^{2}+{\left(r\cos b\right)}^{2}-2{R}_{E}r\cos b\cos l\right)}^{1/2}$ and $z=r\sin b$, where R_E = 8.5 kpc is the distance of Earth from the Galactic center. τ(E_γ , r) is the opacity of the pair production of high-energy gamma-ray photons on the cosmic microwave background and the interstellar radiation field. The latter is based on the model by Popescu et al. (2017). Note that this term will not appear in the equation if we calculate the neutrino intensity. To compare with measured DGE spectra from a certain region of the Galactic plane, we need to average the intensity over the region of interest, i.e.,

$$\begin{eqnarray}&&{\bar{I}}_{\gamma }={{\rm{\Omega }}}_{\gamma }^{-1}\int \int {I}_{\gamma }\sin (90^\circ -b){dbdl}.\end{eqnarray} \tag{ 2 }$$

For example, the DGE spectra shown in Amenomori et al. (2021) are extracted from 25° < l < 100° ∣b∣ < 5° (corresponding to Ω_γ = 0.228 sr), and 50° < l < 200° ∣b∣ < 5° (corresponding to Ω_γ = 0.456 sr), respectively. The neutrino analysis by IceCube uses the data including the entire Galactic plane, but is primarily sensitive to the northern hemisphere. Therefore, to compare with the neutrino upper limit, we only average the intensity over the Galactic plane in the northern hemisphere (or with decl. δ ≥ 0), i.e.,

$$\begin{eqnarray}&&{\bar{I}}_{\nu }={{\rm{\Omega }}}_{\nu }^{-1}\int \int {I}_{\nu }\theta \left(\delta (l,b)\right)\sin (\pi /2-b){dbdl},\end{eqnarray} \tag{ 3 }$$

where θ is the Heaviside function and δ can be found by $\sin \delta =\sin (27\buildrel{\circ}\over{.} 13)\sin b$ + $\cos (27\buildrel{\circ}\over{.} 13)\cos b\cos (122\buildrel{\circ}\over{.} 93-l)$, where 2713 is the decl. of the north Galactic pole and 12293 is the longitude of the northern equatorial pole. The corresponding solid angle of the region of interest for the neutrino emission is Ω_ν = 0.553 sr.

We compare the gamma-ray intensity and neutrino intensity predicted by the model with the observations in Figure 1. For the predicted gamma-ray intensity, our result is generally the same as that shown in Figure 20 of Lipari & Vernetto (2018). The model intensity is slightly lower than the measured intensity by the Tibet AS+MD array, in particular at the highest-energy bin, for both 25° < l < 100° and 50° < l < 200°. Such a difference is not significant considering the systematic error, as pointed out by Amenomori et al. (2021). However, we also see that the accompanying neutrino intensity exceeds the 90% confidence level (C.L.) upper limit of IceCube ⁴ below 100 TeV by about a factor of 2. In other words, there might be a tension between the diffuse sub-PeV gamma-ray data and the neutrino data, and if we want to reconcile the model neutrino intensity with the IceCube's upper limit, we need to multiply a factor of 0.5 to the neutrino intensity. In the meantime, it would reduce the predicted gamma-ray intensity by the same factor and make it insufficient to explain the DGE data, especially for the highest-energy bin. Note that there are some uncertainties of the model, such as the metallicity of both CRs and the interstellar medium that would influence the predicted gamma-ray or neutrino intensity up to a factor of 2, as well as the spatial variation of CR density and spectra in the Galactic plane, etc. Therefore, it is acceptable to have the flexibility of a factor of a few for the model intensity.

Zoom In Zoom Out Reset image size

In fact, the model uncertainty influencing the amplitude of the gamma-ray intensity would not affect our result significantly because we rescale the gamma-ray intensity according to the coproduced neutrino intensity and the IceCube upper limit. On the other hand, the spectral variation of CRs across the Galaxy is important. The factorized model of Lipari & Vernetto (2018) employed here considers the hardening of the CR spectrum toward the inner Galaxy, from a slope of 2.8 in the periphery of the Galaxy to 2.4 at the Galactic center. Lipari & Vernetto (2018) also constructed a "no-factorized model" with a uniform CR spectral slope, identical to the locally measured one, throughout the Galaxy. If we adopt the no-factorized model, the resulting diffuse gamma-ray spectrum would be softer than the present one for 25° < l < 100°. Consequently, it will be more insufficient to explain the measured intensity of the highest-energy bin solely with the diffusive CRs and strengthen our motivation to investigate additional contributions from extended sources. In contrast, if the CR spectral slope turns out to be harder than that in the factorized model in the inner Galaxy, the tension between the sub-PeV DGE and the Galactic neutrino flux upper limit can be alleviated.

It is reported that the total event number in the 398 TeV–1000 TeV bin observed by the Tibet AS+MD array is 10 in each of the regions 25° < l < 100° and 50° < l < 200°. Among both of them, 4 out of 10 photons originate from the region within 4° of the center of Cygnus Cocoon at l ≈ 80° (Amenomori et al. 2021). The Cygnus Cocoon region contains some potential sources for PeV CRs such as the massive star cluster Cygnus OB2 and the supernova remnant Gamma Cygni, so the energetic, fresh CRs may permeate throughout the entire region with a higher density than the average one in the interstellar medium. Although the analysis of the diffuse emission has excluded the region within 05 of the known TeV sources (Amenomori et al. 2021), the radial profile of the TeV gamma ray of Cygnus Cocoon is quite extended, which can be described by a Gaussian profile with a width of 21 (Abeysekara et al. 2021). Therefore, most of the emission of Cygnus Cocoon probably has been counted in the diffuse emission. As a result, we should also consider an additional component from Cygnus Cocoon when modeling the intensity in both regions of 25° < l < 100° and 50° < l < 200°.

The TeV spectrum of the Cygnus Cocoon region has been measured by the High-Altitude Water Cherenkov (HAWC) observatory up to 100 TeV, which can be described by a power-law function with an index of 2.6 (Abeysekara et al. 2021). If part of the DGE between 398 TeV–1000 TeV measured by the Tibet AS+MD array can be attributed to Cygnus Cocoon, it implies that its spectrum should continue somehow up to 1 PeV and indicate the existence of an extreme accelerator of CRs in that region processing protons up to at least 10 PeV. To evaluate the sub-PeV flux of Cygnus Cocoon, however, we need to know the relative exposure time of the instrument on the Cygnus region and on the other region of the Galactic plane, which is not given. If we simply assume a uniform exposure of the Tibet AS+MD array over the observed Galactic plane, the fact that 4 out of 10 total events in both 25° < l < 100° and 50° < l < 200° from Cygnus Cocoon would mean that a fraction f_cyg ≃ 40% of the total DGE flux (${{\rm{\Omega }}}_{\gamma }{\bar{I}}_{\gamma }$) at this energy originates from the source. Comparing it with the HAWC data, we see a flattening of the spectrum above 30–100 TeV (see the filled red and blue squares in Figure 2) and this would imply a different origin of the emission below 30 TeV. We note that this is most likely an overestimation of the sub-PeV flux since the exposure for Cygnus Cocoon would be longer than the average value for the Galactic plane, ⁵ but we take it as a reference case and refer to it as CASE I. On the other hand, if we assume the flux of Cygnus Cocoon at 398 TeV–1000 TeV is consistent with the extrapolation of HAWC's spectrum, Cygnus Cocoon would account for only about f_cyg ≃ 5% of the DGE intensity at the highest-energy bin measured by the Tibet AS+MD array (see the open red and blue squares in Figure 2). We refer to this case as CASE II.

Zoom In Zoom Out Reset image size

3.1. The Case of Hadronic Origin

The HAWC observation on Cygnus Cocoon favors a hadronic origin of the TeV gamma-ray emission (Abeysekara et al. 2021). For a phenomenological modeling of the gamma-ray emission of Cygnus Cocoon, we consider a total proton energy of W_p contained in the Cygnus region with a spectrum following ${N}_{p,\mathrm{cyg}}\equiv {{dN}}_{p}/{{dE}}_{p}={N}_{0}{E}^{-2}\exp (-{E}_{p}/{E}_{c})$ to explain the emission above 30 TeV of Cygnus Cocoon in CASE I, and we assume the proton spectrum to be ${N}_{p,\mathrm{cyg}}\equiv {{dN}}_{p}/{{dE}}_{p}\,={N}_{0}{E}_{p}^{-2}{\left(1+{E}_{p}/{E}_{b}\right)}^{-1}$ to account for also the spectrum below 1 TeV in CASE II, where N₀ can be determined by $\int {E}_{p}({{dN}}_{p}/{{dE}}_{p}){{dE}}_{p}={W}_{p}$.

We then can obtain the gamma-ray flux related to the sub-PeV DGE from Cygnus Cocoon by

$$\begin{eqnarray}\begin{array}{rcl}{F}_{\mathrm{cyg}},\gamma ({E}_{\gamma }) & = & \displaystyle \frac{1}{4\pi {r}_{\mathrm{cyg}}^{2}}\displaystyle \int {{ \mathcal F }}_{\gamma }\left\{{E}_{\gamma },{N}_{p,\mathrm{cyg}}({E}_{p}),{n}_{\mathrm{cyg}}\right\}\\ & & \times \exp \left[-\tau ({E}_{\gamma },{r}_{\mathrm{cyg}})\right]{{dE}}_{p},\end{array}\end{eqnarray} \tag{ 4 }$$

given the average gas density of n_cyg = 20 cm⁻³ in Cygnus Cocoon (Aharonian et al. 2019) and its distance r_cyg = 1.4 kpc. We find the derived gamma-ray flux can match both the HAWC data and the postulated sub-PeV flux, with E_b = 20 TeV and W_p = 2.5 × 10⁴⁹ erg in CASE I, and with E_c = 30 PeV and W_p = 1.5 × 10⁴⁸ erg in CASE II as shown in Figure 2. In this case, high-energy neutrinos are naturally expected from Cygnus Cocoon, as also predicted in some previous literature (e.g., Beacom & Kistler 2007; Evoli et al. 2007; Bi et al. 2009; Fox et al. 2013; Tchernin et al. 2013; Yoast-Hull et al. 2017). We calculated the coproduced neutrino flux and find it consistent with the sum of the upper limits for 2HWC J2031+415 and Gamma Cygni (Aartsen et al. 2020), which are two potential CR sources related to Cygnus Cocoon.

In addition to Cygnus Cocoon, there may be other sources injecting CRs with energy beyond 1 PeV into our Galaxy, although the others may be less powerful than Cygnus Cocoon. Similar to the case of Cygnus Cocoon, these super-PeV CRs diffuse quite fast, distributing over an extended region around their source, and produce gamma rays with a harder spectrum than that from the bulk of the interstellar medium. Without known the properties of these extended sources, we use Cygnus Cocoon's emission as a proxy of all these extended sources and parameterize the contribution of other sources by a coefficient ξ, i.e.,

$$\begin{eqnarray}&&{I}_{\mathrm{Ext}.\mathrm{Src},\gamma }({E}_{\gamma })={{\rm{\Omega }}}_{\gamma }^{-1}(1+\xi ){F}_{\mathrm{cyg},\gamma }({E}_{\gamma }).\end{eqnarray} \tag{ 5 }$$

The result is shown in Figure 3. In CASE I, it needs ξ = 0.5 to make the sum of the contribution by diffuse CRs and extended sources fit the diffuse gamma-ray spectra measured by the Tibet AS+MD array. The coproduced neutrino flux is consistent with the neutrino upper limit. Note that neutrinos now come from both discrete sources and diffusive CRs in the Galactic disk, so we compare it with the sum of the neutrino upper limits for diffuse Galactic emission as shown in Figure 1 and the sources. For the latter, we employ the results in Aartsen et al. (2017) for five catalogs of potential Galactic high-energy CR sources and sum up the upper limits of all the catalogs. In CASE II, since the contribution of Cygnus Cocoon is low, it needs ξ = 5 to fit the diffuse gamma-ray spectra. However, due to the soft proton spectrum in this case, the coproduced neutrino flux exceeds the upper limit. Therefore, the value of f_cyg, the fraction of the sub-PeV DGE that can be attributed to Cygnus region, is probably larger than 5% (but lower than 40% as CASE I is most likely an overestimation). With a larger f_cyg, a harder proton spectrum will be inferred and subsequently the neutrino flux at low energies can be reduced. A larger f_cyg also implies that a hardening in the gamma-ray spectrum of Cygnus Cocoon would appear beyond 30–100 TeV, and this might suggest a different origin for the TeV emission below 30 TeV. Note that the spectra of other sources are not necessarily the same as that of Cygnus Cocoon. Thus, given a soft spectrum of Cygnus Cocoon with a low flux at sub-PeV such as in CASE II, the DGE at sub-PeV energy could be still explained if other sources have harder spectra.

Zoom In Zoom Out Reset image size

3.2. The Case of Leptonic Origin

Although HAWC's observation on Cygnus Cocoon favors a hadronic origin, it does not exclude a leptonic origin completely for TeV gamma-ray emission above 1 TeV. Indeed, if there exists a powerful CR proton accelerator in the center of Cygnus Cocoon, it should be able to produce high-energy electrons too. The pulsar wind nebulae powered by PSR J2032+4127 and/or other unresolved pulsars may also inject high-energy electrons into ambient medium. Therefore, we also discuss a possible source contribution of leptonic origin to the sub-PeV DGE and we refer to this case as CASE III.

The sub-PeV gamma rays in the case of the leptonic origin can arise from the inverse Compton scattering (IC) of PeV electrons off the CMB radiation, whereas the IC process off radiation field of higher temperature such as the dust infrared radiation and the stellar optical–UV radiation are severely suppressed by the Klein–Nishina (KN) effect. In fact, even for CMB radiation with a typical energy of 6 × 10⁻⁴ eV, their scatterings with PeV electrons have also entered the KN regime, and thus we may expect the upscattered photon energy to be close to the electron energy, i.e., E_γ ≃ E_e . Electrons lose their energies as they radiate. The dominant cooling channel of electrons in the considered scenario is the synchrotron radiation with a cooling timescale of ${t}_{\mathrm{syn}}\simeq 500{\left({E}_{e}/1\mathrm{PeV}\right)}^{-1}{\left(B/5\mu {\rm{G}}\right)}^{-2}\,\mathrm{yr}$, as long as the magnetic field in Cygnus Cocoon is not much weaker than the typical interstellar magnetic field 3–5 μG (compare the synchrotron cooling timescale t_syn and the IC cooling timescale t_IC in Figure 5). Since the cooling timescale is much shorter than the dynamical timescale of the system, the IC spectrum can be estimated by ${E}_{\gamma }^{2}{N}_{\gamma }={E}_{e}^{2}{Q}_{e,\mathrm{inj}}{t}_{\mathrm{syn}}/{t}_{\mathrm{IC}}\propto {E}_{e}^{1-p}\propto {E}_{\gamma }^{1-p}$. Here, ${Q}_{e,\mathrm{inj}}\propto {E}_{e}^{-p}$ is the injection spectrum, ${t}_{\mathrm{syn}}\propto {E}_{e}^{-1}$ and ${t}_{\mathrm{IC}}\propto {E}_{e}^{0}$ approximately in the range of E_e = 0.1–1 PeV. Thus, we obtain ${N}_{\gamma }\propto {E}_{\gamma }^{-p-1}$. To explain the DGE measured by the Tibet AS+MD array, the spectrum should be harder than ${E}_{\gamma }^{-2.7}$, because otherwise the source contribution would overshoot the measured DGE intensity at 0.1–0.4 PeV while it fits the 0.4–1 PeV intensity. This leads to a requirement of a hard electron injection spectrum with p < 1.7.

On the other hand, a major difference of the leptonic case from the hadronic case is that electrons cool efficiently as they propagate while protons almost do not cool. Because the analysis of the Tibet AS+MD array masked the region within 05 of each known TeV source, those PeV electrons should at least diffuse to the region beyond 05 of the source, which corresponds to a radius of 12 pc at the nominal distance of 1.4 kpc for Cygnus Cocoon. It then puts a constraint on the diffusion coefficient by $\sqrt{4{{Dt}}_{\mathrm{syn}}}\gg 12\,\mathrm{pc}$, or $D(1\,\mathrm{PeV})\gg 2\,\times {10}^{28}{\left(B/5\mu {\rm{G}}\right)}^{2}\ {\mathrm{cm}}^{2}\ {{\rm{s}}}^{-1}$. This translates to $\beta \gg 0.015{\left(B/5\mu {\rm{G}}\right)}^{2}$, if we follow Abeysekara et al. (2021), which considers an average diffusion coefficient in the Galaxy to be ${D}_{\mathrm{ISM}}(E)\,=3\times {10}^{28}{\left(E/10\mathrm{GeV}\right)}^{1/3}{\mathrm{cm}}^{2}\,{{\rm{s}}}^{-1}$ and denote the diffusion coefficient in the Cygnus region by D_cyg(E) = β D_ISM(E). For a stronger magnetic field, a faster diffusion or a larger β will be needed. We show an example of the expected DGE and SED of Cygnus Cocoon under CASE III in Figure 4, where we employ p = 1.5, β = 0.1, and an electron luminosity of L_e,inj = 10³⁵ergs⁻¹. Note that high-energy neutrinos are not expected from the sources in this case. We see that the DGE measured by the Tibet AS+MD array is marginally explained although the model overpredicts the 200 TeV DGE intensity in 25° < l < 100°.

Zoom In Zoom Out Reset image size

To summarize, we calculated the diffuse gamma-ray intensity and neutrino intensity of the Galactic plane, which arise from hadronic interactions between CRs and the interstellar medium. We found that if we require the model neutrino intensity to be consistent with the 90% C.L. upper limit of Galaxy, the model gamma-ray intensity is lower than the measured flux by the Tibet AS+MD array, especially at the highest-energy bin in 398–1000 TeV. We speculated that an additional contribution by CRs with a relatively hard spectrum from the Cygnus region and/or some other extended sources of CRs beyond PeV is needed. Provided that these sources can accelerate electrons to energy beyond 1 PeV, the sub-PeV diffuse gamma-ray emissions may also be contributed by the inverse Compton scatterings of PeV electrons off CMB, as long as the injected electron spectra are harder than E^−1.7 and the diffusion coefficient in the ambient source region is much larger than $0.015{\left(B/5\mu {\rm{G}}\right)}^{2}$ times the average diffusion coefficient of the interstellar medium. In such a scenario, Cygnus Cocoon would harbor one (or more) extreme particle accelerators of either protons or electrons beyond PeV, and its gamma-ray spectrum might show a hardening above 30–100 TeV. If the sub-PeV gamma-ray flux from the Cygnus region turns out be low, other extreme particle accelerators in our Galaxy would make a more important contribution to the diffuse sub-TeV gamma rays measured by the Tibet AS+MD array. The upcoming TeV–PeV gamma-ray instrument LHAASO (Bai et al. 2019) is promising to test our speculation.

Finally, we caveat that the upper limits of the neutrino flux employed in this work are only in 90% C.L. Therefore, the constraint from the neutrino flux is not so strict and hence the contribution needed to be attributed to extended sources may be lower. Aartsen et al. (2017) showed a 2D likelihood scan of the normalization and spectral slope of Galactic neutrino flux. For the benchmark spectral index, i.e., 2.5, the 99% C.L. flux upper limit is about 1.5 times the 90% C.L. upper limit as employed in this study. Adopting the 99% C.L. upper limit would relax the tension between the diffuse neutrinos and diffuse gamma rays. Also, the analysis for the neutrino upper limits of 2HWC J2031+415 and Gamma Cygni, the two sources related to the Cygnus region, assumes a pointlike source. A dedicated analysis optimized for the extended source could result in a different upper limit. Therefore, we encourage a detailed analysis of the neutrino emission on the Cygnus region with the spatial template of gamma-ray emission above 30 TeV.

We thank Christian Haack for helpful discussion on the IceCube data and the anonymous referee for the constructive comments. This work is supported by NSFC grants U2031105, 11625312, and 11851304, and the National Key R & D program of China under the grant 2018YFA0404203.

We consider a simplified case in which high-energy electrons are injected at a constant rate at r = 0 and diffuse isotropically to a larger radius. The injection spectrum is given by ${Q}_{e,\mathrm{inj}}={Q}_{0}{E}_{e}^{-p}$ for E_e ≤ 2 PeV. The prefactor Q₀ can be found by ∫E_e Q_inj(E_e )dE_e = L_e,inj with L_e,inj being the total electron luminosity, which will be treated as an input parameter. Following Syrovatskii (1959), the electron distribution can be given by

$$\begin{eqnarray}&&\begin{array}{l}N({E}_{e},r)={\displaystyle \int }_{0}^{{t}_{\mathrm{age}}}\displaystyle \frac{{Q}_{e,\mathrm{inj}}({E}_{g}){dt}}{{\left(4\pi \lambda ({E}_{e},t)\right)}^{3/2}}\\ \quad \exp \left[-\displaystyle \frac{{r}^{2}}{4\lambda ({E}_{e},t)}\right]\displaystyle \frac{{{dE}}_{g}}{{{dE}}_{e}},\end{array}\end{eqnarray} \tag{ 6 }$$

where t_age is taken to be 1 Myr and $\lambda ({E}_{e},t)={\int }_{t}^{{t}_{\mathrm{age}}}D({E}_{e}^{\prime} (t^{\prime} )){dt}^{\prime}$. ${E}_{e}^{\prime} (t^{\prime} )$ represents the trajectory in the energy space of an electron, the energy of which is E_e at present, and E_g is the initial energy of the electron at the generation (injection) time t. The relation between E_e and E_g as well as dE_g /dE_e can be obtained by the energy evolution of the electron, i.e.,

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{{dE}}_{e}}{{dt}}=-\displaystyle \frac{4}{3}{\sigma }_{T}c{\left(\displaystyle \frac{{E}_{e}}{{m}_{e}{c}^{2}}\right)}^{2}\\ \quad \times \ \left[{U}_{B}+{{\rm{\Sigma }}}_{i}{U}_{i,\mathrm{ph}}/{\left(1+4\displaystyle \frac{{E}_{e}{\epsilon }_{i}}{{m}_{e}^{2}{c}^{4}}\right)}^{3/2}\right],\end{array}\end{eqnarray} \tag{ 7 }$$

where σ_T is the Thomson cross section, m_e is the electron mass, and c is the speed of light. U_B = B²/8π is the magnetic field energy density. U_i,ph is the energy density of the ith radiation field. _i = 2.82kT_i is the typical photon energy of the ith radiation field of a blackbody or a graybody radiation field with a temperature T_i and k is the Boltzmann constant (Moderski et al. 2005). U_i,ph here includes the CMB radiation and the mean interstellar radiation field of the Galaxy that is composed of an infrared radiation field with T = 30 K and U = 0.3 eVcm⁻³ and an optical radiation field with T = 5000 K and U = 0.3 eVcm⁻³, although the interstellar radiation field has a negligible influence on PeV electrons due to the KN effect. The electron cooling timescale as a function of energy is shown in Figure 5.

Zoom In Zoom Out Reset image size

With the electron's spatial distribution, we can calculate the gamma-ray surface brightness profile I_cyg(E_γ , θ) by integrating the IC radiation over the line of sight toward an arbitrary angle θ from the center of Cygnus Cocoon, following the method detailed in Liu et al. (2019). The total unmasked gamma-ray flux related to Cygnus Cocoon then reads

$$\begin{eqnarray}&&{F}_{\mathrm{cyg},\gamma }({E}_{\gamma })=2\pi {\int }_{{0.5}^{^\circ }}^{\infty }{I}_{\mathrm{cyg}}({E}_{\gamma },\theta )\sin \theta d\theta .\end{eqnarray} \tag{ 8 }$$

The contribution of extended sources in the leptonic case to the DGE can be obtained by Equation (5).