Galactic cosmic ray propagation: sub-PeV diffuse gamma-ray and neutrino emission

Abstract

The Tibet ASγ experiment just reported their measurement of sub-PeV diffuse gamma-ray emission from the Galactic disk, with the highest energy up to 957 TeV. These diffuse gamma rays are most likely the hadronic origin by cosmic ray (CR) interaction with interstellar gas in the galaxy. This measurement provides direct evidence to the hypothesis that the Galactic Cosmic Rays (GCRs) can be accelerated beyond PeV energies. In this work, we try to explain the sub-PeV diffuse gamma-ray spectrum with different CR propagation models. We find that there is a tension between the sub-PeV diffuse gamma-ray and the local CR spectrum. To describe the sub-PeV diffuse gamma-ray flux, it generally requires larger local CR flux than measurement in the knee region. We further calculate the PeV neutrino flux from the CR propagation model. Even all of these sub-PeV diffuse gamma rays originate from the propagation, the Galactic Neutrinos (GNs) only account for less than ~15% of observed flux, most of which are still from extragalactic sources.

Graphical abstract

Keywords

galactic cosmic ray / diffuse gamma ray / neutrino

1 Introduction

In the all-particle spectrum, the most prominent feature is the so-called “knee” structure at ~4 PeV, where the spectrum exhibits a slight steepening, with the change in slope [1-5]. Meanwhile above the knee, a tendency to the heavier nuclei has been noted [6, 7]. It is believed that the cosmic rays (CRs) less than PeV energies are principally of galactic origin. However, which Galactic objects could accelerate PeV CRs, has long been in dispute. The supernova remnants (SNRs) had long been proposed as the sources of Galactic Cosmic Rays (GCRs) [8]. But long gazes into most SNRs associated with the molecular clouds indicate that the cut-off energy in the gamma-ray spectrum is less than 10 TeV, which means the inferred maximum cosmic ray (CR) energy is below ~100 TeV.

Massive strides about the PeV sources have been made in recent years. The Galactic center, which embodies a supermassive black hole, Sagittarius A*, with mass up to 4 × 10⁶ solar mass, is regarded as a PeVatron candidate. The observation of anisotropy phase also seems to hint at this [9-11]. The H.E.S.S. Collaboration performed deep observations of the central molecular zone surrounding the Galactic center. They detected a single power-law gamma-ray spectrum up to tens of TeV, without a break or a cutoff. This partly indicates a PeV proton accelerator at the Galactic center [12]. Early this year, the Tibet AS$\gamma$ experiment reported their observation of SNR G106.3+2.7 [13]. They successfully detected the 100 TeV gamma-ray photons, which are well correlated with the nearby molecular cloud. Therefore the morphological feature favours a hadronic origin of such very-high energy gamma rays. The young massive star clusters are another possible acceleration sites. Most recently, HAWC published their gamma-ray observation of the Cyngnus Cocoon, a well-known superbubble, with energy ranging from 1 to 100 TeV [14]. The 100 TeV gamma-ray photons may originate from the accelerated CRs in the Cyg OB2, an active star-forming region. These recent discoveries provide hard evidence concerning the galactic origin of PeV CRs.

Apart from searching for point sources, the observations of 100 TeV diffuse gamma-ray emission could also be served as the evidence of PeV CRs origin. At such high energy, the gamma-ray photons are expected to have interactions with the low energy background radiation, so that the extragalactic flux is strongly suppressed. On the other hand, the all-electron spectrum shows a clear steepening above 1 TeV, where the spectral index changes to ~ −3.9 [15, 16]. Coupled with the energy loss during propagation of electrons, the 100 TeV gamma-ray flux from the inverse Compton scattering off electrons could be safely neglected. Such energetic gamma rays are clearly generated from the ${\pi }^0$-decay process during the hadronic interaction of the CR nuclei with the interstellar medium (ISM). Previously, only CASA-MIA [17] and KASCADE [18] experiments set upper limits on the diffuse gamma rays above 20 TeV. At present, the Tibet AS$\gamma$ experiment managed to detect the sub-PeV diffuse gamma rays in the Galactic disk, with energy ranging from 100 TeV to 1 PeV [19]. The observed highest energy is up to 957 TeV, very close to 1 PeV. This discovery indicates that the Galactic sources could accelerate CRs to at least ~10 PeV. Some subsequent works had been carried out to study the origin of sub-PeV diffuse gamma-ray emission [20-23].

During the p−p interactions of CR with ISM, the generated gamma-ray photons are simultaneously accompanied with the neutrinos, which are the decay products of the ${\pi }^{\pm }$s. Thereupon, the high-energy neutrinos are also regarded as a good probe to the hadronic interaction. The astrophysical neutrinos have been detected by the IceCube Neutrino Observatory [24-26]. The overall distribution of the neutrino event samples is consistent with an isotropic distribution, which demonstrates that they are held to be predominantly extragalactic. The directional searches also have found an excess from a starburst galaxy [27] and the neutrino emission associated with a blazar [28, 29]. However the proportion of Galactic contribution is still uncertain. A separate fit of the Northern and Southern hemisphere signals in the four-year signal shows a preference to a harder spectrum in the Northern hemisphere [30] which could potentially be due to the presence of a softer contribution of the flux from the inner galaxy in the Southern hemisphere.

The sub-PeV diffuse gamma-ray emission effectively traces the spatial distribution of remote CRs, so they could be well applied to testify the available propagation models. Besides, it is helpful to unveil the origin of knee region [31]. In this work, we investigate the propagation origin of sub-PeV diffuse gamma rays. We find, when the sub-PeV diffuse gamma-ray flux is entirely attributed to p−p interactions of CRs, and, based on current understandings of CR transport and gas distribution, there is a tension between sub-PeV diffuse gamma rays and local CRs measurement. To explain the sub-PeV diffuse gamma rays measured by the Tibet AS$\gamma$ experiment, the calculated local CR flux inevitably exceeds the observed CR flux above PeV energies. One possibility is that the propagated spectrum close to the Galactic center is harder as for the current propagation models, whereas the propagation around the solar system is still unchanged. Meanwhile these sub-PeV gamma-ray photons from point sources, especially unresolved, may not be entirely subtracted. We further evaluate the Galactic contribution to the diffuse neutrino flux and find at most $15\mathrm{\%}$ of observed flux come from GCRs propagation. It should be noted that Lipari et al. [32] have made the prediction for the diffuse gamm rays above hundreds of TeV. But the spatial distribution of CRs is extrapolated according to the local observations. In this work, the whole CR spatial distribution is evaluated by solving the diffusion equation numerically.

2 Model description

2.1 Homogeneous diffusion

In the conventional propagation model, the diffusion process is supposed to be homogeneous and isotropic, so the diffusion coefficient is only a function of rigidity $\mathcal{R}=p/(Ze)$, namely,

(1)$\begin{array}{r}D(\mathcal{R})=D_0{\beta }^{\eta }{\left({\displaystyle \frac{\mathcal{R}}{{\mathcal{R}}_0}}\right)}^{{\delta }_0},\end{array}$

where ${\mathcal{R}}_0=4$ GV, $\beta$ is the particle’s velocity in unit of light speed. $\eta$ is a phenomenological constant in order to fit the low-energy data. $D_0$ and ${\delta }_0$ are constants representing the diffusion coefficient and its high-energy rigidity dependence in the outer halo respectively. Usually, the power index ${\delta }_0$ is taken from $0.3$ to $0.6$, as inferred from the fitting of boron-to-carbon ratio [33]. After propagation, the CR spectrum falls off as a single power-law, $\varphi \propto {\mathcal{R}}^{-\nu -\delta }$, where $\nu$ is the power index of CR spectrum at source.

2.2 Spatial-dependent propagation

However the spectral hardening of CR nuclei above $\sim 200$ GV [34] as well as the anisotropy observation [11] severely challenge the conventional homogeneous propagation model [35, 36]. The spatial-dependent propagation (SDP) was intially introduced to account for the excess of CR nuclei [37]. Afterwards, it is further applied to large-scale anisotropy [38, 39] and diffuse gamma-ray [40] observations. For a comprehensive introduction, one can refer to [41] and [42].

In contrast to the homogeneous diffusion, the entire diffusive halo in the SDP model is split into two zones characterized by diverse diffusion properties, i.e., inner halo (IH) and outer halo (OH). The Galactic disk and its surrounding areas within a few hundred parsecs are called IH, where the diffusion is slower and relevant to the radial distribution of sources. The diffusion in extended regions outside of IH, i.e., OH, is faster and approaches to the conventional propagation. The diffusion coefficient $D$ in the whole region is thus parameterized as

(2)$\begin{array}{r}{D}_{xx}(r,z,\mathcal{R})=D_{0}F(r,z){\beta }^{\eta }{\left({\displaystyle \frac{\mathcal{R}}{{\mathcal{R}}_0}}\right)}^{{\delta }_{0}F(r,z)},\end{array}$

where $r$ and $z$ are cylindrical coordinate, the other parameters are same as the conventional propagation model. The spatial dependence function $F(r,z)$ is shown as follows:

(3)$\begin{array}{r}F(r,z)=\{\begin{array}{ll}g(r,z)+\left[1-g(r,z)\right]{\left({\displaystyle \frac{z}{\xi L}}\right)}^n,& |z|\le \xi L\\ 1,& |z|>\xi L\end{array},\end{array}$

in which $g(r,z)$ is $N_m/[1+f(r,z)]$, and $f(r,z)$ is the source density distribution. $N_m$ is a constant characterized the particle’s diffusion velocity in the inner halo. $L$ is the half-thickness of the propagation cylinder, and $\xi L$ is the half-thickness of the inner halo. The factor ${\left(\frac{z}{\xi L}\right)}^n$ describes the smoothness of the parameters at the transition between the two halos. Note that the spatial dependence of the diffusion coefficient is phenomenologically assumed. Physically it may be related with the magnetic field distribution, or possibly the turbulence driven by CRs [43]. We adopt the diffusion reacceleration model in this work, with the reacceleration being described by a diffusion in the momentum space. The momentum diffusion coefficient, $D_{pp}$, correlates with $D_{xx}$ via $D_{pp}{D}_{xx}=$ $\frac{4p^2{v}_A^2}{3\delta (4-{\delta }^2)(4-\delta )}$, where $v_A$ is the Alfvén velocity, $p$ is the momentum, and $\delta$ is the rigidity dependence slope of the spatial diffusion coefficient [44]. The spatial distribution of CRs is obtained by numerically solving the diffusion equation with the DRAGON package. The diffuse gamma-ray distribution around the Galactic disk is calculated based on the CRs spatial distribution obtained with DRAGON. The employed gas distribution follows the default setup in the GALPROP.

2.3 Background source distribution

Supernova remnants (SNRs) are considered to be the most plausible candidates for the acceleration of CRs. The spatial distribution of SNRs are approximated as an axisymmetric form parameterized as

(4)$\begin{array}{r}f(r,z)={\left({\displaystyle \frac{r}{r_{\odot }}}\right)}^{\alpha }\exp \left[-{\displaystyle \frac{\beta (r-r_{\odot })}{r_{\odot }}}\right]\exp \left(-{\displaystyle \frac{|z|}{z_s}}\right),\end{array}$

where $r_{\odot }\equiv 8.5$ kpc represents the distance from the Galactic center to the solar system. Parameters $\alpha$ and $\beta$ are taken to be $1.69$ and $3.33$ [45]. The density of the SNR distribution decreases exponentially along the vertical height from the Galactic plane, with $z_s=200$ pc.

In the calculations below, the injection spectra of all propagation models are assumed to have a power-law plus a high-energy exponential cutoff, i.e.,

(5)$\begin{array}{r}{q}_i(\mathcal{R})=q_0^i{\left({\displaystyle \frac{\mathcal{R}}{{\mathcal{R}}_{\mathrm{b}\mathrm{r}}}}\right)}^{{\nu }_i}\exp \left[-{\displaystyle \frac{\mathcal{R}}{{\mathcal{R}}_{\mathrm{c}}}}\right],\end{array}$

for the $i$-th composition. where $q_0$ is the normalization factor, $\nu$ is the spectral incides, ${\mathcal{R}}_{\mathrm{b}\mathrm{r}}$ is break rigidity fixed at 2.2 GV, ${\mathcal{R}}_{\mathrm{c}}$ is the cutoff rigidity.

2.4 Nearby source

The propagation of particles from the nearby source is calculated using the Green’s function method, assuming a spherical geometry with infinite boundary conditions. Assuming instantaneous injection from a point source, the CR density as a function of space, rigidity, and time can be calculated as

(6)$\begin{array}{r}\varphi (r,\mathcal{R},t)={\displaystyle \frac{q_{\mathrm{i}\mathrm{n}\mathrm{j}}(\mathcal{R})}{{\left(\sqrt{2\pi }\sigma \right)}^3}}\exp \left(-{\displaystyle \frac{r^2}{2{\sigma }^2}}\right),\end{array}$

where $q_{\mathrm{i}\mathrm{n}\mathrm{j}}(\mathcal{R})$ is the injection spectrum as a function of rigidity, $\sigma (\mathcal{R},t)=\sqrt{2D(\mathcal{R})t}$ is the effective diffusion length within time $t$. The diffusion coefficient $D(\mathcal{R})$ takes the solar system value of Eq. (2). The function form of $q_{\mathrm{i}\mathrm{n}\mathrm{j}}(\mathcal{R})$ is assumed to be power-law with an exponential cutoff, $q_{\mathrm{i}\mathrm{n}\mathrm{j}}(\mathcal{R})=q_0{\mathcal{R}}^{-\alpha }\exp (-\mathcal{R}/{\mathcal{R}}_{\mathrm{c}}^{\prime })$. The normalization $q_0$ is determined through fitting to the CR energy spectra. The distance and age of the local source are set to be $d=330$ pc and $\tau =3.4\times {10}^5$ years [38], respectively.

3 Results

In this work, to make a detailed study of the propagation origin of sub-PeV diffuse gamma rays, we compare two kinds of common propagation scenarios, i.e. homogeneous diffusion (HD) and spatial-dependent propagation (SDP). We have introduced two SDP models according to the origin of spectral harding. In the model SDP-A, the excess of CR nuclei above $200$ GV is regarded as a local effect, which originate from a local SNR. A nearby SNR is also favored in order to describe the evolution of anisotropy amplitude and phase with energy [38, 39, 46]. For model SDP-B, the excesses chiefly originate from the spatial variation of diffusion coefficient. The excess of diffuse gamma-ray emission at Galactic plane above a few GeV has also been well accounted for in model SDP-B [40], due to the spectral hardening of CRs in the whole diffusive halo.

3.1 CR energy spectra

Before calculating the diffuse gamma-ray distribution from the ${\pi }^0$-decay, the CRs spatial distribution in the Galactic halo have to be evaluated by certain propagation and injection parameters. For this purpose, the local observations of CR energy spectrum have to be fitted first, in order to obtain the propagation and injection parameters. In the HD propagation, the essential propagation parameters includes $D_0,{\delta }_0$ and $L$. As for the SDP models, three additional parameters have to be involved, i.e., $N_m,\xi$ and $n$. Fig.1 shows the fitting to the latest B/C ratio published by AMS-02 experiment [47], and the corresponding propagation parameters are listed in Tab.1. Compared with the conventional propagation, the SDP models anticipates a flattening for the B/C ratio above hundreds of GeV.

Fig.1 Calculation of B/C ratio to obtain the propagation parameters in HD, SDP-A and SDP-B models, with B/C data taken from the AMS-02 measurement [47].

Full size|PPT slide

Tab.1 Diffusion parameters of three propagation models.

Model	$D_0({\mathrm{c}\mathrm{m}}^2\cdot {\mathrm{s}}^{-1})$	${\delta }_0$	$N_m$	$\xi$	$n$	$v_A(\mathrm{k}\mathrm{m}\cdot {\mathrm{s}}^{-1})$	$L(\mathrm{k}\mathrm{p}\mathrm{c})$
HD	4.75 × 1028	0.46				22	5
SDP-A	5.64 × 1028	0.56	0.51	0.1	3.5	6	5
SDP-B	5.84 × 1028	0.55	0.49	0.1	3.5	6	5

Given the propagation parameters, we further fit the local CR energy spectra to obtain the source’s injection parameters. The calculated proton and helium spectra are illustrated in Fig.2, both of which are the major elements of CRs less than $10$ PeV. The injection parameters are listed in Tab.2. The green lines are the propagated proton and helium spectra in the HD model. For the HD model, we do not taken into account the local source, since the nuclei excess between $200$ GV and tens of TeV does not affect the sub-PeV gamma-ray photon productivity. Compared with the observations less than 100 TeV, both proton and helium spectra measured by the KASCADE experiment indicate a visible steepening. In this work, the cutoff rigidity of different elements is assumed to have $Z$-dependence. To fit KASCADE observations, the cutoff rigidity of proton and helium are set to $7$ PV for HD and SDP-A model, but $3.5$ PV for SDP-B model.

Fig.2 Calculated proton and helium energy spectra in propagation models of HD, SDP-A and SDP-B. The proton and helium data are taken from AMS-02 [50, 51], CREAM-II [48], NUCLEON [52], DAMPE [49] and KASCADE [53], IceTop [54] and HAWC [55].

Full size|PPT slide

Tab.2 Injection parameters of the background and local sources in three propagation models.

Model	Parameters	p	He	C	N	O	Ne	Mg	Si	Fe
HD	Normalization*	4420	401	16.8	2.26	19.3	2.24	2.78	4.20	2.92
	${10}^{-5}[({\mathrm{m}}^2\cdot \mathrm{s}\mathrm{r}\cdot \mathrm{s}\cdot \mathrm{G}\mathrm{e}\mathrm{V}{)}^{-1}]$
	$\nu$	2.32	2.26	2.31	2.37	2.32	2.31	2.31	2.41	2.27
	${\mathcal{R}}_c(\mathrm{P}\mathrm{V})$	7	7	7	7	7	7	7	7	7
SDP-A background	Normalization	4090	259		1.40	13.3	1.50	1.86	2.05	2.21
	${10}^{-5}[({\mathrm{m}}^2\cdot \mathrm{s}\mathrm{r}\cdot \mathrm{s}\cdot \mathrm{G}\mathrm{e}\mathrm{V}{)}^{-1}]$
	$\nu$	2.37	2.29	2.32	2.37	2.34	2.32	2.33	2.45	2.32
	${\mathcal{R}}_c(\mathrm{P}\mathrm{V})$	7	7	7	7	7	7	7	7	7
SDP-A local	$q_0$	2800	1400	36	5.2	40	6.2	6.6	5.7	5.2
	${10}^{49}(\mathrm{G}\mathrm{e}{\mathrm{V}}^{-1})$
	$\alpha$	2.10	2.10	2.05	2.05	2.05	2.05	2.05	2.05	2.05
	${\mathcal{R}}_c^{\prime }(\mathrm{T}\mathrm{V})$	28	28	28	28	28	28	28	28	28
SDP-B	Normalization	4500	282	12.8	1.84	15.8	1.93	2.59	2.64	2.89
	${10}^{-5}[({\mathrm{m}}^2\cdot \mathrm{s}\mathrm{r}\cdot \mathrm{s}\cdot \mathrm{G}\mathrm{e}\mathrm{V}{)}^{-1}]$
	$\nu$	2.34	2.26	2.33	2.39	2.34	2.35	2.37	2.41	2.35
	${\mathcal{R}}_c(\mathrm{P}\mathrm{V})$	3.5	3.5	3.5	3.5	3.5	3.5	3.5	3.5	3.5

*The normalization is set at kinetic energy per nucleon $E_{\mathrm{k}}=100$ GeV/n.

The SDP-A model is marked as the red lines, in which the dash-dot, dash and solid lines denote the fluxes from background sources, local SNR and sum of them respectively. To reproduce the softening at 20 TeV [48, 49], the cutoff rigidities of local CRs are set to 28 TV. It has been shown that this break is relevant to the peak in the anisotropy amplitude at $\sim 10$ TeV and the flip of anisotropy phase at $\sim 100$ TeV. In this SDP model, energetic CRs principally propagate within the IH. Therefore compared with HD, the background spectra (i.e., red dash-dot lines) gradually harden above ${10}^4$ GeV, so that the calculated fluxes above ${10}^4$ GeV are higher than HD model. The corresponding sub-PeV diffuse gamma-ray flux is expected to be higher in the SDP model.

The blue lines are the calculated proton and helium fluxes in SDP-B model. Here the parameter $N_m$ is $0.49$, which is smaller than $0.51$ of SDP-A model, so the propagated spectrum and corresponding B/C ratio are harder than those of SDP-A model. The cutoff rigidity is set to $3.5$ PV in order to fit the KASCADE data. But the observed softening structure at around tens of TeV could not be explained, and both proton and helium fluxes inevitably exceed those of HD and SDP-A models above $100$ TeV. As we will show below, these extra fluxes are important to account for the observed sub-PeV diffuse gamma rays.

3.2 All-particle spectrum

The diffuse sub-PeV photons principally originate from the CRs at PeV energies, and an amount of ground-based observations of all-particle spectrum have been performed at this energy range with high precision. To calculate all-particle spectrum, we fit other major components, namely C, N, O, Ne, Mg, Si and Fe respectively. All of them conform with the observations at lower energy separately and are extrapolated the model to the knee region. The injection parameters are summarized in Tab.2. The cutoff rigidity of nuclei is $Z$-dependent in order to fit the knee region. It is worth noting that at knee region, which is most relevant to the Tibet sub-PeV gamma-ray observation, the major composition of CR flux is proton and helium, so the fluxes of heavier composition do not take leading role.

In Fig.3, we compare the calculated all-particle spectra with observations. We could see that the expected all-particle spectrum in HD model is lower than the observations by simply extrapolating the low-energy spectrum assuming a single power-law injection spectrum. This would make the expected diffuse sub-PeV flux well below the observations, as we would see. The all-particle spectrum could be accounted for until $\sim 10$ PeV in both SDP models. Beyond that energy, the Galactic contribution drops significantly, and these CRs may originate from the extra-galactic sources. Furthermore the CRs above $10$ PeV essentially do not contribute significantly to the sub-PeV gamma-ray photons during the interaction with the interstellar medium. Therefore we do not further study the origin of such high energy CRs. It is worth noting that the all-particle flux in SDP-B model approaches to the upper limits of measurements above PeV energies.

Fig.3 Calculated all-particle spectra in three propagation models. The all-particle data are taken from Horandel [56], IceTop [54, 57], Tibet [5], HAWC [58], ARGO [59], ATIC [60], NUCLEON [61] and KASCADE [4].

Full size|PPT slide

3.3 Diffuse gamma-ray emission

After reproducing the local CR observations, we further calculate the spatial distribution of CRs and the corresponding diffuse gamma-ray spectra. Above $\sim 20$ TeV, the gamma-ray photons would have interactions with low-energy Galactic interstellar radiation field (ISRF) and the observed gamma-ray flux is expected to be strongly suppressed [62, 63]. Fig.4 shows the diffuse gamma-ray spectra calculated by the three propagation models. The solid and dash lines are the diffuse gamma-ray spectra with/without gamma-ray attenuation from pair production respectively. As illustrated in the figure, above $100$ TeV, the gamma-ray fluxes show observable attenuation.

Fig.4 Calculated diffuse gamma-ray spectra in three propagation models. The gamma ray data are taken from ARGO-YBJ [64] and Tibet AS+MD [19] experiments.

Full size|PPT slide

In HD model, the calculated diffuse gamma-ray fluxes are well below the observations. Compared with HD, the gamma-ray spectra above 1 TeV show noticeable bumps in both SDP models so that the flux has been greatly enhanced correspondingly. This results from the spectral hardening of CRs between $100$ TeV and $10$ PeV in the entire halo. Meanwhile compared with SDP-B model, the diffuse gamma-ray flux in SDP-A model is lower, which means the CR flux in the whole Galaxy is always lower than SDP-B model above $100$ TeV. This is because the excess of CR flux become softening at $\sim 20$ TeV in SDP-A model, and is attributed to a local effect. Especially close to the direction of anti-Galactic center, i.e., ${50}^{\circ }<l<{200}^{\circ }$, the gamma-ray flux above $100$ TeV is significantly lower than Tibet observation. Furthermore, the flux at 1 TeV is slightly lower than ARGO measurements.

Only the expectations by the SDP-B model could marginally reaches the lower limits of both ARGO-YBJ and Tibet AS+MD observations at ${25}^{\circ }<l<{100}^{\circ }$. This is due to the enhanced CR flux from $10$ TeV to $10$ PeV. Meanwhile in SDP-B model, the hardening is regarded as the propagation effect, which enables the CR flux above $10$ PeV in the entire diffusive halo to be augmented overall. Close to the Galactic center, the propagated CR energy spectra is harder. But as we have seen that, the CR proton flux exceeds the available observations and the softening at tens of TeV could not be reproduced.

As can be seen, the origin of sub-PeV diffuse gamma-ray flux has a tension with local CR observations more or less. When the local CR observations are reproduced, the calculated gamma-ray flux is inadequate to explain the Tibet observations. But if the sub-PeV gamma-ray flux is accounted for, the required CR flux surely exceeds the local CR energy spectra. This is obvious at ${50}^{\circ }<l<{200}^{\circ }$, which requires an unexpectedly large CR flux at anti-Galactic direction.

3.4 Diffuse neutrino spectrum

Based on the sub-PeV gamma-ray observation, we could reckon the Galactic contribution below PeV. On average, the p−p collision produces nearly one-third neutral pions and two-thirds charged pions. Each neutral pion decays into a pair of gamma rays, whereas each charged pion decays into two muon neutrinos and one electron neutrino (here we do not distinguish between neutrinos and anti-neutrinos). The initial neutrino flavor ratio is approximately ${\nu }_e:{\nu }_{\mu }:{\nu }_{\tau }=1:2:0$ from charged pion decay. After travelling, The flavor ratio is transformed to ${\nu }_e:{\nu }_{\mu }:{\nu }_{\tau }=1:1:1$ due to the vacuum neutrino oscillation. The typical neutrino energy from charged pion decay is approximately half of the gamma-ray photon from neutral pion decay. The gamma-ray spectrum at source is $\text{d}{\mathrm{\Phi }}_{\gamma }/\text{d}{E}_{\gamma }={\varphi }_{\gamma }{E}_{\gamma }^{-\mathrm{\Gamma }}$, the resulting neutrino spectrum is shifted relative to the gamma-ray spectrum [65], i.e.,

(7)$\begin{array}{r}{\displaystyle \frac{\text{d}{\mathrm{\Phi }}_{\nu }}{\text{d}{E}_{\nu }}}={\varphi }_{\nu }{E}_{\nu }^{-\mathrm{\Gamma }}={\left({\displaystyle \frac{1}{2}}\right)}^{\mathrm{\Gamma }-1}{\varphi }_{\gamma }{E}_{\nu }^{-\mathrm{\Gamma }}.\end{array}$

In Fig.5, we evaluate the diffuse neutrino flux in three propagation models respectively. Compared with the fitting result of IceCube observations, which has a power index of $\sim 2.87$, the spectrum of galactic diffuse neutrinos is softer above hundreds of TeV. Because of the cut-off energy of background injection spectrum of SDP-B model is smaller than another two propagation models, it can be seen that the calculated diffuse neutrino flux of SDP-B model is lower than that of the other two models above $\sim 500$ TeV. While it is larger than the other two models with energy less than $\sim 500$ TeV. Even in this case, the galactic diffuse neutrinos generated during CRs propagation, only take up to at most $15\mathrm{\%}$ of latest observations and are less than the $90\mathrm{\%}$ C.L. upper limit for the neutrino flux of the Galactic plane measured by the ANTARES neutrino telescope and IceCube.

Fig.5 Diffuse neutrino flux calculated by the three propagation models. The data are taken from the ICECUBE $7.5$ years’ observation [66]. The violet line is the power-law fitting to the data, with normalization $\mathrm{\Phi }=6.37\times {10}^{18}$ GeV−1·cm−2·s−1·sr−1 at $100$ TeV and power index $\gamma =-2.87$. The orange bar with arrows exhibits the 90% C.L. upper limit for the neutrino flux of the Galactic plane using the ANTARES neutrino telescope[67], while the gray bar with arrows is the same as above except for using the IceCube 7 years data [68].

Full size|PPT slide

4 Conclusion

The argument that the CRs below knee region are accelerated by the Galactic sources, has long been the lack of clear evidence. Most recently, the Tibet AS$\gamma$ experiment reported their findings of the diffuse gamma ray photons between 100 TeV and 1 PeV, nearby the Galactic disk. These sub-PeV photons are likely the ${\pi }^0$-decay products from the hadronic interactions between the CR nuclei and interstellar materials. The Galactic origin of CRs with energy up to PeV have a definite conclusion along with observations of Galactic center [12, 69], SNR G106.3+2.7 [13, 70], the massive young cluster Westerlund 1 [69, 71], Cygnus OB2 [69], gamma-ray binary source (LS 5039, PSR B1259-63, LSI +61°303, PSR J2032+4127) [72] as well as Cyngnus Cocoon [14] and so on.

In this work, we study the propagation origin of these sub-PeV diffuse gamma rays by considering three competing propagation models. One is homogeneous diffusion and the other two are SDP models. Compared with the homogeneous diffusion, the SDP-A model could generate more gamma rays above $10$ TeV around the disk due to the flattenning propagated CR spectra above tens of TeV. Nonetheless, to meet the observed sub-PeV gamma-ray flux at ${25}^{\circ }<l<{100}^{\circ }$, the local CR flux in the SDP plus local source model has to slightly exceed the observation of the all-particle spectrum at around $10$ PeV. And the calculated gamma rays are still inadequate for the measurements of 1 TeV at ${25}^{\circ }<l<{100}^{\circ }$ and above $100$ TeV at ${50}^{\circ }<l<{200}^{\circ }$. As for the SDP-B model, all above gamma-ray observations could be accounted for. However, such a model meets up with severe challenges. The local CR flux would be greatly enhanced, which far exceeds the all-particle measurements. Meanwhile the softening of CR nuclei at $\sim 20$ TeV could not be reproduced, much less the observations of large-scale anisotropies.

In the current propagation models, there is a tension between sub-PeV gamma rays observation and local CRs measurements. One possibility is that the propagated spectrum is harder close to the Galactic center, whereas the propagation around the solar system is still unchanged, as indicated by the analysis of [73] for the radial distribution of diffuse gamma-ray. Therefore, the observed sub-PeV gamma rays and local CRs could be simultaneously satisfied by the SDP model. Meanwhile, these sub-PeV gamma-ray photons from the Galactic disk may not entirely come from the CR propagation. Most of CR sources and interstellar medium are located around the Galactic disk. There may be some of gamma-ray photons from CR sources, despite that the events within ${0.5}^{\circ }$ from the known TeV sources have been subtracted. This has also been confirmed by the Tibet observation in fact. Above 398 TeV, 4 events at ${50}^{\circ }<l<{200}^{\circ }$ are detected less than 4° from the center of the Cygnus cocoon, which is just proved to be a PeV source by HAWC experiment. We hope more extensive research of diffuse gamma rays and the observation of single CR composition between $100$ TeV and $10$ PeV with high precision could testify our conclusion.

We evaluate the corresponding neutrino flux from the Galactic comic ray propagation. We find that even if the observed sub-PeV diffuse gamma rays could be accounted for, the Galactic plane could contribute $\sim 15\mathrm{\%}$ of observed neutrino flux.

References

[1] M. Nagano, T. Hara, Y. Hatano, N. Hayashida, S. Kawaguchi, K. Kamata, T. Kifune, Y. Mizumoto. Energy spectrum of primary cosmic rays between 10^14.5 and 10¹⁸ eV. J. Phys. G Nucl. Phys., 1984, 10(9): 1295

[2] M. A. K. Glasmacher, M. A. Catanese, M. C. Chantell. . The cosmic ray energy spectrum between 10¹⁴ and 10¹⁶ eV. Astropart. Phys., 1999, 10(4): 291

[3] M. Aglietta, B. Alessandro, P. Antonioli, F. Arneodo, L. Bergamasco. . The cosmic ray primary composition in the “knee” region through the EAS electromagnetic and muon measurements at EAS-TOP. Astropart. Phys., 2004, 21(6): 583

[4] T.AntoniW. D. ApelA.F. BadeaK.BekkA.Bercuci, ., KASCADE measurements of energy spectra for elemental groups of cosmic rays: Results and open problems, Astropart. Phys. 24(1–2), 1 (2005)

[5] M. Amenomori, X. J. Bi, D. Chen, S. W. Cui. . The all-particle spectrum of primary cosmic rays in the wide energy range from 10¹⁴ to 10¹⁷ eV observed with the Tibet-III air-shower array. Astrophys. J., 2008, 678(2): 1165

[6] K. H. Kampert, M. Unger. Measurements of the cosmic ray composition with air shower experiments. Astropart. Phys., 2012, 35(10): 660

[7] R.AloisioP. BlasiI.De MitriS.Petrera, Selected Topics in Cosmic Ray Physics, page 1, 2018

[8] W. Baade, F. Zwicky. Cosmic rays from supernovae. Contributions from the Mount Wilson Observatory, 1934, 3: 79

[9] R. Abbasi, Y. Abdou, T. Abu-Zayyad, M. Ackermann, J. Adams. . Observation of anisotropy in the galactic cosmic-ray arrival directions at 400 TeV with IceCube. Astrophys. J., 2012, 746(1): 33

[10] M. G. Aartsen, K. Abraham, M. Ackermann, J. Adams, J. A. Aguilar. . Anisotropy in cosmic-ray arrival directions in the southern hemisphere based on six years of data from the IceCube detector. Astrophys. J., 2016, 826(2): 220

[11] M. Amenomori, X. J. Bi, D. Chen, T. L. Chen, W. Y. Chen. . Northern sky galactic cosmic ray anisotropy between 10 and 1000 TeV with the Tibet air shower array. Astrophys. J., 2017, 836(2): 153

[12] HESS Collaboration. Acceleration of petaelectronvolt protons in the Galactic Centre. Nature, 2016, 531(7595): 476

[13] TibetASγ CollaborationM.AmenomoriY.W. Bao, ., Potential PeVatron supernova remnant G106.3+2.7 seen in the highest-energy gamma rays, Nat. Astron. (2021)

[14] A.U. AbeysekaraA.AlbertR.Alfaro, ., HAWC observations of the acceleration of very-high-energy cosmic rays in the Cygnus Cocoon, arXiv: 2103.06820 (2021)

[15] DAMPE Collaboration. Direct detection of a break in the teraelectronvolt cosmic-ray spectrum of electrons and positrons. Nature, 2017, 552(7683): 63

[16] D. Kerszberg for the H.E.S.S. Collaboration, The cosmic-ray electron spectrum measured with H.E.S.S. 2017

[17] A. Borione, M. A. Catanese, M. C. Chantell, C. E. Covault, J. W. Cronin, B. E. Fick, L. F. Fortson, J. Fowler, M. A. K. Glasmacher, K. D. Green, D. B. Kieda, J. Matthews, B. J. Newport, D. Nitz, R. A. Ong, S. Oser, D. Sinclair, J. C. van der Velde. Constraints on gamma-ray emission from the galactic plane at 300 TeV. Astrophys. J., 1998, 493(1): 175

[18] W.D. ApelJ. C. Arteaga-VelázquezK.Bekk, ., KASCADE-Grande limits on the isotropic diffuse gamma-ray flux between 100 TeV and 1 EeV, arXiv: 1710.02889 (2017)

[19] M. Amenomori, Y. W. Bao, X. J. Bi, D. Chen, T. L. Chen. . First detection of sub-PeV diffuse gamma rays from the galactic disk: Evidence for ubiquitous galactic cosmic rays beyond PeV energies. Phys. Rev. Lett., 2021, 126(14): 141101

[20] R.-Y.LiuX.-Y. Wang, Origin of galactic sub-PeV diffuse gamma-ray emission: Constraints from high-energy neutrino observations, arXiv: 2104.05609 (2021)

[21] V.VecchiottiF.ZuccariniF.L. VillanteG.Pagliaroli, Unresolved sources naturally contribute to PeV γ-ray diffuse emission observed by Tibet ASγ, arXiv: 2107.14584 (2021)

[22] S. Koldobskiy, A. Neronov, D. Semikoz. Pion decay model of the Tibet-AS γ PeV gamma-ray signal. Phys. Rev. D, 2021, 104(4): 043010

[23] P. Zhang, B. Qiao, Q. Yuan, S. Cui, Y. Guo. Ultrahigh-energy diffuse gamma ray emission from cosmic-ray interactions with the medium surrounding acceleration sources. Phys. Rev. D, 2022, 105(2): 023002

[24] IceCube Collaboration. Evidence for high-energy extraterrestrial neutrinos at the IceCube detector. Science, 2013, 342(6161): 1242856

[25] M. G. Aartsen, R. Abbasi, Y. Abdou, M. Ackermann, J. Adams. . First observation of PeV-energy neutrinos with IceCube. Phys. Rev. Lett., 2013, 111(2): 021103

[26] M. G. Aartsen, M. Ackermann, J. Adams, J. A. Aguilar, M. Ahlers. . Observation of high-energy astrophysical neutrinos in three years of IceCube data. Phys. Rev. Lett., 2014, 113(10): 101101

[27] M. G. Aartsen, M. Ackermann, J. Adams, J. A. Aguilar, M. Ahlers. . Time-integrated neutrino source searches with 10 years of IceCube data. Phys. Rev. Lett., 2020, 124(5): 051103

[28] M. Aartsen, M. Ackermann, J. Adams, J. A. Aguilar, M. Ahlers. . Multimessenger observations of a flaring blazar coincident with high-energy neutrino IceCube-170922A. Science, 2018, 361(6398): eaat1378

[29] M. Aartsen, M. Ackermann, J. Adams, J. A. Aguilar, M. Ahlers. . Neutrino emission from the direction of the blazar TXS 0506+056 prior to the IceCube-170922A alert. Science, 2018, 361(6398): 147

[30] M.G. AartsenK.AbrahamM.Ackermann, ., A combined maximum-likelihood analysis of the high energy astrophysical neutrino flux measured with IceCube, arXiv: 1507.03991 (2015)

[31] Y. Q. Guo, H. B. Hu, Q. Yuan, Z. Tian, X. J. Gao. Pinpointing the knee of cosmic rays with diffuse PeV γ-rays and neutrinos. Astrophys. J., 2014, 795(1): 100

[32] P. Lipari, S. Vernetto. Diffuse galactic gamma-ray flux at very high energy. Phys. Rev. D, 2018, 98(4): 043003

[33] Q. Yuan, S. J. Lin, K. Fang, X. J. Bi. Propagation of cosmic rays in the AMS-02 era. Phys. Rev. D, 2017, 95(8): 083007

[34] O.AdrianiG. C. BarbarinoG.A. BazilevskayaR.BellottiM.Boezio, ., PAMELA measurements of cosmic-ray proton and helium spectra, Science 332(6025), 69 (2011)

[35] P.BlasiE. Amato, Diffusive propagation of cosmic rays from supernova remnants in the Galaxy (II): anisotropy, arXiv: 1105.4529 (2011)

[36] W. Liu, X. J. Bi, S. J. Lin, B. B. Wang, P. F. Yin. Excesses of cosmic ray spectra from a single nearby source. Phys. Rev. D, 2017, 96(2): 023006

[37] N. Tomassetti. Origin of the cosmic-ray spectral hardening. Astrophys. J. Lett., 2012, 752(1): L13

[38] W. Liu, Y.-Q. Guo, Q. Yuan. Indication of nearby source signatures of cosmic rays from energy spectra and anisotropies. J. Cosmol. Astropart. Phys., 2019, 10: 010

[39] B.-Q. Qiao, W. Liu, Y.-Q. Guo, Q. Yuan. Anisotropies of different mass compositions of cosmic rays. J. Cosmol. Astropart. Phys., 2019, 12: 007

[40] Y. Q. Guo, Q. Yuan. Understanding the spectral hardenings and radial distribution of Galactic cosmic rays and Fermi diffuse γ rays with spatially-dependent propagation. Phys. Rev. D, 2018, 97(6): 063008

[41] Y. Q. Guo, Z. Tian, C. Jin. Spatial-dependent propagation of cosmic rays results in the spectrum of proton, ratios of P/P, and B/C, and anisotropy of nuclei. ApJ, 2016, 819: 54

[42] W. Liu, Y. Yao, Y. Q. Guo. Revisiting the spatially dependent propagation model with the latest observations of cosmic-ray nuclei. Astrophys. J., 2018, 869(2): 176

[43] P. Blasi, E. Amato, P. D. Serpico. Spectral breaks as a signature of cosmic ray induced turbulence in the galaxy. Phys. Rev. Lett., 2012, 109(6): 061101

[44] E. S. Seo, V. S. Ptuskin. Stochastic reacceleration of cosmic rays in the interstellar medium. Astrophys. J., 1994, 431: 705

[45] G. Case, D. Bhattacharya. Revisiting the galactic supernova remnant distribution. Astron. Astrophys. Suppl. Ser., 1996, 120: 437

[46] M. Ahlers. Deciphering the dipole anisotropy of galactic cosmic rays. Phys. Rev. Lett., 2016, 117(15): 151103

[47] M. Aguilar, L. Ali Cavasonza, G. Ambrosi, L. Arruda, N. Attig. . Precision measurement of the boron to carbon flux ratio in cosmic rays from 1.9 GV to 2.6 TV with the alpha magnetic spectrometer on the international space station. Phys. Rev. Lett., 2016, 117(23): 231102

[48] Y. S. Yoon, T. Anderson, A. Barrau, N. B. Conklin, S. Coutu. . Proton and helium spectra from the CREAM-III flight. Astrophys. J., 2017, 839(1): 5

[49] Q. An, R. Asfandiyarov, P. Azzarello, P. Bernardini, X. J. Bi. . Measurement of the cosmic ray proton spectrum from 40 GeV to 100 TeV with the DAMPE satellite. Sci. Adv., 2019, 5(9): eaax3793

[50] M. Aguilar, D. Aisa, B. Alpat, A. Alvino, G. Ambrosi. . Precision measurement of the proton flux in primary cosmic rays from rigidity 1 GV to 1.8 TV with the alpha magnetic spectrometer on the international space station. Phys. Rev. Lett., 2015, 114(17): 171103

[51] M. Aguilar, L. Ali Cavasonza, B. Alpat, G. Ambrosi, L. Arruda. . Observation of the identical rigidity dependence of He, C, and O cosmic rays at high rigidities by the alpha magnetic spectrometer on the international space station. Phys. Rev. Lett., 2017, 119(25): 251101

[52] E. Atkin, V. Bulatov, V. Dorokhov. . First results of the cosmic ray NUCLEON experiment. J. Cosmol. Astropart. Phys., 2017, 07: 020

[53] W.D. ApelJ. C. Arteaga-VelázquezK.Bekk M.BertainaJ. Blümer, ., KASCADE-Grande measurements of energy spectra for elemental groups of cosmic rays, Astropart. Phys. 47, 54 (2013)

[54] M. G. Aartsen, M. Ackermann, J. Adams, J. A. Aguilar, M. Ahlers. . Cosmic ray spectrum and composition from PeV to EeV using 3 years of data from IceTop and IceCube. Phys. Rev. D, 2019, 100(8): 082002

[55] J. C. Arteaga-Velázquez, HAWC measurements of the energy spectra of cosmic ray protons, helium and heavy nuclei in the TeV range, arXiv: 2108.03208 (2021)

[56] J. R. Hörandel. On the knee in the energy spectrum of cosmic rays. Astropart. Phys., 2003, 19(2): 193

[57] M. G. Aartsen, R. Abbasi, M. Ackermann, J. Adams, J. A. Aguilar. . Cosmic ray spectrum from 250 TeV to 10 PeV using IceTop. Phys. Rev. D, 2020, 102(12): 122001

[58] R. Alfaro, C. Alvarez, J. D. Álvarez, R. Arceo, J. C. Arteaga-Velázquez. . All particle cosmic ray energy spectrum measured by the HAWC experiment from 10 to 500 TeV. Phys. Rev. D, 2017, 96(12): 122001

[59] G.Di Sciascio, Measurement of the cosmic ray energy spectrum with ARGO-YBJ, arXiv: 1408.6739 (2014)

[60] A. D. Panov, J. H. Jr Adams, H. S. Ahn, G. L. Bashinzhagyan, J. W. Watts, J. P. Wefel, J. Wu, O. Ganel, T. G. Guzik, V. I. Zatsepin, I. Isbert, K. C. Kim, M. Christl, E. N. Kouznetsov, M. I. Panasyuk, E. S. Seo, N. V. Sokolskaya, J. Chang, W. K. H. Schmidt, A. R. Fazely. Energy spectra of abundant nuclei of primary cosmic rays from the data of ATIC-2 experiment: Final results. Bull. Russ. Acad. Sci. Phys., 2009, 73(5): 564

[61] E. V. Atkin, V. L. Bulatov, O. A. Vasiliev, A. G. Voronin, N. V. Gorbunov. . Energy spectra of cosmic-ray protons and nuclei measured in the NUCLEON experiment using a new method. Astron. Rep., 2019, 63(1): 66

[62] J.-L. Zhang, X.-J. Bi, H.-B. Hu. Very high energy γ ray absorption by the galactic interstellar radiation field. Astron. Astrophys., 2006, 449(2): 641

[63] I. V. Moskalenko, T. A. Porter, A. W. Strong. Attenuation of very high energy gamma rays by the Milky way interstellar radiation field. Astrophys. J., 2006, 640(2): L155

[64] B. Bartoli, P. Bernardini, X. J. Bi, P. Branchini, A. Budano. . Study of the diffuse gamma-ray emission from the galactic plane with ARGO-YBJ. Astrophys. J., 2015, 806(1): 20

[65] M. D. Kistler, J. F. Beacom. Guaranteed and prospective galactic TeV neutrino sources. Phys. Rev. D, 2006, 74(6): 063007

[66] R.AbbasiM. AckermannJ.Adams, ., The IceCube high-energy starting event sample: Description and flux characterization with 7.5 years of data, arXiv: 2011.03545 (2020)

[67] S. Adrián-Martínez, A. Albert, M. André, M. Anghinolfi, G. Anton. . Constraints on the neutrino emission from the galactic ridge with the ANTARES telescope. Phys. Lett. B, 2016, 760: 143

[68] M. G. Aartsen, M. Ackermann, J. Adams, J. A. Aguilar, M. Ahlers. . Constraints on galactic neutrino emission with seven years of IceCube data. Astrophys. J., 2017, 849(1): 67

[69] F. Aharonian, R. Yang, E. de Oña Wilhelmi. Massive stars as major factories of galactic cosmic rays. Nat. Astron., 2019, 3(6): 561

[70] P. Cristofari. The hunt for pevatrons: The case of supernova remnants. Universe, 2021, 7(9): 324

[71] A. M. Bykov, D. C. Ellison, P. E. Gladilin, S. M. Osipov. Ultrahard spectra of PeV neutrinos from super novae in compact star clusters. Mon. Not. R. Astron. Soc., 2015, 453(1): 113

[72] A. M. Bykov, A. E. Petrov, M. E. Kalyashova, S. V. Troitsky. PeV photon and neutrino flares from galactic gamma-ray binaries. Astrophys. J. Lett., 2021, 921(1): L10

[73] R. Yang, F. Aharonian, C. Evoli. Radial distribution of the diffuse γ-ray emissivity in the galactic disk. Phys. Rev. D, 2016, 93(12): 123007

[74] A. W. Strong, I. V. Moskalenko. Propagation of cosmic-ray nucleons in the galaxy. Astrophys. J., 1998, 509(1): 212

[75] A. W. Strong, I. V. Moskalenko, O. Reimer. Diffuse continuum gamma rays from the galaxy. Astrophys. J., 2000, 537(2): 763

[76] C. Evoli, D. Gaggero, D. Grasso, L. Maccione. Cosmic ray nuclei, antiprotons and gamma rays in the galaxy: A new diffusion model. J. Cosmol. Astropart. Phys., 2008, 10: 018

[77] C. Evoli, D. Gaggero, A. Vittino. . Cosmic-ray propagation with DRAGON2 (I): Numerical solver and as trophysical ingredients. J. Cosmol. Astropart. Phys., 2017, 02: 015

Acknowledgements

This work was supported by the National Key Research and Development Program of China (No. 2016YFA0400200), the National Natural Science Foundation of China (Nos. U1738209, 11875264, 11635011, and U2031110). Software: GALPROP ([74, 75]) available at https://galprop.stanford.edu. DRAGON ([76, 77]) available at https://github.com/cosmicrays.

RIGHTS & PERMISSIONS

2022 Higher Education Press