Neutrino yield from Galactic cosmic rays

# Neutrino yield from Galactic cosmic rays

## Abstract

We calculate the neutrino yield from collisions of cosmic ray (CR) nuclei on gas using the event generator QGSJET-II. We present first the general characteristics and numerical results for the neutrino yield assuming power-law fluxes for the primary CR nuclei. Then we use three parameterisations for the Galactic CR flux to derive the neutrino yield for energies around and above the knee. The shape and the normalization of the resulting neutrino flux above  eV depends on the composition of the Galactic CR flux employed, but is generally dominated by its proton component. The spectral shape and magnitude of the neutrino flux suggest that the IceCube excess is not connected to interactions of Galactic sea CRs. If a fraction of these events has a Galactic origin, then they may be caused by CR overdensities around recent close-by Galactic sources.

###### pacs:
98.70.Sa, 98.35.Eg, 98.70.Rz

## I Introduction

The IceCube Collaboration recently announced evidence for the first detection of extraterrestrial neutrinos at the confidence level (1). This announcement followed the observation of two PeV neutrino cascades (2). The combined data set consists of 28 events with detected energy in the range between 30 TeV and 2 PeV, while background events are expected from atmospheric muons and neutrinos. This excess of events (denoted “IceCube excess” in the following) is consistent with a diffuse intensity (summed over flavors) at the level of

 E2νIν≃(36±12)×eVcm−2s−1sr−1, (1)

based on 17 events in the 60 TeV to 2 PeV energy range.

The implications of this ground-breaking discovery have been discussed widely (4); (3); (5); (6); (7); (8). In particular, the non-observation of events beyond 2 PeV has been interpreted as a break or an exponential cutoff in the neutrino flux, but a single power-law with is compatible with the data (3). Such a suppression, if present, is not typical for many models of extragalactic neutrino sources: For instance, if active galactic nuclei or -ray bursts were sources of ultrahigh energy cosmic rays, then they should produce neutrinos with energies up to  eV. The same holds for the cosmogenic neutrino flux which peaks around  eV for proton primaries (8); (9). Note also that the intensity (1) practically saturates the cascade bound on extragalactic neutrino intensity derived in Ref. (9) from the Fermi-LAT observations of the diffuse gamma-ray background1.

Alternatively, the neutrino events may have a Galactic origin. The neutrino excess has been associated e.g. with unidentified TeV -ray sources (6), CR pevatrons (7) or PeV dark matter (10). The Galactic neutrino flux contains a guaranteed component which is produced by CRs interacting with gas during their confinement in a CR halo. This minimal Galactic neutrino flux has been discussed since the late 1970s (11).

The maximal energy of neutrinos produced via pion production in proton-photon () or proton-gas () interactions is approximately 10% of the energy of the proton primary. Thus the observed 2 PeV neutrino event requires proton energies above 20 PeV. However, already at  eV protons represent only a subdominant fraction of the primary CR flux compared to helium and heavier nuclei (12); (13). Moreover, the composition of the CR flux becomes increasingly heavier in the energy range between the knee at  PeV and  eV (12); (14); (15). Since the maximal neutrino energy in nucleus-proton collisions is a factor lower than in interactions ( being the nuclear mass number), the required minimal CR energy to explain the IceCube events increases compared to processes. This implies in turn that the number of potential scattering events, and thus secondary fluxes, is drastically reduced, because the CR spectrum is steeply falling, above the knee (12). Therefore it is essential to account correctly for the elemental composition of the Galactic CR flux, if one aims at relating the neutrino intensity required to explain the IceCube excess to the primary CR intensity (16).

Aim of this work is to quantify the neutrino yield from nucleus-proton collisions using up-to-date simulation tools for the relevant hadron production processes and including information on the elemental composition of the CR flux around and above the knee region. Our simulations are based on the event generator QGSJET-II-04 (17) which includes relevant experimental information from run I of LHC (18). We present in Sec. II the general characteristics of the neutrino yields for the case of power-law fluxes of primary nuclei. Then we calculate in Sec. III the neutrino yields for three parametrisations of the CR flux, suitable for the energy region of the knee. Finally, we comment in Sec. IV on the IceCube neutrino excess in view of our findings before we conclude.

## Ii Neutrino yield for power-law CR fluxes

The neutrino intensity produced by CR interactions with the interstellar gas is given by

 Iν(E)=~εMdCRngas∑i∫∞AiEdE′Ii(E′)σinelAip(E′/Ai) ×dnAip→ν(E′/Ai,E)dE, (2)

where the sum goes over the primary CR mass groups , is the partial intensity for the -th group, is the gas density, and is the path-length travelled by the CRs. For simplicity, we assume here a uniform gas density and neglect the dependence of the confinement time on the charge of the nuclei. We account for the helium contribution in the interstellar medium by means of an enhancement factor2, around , as explained in the following. Particle physics enters via the inelastic cross section for an interaction of a nucleus of mass number and energy per nucleon with a proton, and the neutrino production spectrum per inelastic event. The latter is defined as the convolution of the production spectra for different hadron species and the spectra for their decays into neutrinos,

 dnAp→ν(EA,E)dE=∑h∫dEhdnAp→h(EA,Eh)dEh ×dndech→ν(Eh,E)dE. (3)

Introducing the energy fraction of the produced neutrinos in Eq. (2), we can rewrite it in the case of power-law energy spectra of the CRs, , as

 Iν(E)=~εMdCRngas∑iIi(E)ZνAi(E,αi), (4)

where the so-called -factors3 (20) for neutrino production are defined as

 ZνA(E,α)=A−α∫10dzzα−1σinelAp(E/z) ×dnAp→ν(E/z,z)dz. (5)

Using (3), it is convenient to express via the -factors for hadron production as follows

 ZνA(E,α)=∑h∫10dzνzα−1νfdech→ν(zν)ZhA(E/zν,α), (6)

where are defined similarly to as

 ZhA(E,α)=A−α∫10dzzα−1σinelAp(E/z) ×dnAp→h(E/z,z)dz (7)

and are the (Lorentz-invariant) spectra for hadron decays into neutrinos,

 fdech→ν(z+ν)=dndech→ν(E,z+ν)dz+ν. (8)

In deriving (6), we used the high-energy limit

 z+ν≡(Eν+pzν)/(Eh+pzh)≃Eν/Eh=zν. (9)

Thus, all the dependence on the properties of proton-proton and nucleus-proton interactions is contained in the hadronic -factors . They also define how strongly the contribution to neutrino production from primary CR nuclei is suppressed relative to the one of protons. Due to the steep slopes of the primary spectra, these -factors are dominated by the hadron spectra in the very forward direction. This allows one to estimate the suppression of the nuclear contribution, using the well-known relation for the mean number of interacting (“wounded”) projectile nucleons for nucleus – nucleus collisions (21),

 ⟨nwpAB(E)⟩=AσinelpB(E)σinelAB(E). (10)

This relation holds both in the Glauber approach and in Reggeon Field Theory, if one neglects the contribution of target diffraction (22).

Thus, for the forward () spectra of secondary hadrons we get

 dnAp→h(E,z)dz≃⟨nwpAp(E)⟩dnpp→h(E,z)dz =Aσinelpp(E)σinelAp(E)dnpp→h(E,z)dz. (11)

Substituting (11) into (7), we obtain

 ZhA(E,α)≃A1−αZhp(E,α), (12)

 Iν(E)≃~εMdCRngas∑iIi(E)Zνp(E,αi)A1−αii. (13)

Thus the simple rule which is often applied to convert neutrino fluxes from collisions into those of collisions is only modified by the ratio of -factors , if the spectral slopes for CR nuclei differ from the one for protons, .

In Table 1, we present the -factors multiplied by for charged pion and kaon production4 by different primary nuclei, calculated with the QGSJET-II-04 model for the primary CR slope . Additionally, we demonstrate in Table 2 the dependence of these factors on the slopes of the primary spectra for GeV. As one can see from the Tables, Eq. (12) holds here to a very good accuracy: The factors for different primary particles agree with each other to better than 10% accuracy. The relatively weak energy-dependence of these factors stems from the energy rise of the inelastic nucleus-proton cross sections.

We also use the ratios of -factors for hadron production on helium and proton targets , as compiled in Table 3 for GeV, to determine the enhancement factor which accounts for the contribution to neutrino production from CR interactions with the helium component of the ISM. Here is defined by Eq. (7) with the replacements and . For , we have for all CR primaries; typical deviations do not exceed the 10% level5 [c.f. Eq. (12) and Table 3]. For the He/H abundance ratio in the ISM (23), we thus obtain

 ~εM=1+RHe/HZhA|HeZhA≃1.3. (14)

To check the model dependence of our results, we repeated the same calculations as in Table 1 for and GeV, using the EPOS-LHC model (24); (25), the results being collected in Table 4.

For EPOS-LHC, the deviations from Eq. (12) are more significant, reaching 20% in the case of primary iron. More importantly, also the predicitions for the -factors for primary protons by EPOS-LHC are higher than by QGSJET-II-04. Thus this range defines the characteristic uncertainty of our results that arises from the treatment of hadronic interactions.

## Iii Neutrino flux in the knee region

The elemental composition of the CR flux below  eV is relatively well determined and can be be described to first order by power-laws. At higher energies, the low CR flux prevents direct measurements, and the elemental composition of the CR flux becomes rather uncertain; for a review of experimental methods and results see e.g. Ref. (26). While there exist yet substantial uncertainties concerning the partial contributions of different mass groups to the primary CR composition, there is a general agreement that the knee in the total CR spectrum at  PeV coincides with a suppression of the primary proton flux, and that the composition becomes increasingly heavier in the energy range between the knee and  eV (12); (15); (14); (27).

Explanations for the origin of the knee fall in three main categories. First, there have been speculations that interactions may change in the multi-TeV region and the CR flux may be suppressed because of additional energy loss channels. This possibility is now excluded by LHC data (28). Second, the knee may correspond to the maximum rigidity to which CRs can be accelerated by the dominant population of Galactic CR sources (29); (30). Third, the knee energy may correspond to the rigidity at which the CR Larmor radius is of the order of the coherence length of the turbulent magnetic field in the Galactic disk. As a result, a transition from large-angle to small-angle scattering or Hall diffusion is expected, the energy dependence of the confinement time changes which in turn induces a steepening of the CR spectrum (31); Candia et al. (2002); Candia et al. (2003); (34).

Both the second and third possibilities lead to a rigidity-dependent sequence of knees at , a behavior first suggested by Peters Peters (1961). In contrast to models in category 2, those of category 3 predict both the position of the knee and the rigidity dependent suppression of the different CR components for a given model of the Galactic magnetic field (34).

Various models which describe the elemental composition of the total CR flux have been developed. The parametrization of Honda and Gaisser (36) is widely used in calculations of the atmospheric neutrino flux up to energies  GeV. Since it does not attempt to include the CR knee, it cannot be extrapolated into the energy range of our interest. The poly-gonato model (37) is a fit of rigidity-dependent knees at to measurements of the total CR intensity. Below and above , the fluxes of individual CR nuclei are assumed to follow power-laws and which are smoothly interpolated. The fluxes are assumed to steepen by a common amount, with , for more details see (37).

The CR intensity in this model, split into five elemental groups, is shown in Fig. 1. We neglected elements with which give in the poly-gonato model (37) an important contribution to the Galactic CR intensity at the highest energies, because their contribution to the neutrino intensity is—as expected from Eq. (12) and as we will see in the following—negligible. The steepening of the CR intensity around the knee is very pronounced in this model, , and as a result the composition is heavier than suggested by KASCADE-Grande data. The resulting neutrino intensity is shown in Fig. 2, where we assumed that the CR nuclei cross the grammage g/cm; this value corresponds for primary protons with energy 1 PeV to the interaction depth . Above few hundred TeV, the intensity is dominated by the proton contribution and is strongly decreasing with energy as . Therefore the expected neutrino energy distribution disagrees with the neutrino spectrum suggested by the IceCube excess, in particular with the two PeV events.

The Hillas model (30) and its variants belong to the category 2, associating the knee with the maximal rigidity achievable in the dominant population of Galactic CR sources. Moreover, the Hillas model assumes that the ankle signals the transition from Galactic to extragalactic CRs. Therefore, an additional population of Galactic CR sources must exist (“the component B” of Ref. (30)) which fills the gap between the knee and the ankle. Thus the Hillas model contains two Galactic components. Each population is assumed to contain five elemental groups and cuts off at a characteristic rigidity.

Variations of the original Hillas model were presented in Refs. (38); (39). We use here the new parametrisation H3a given in Tab. 3 of Ref. (39), where the first population is fitted only above  GeV, i.e. above the hardening observed by CREAM (13) and PAMELA (40). The intensity of Galactic CRs in this model is shown in Fig. 3, the resulting neutrino intensity in Fig. 4. The neutrino intensity is again dominated by the proton contribution; its shape is very similar to the neutrino intensity of the polygonato model.

Finally, we consider a parametrisation of the CR flux motivated by the recent results of Ref. (34): There, the escape of CRs from our Galaxy was studied calculating trajectories of individual CRs in models of the regular and turbulent Galactic magnetic field. For a coherence length  pc of the turbulent field and a reduced turbulent magnetic field, a knee-like structure at  eV was found, which is sufficiently strong to explain the proton knee observed by KASCADE. The resulting intensity of four other elemental groups are shown in Fig. 5. They are consistent with the energy spectra of CR nuclei determined by KASCADE and KASCADE-Grande. The resulting neutrino intensity for g/cm is shown in Fig. 6. The suppression of the neutrino intensity above the knee is less pronounced as in the previous models, since the decrease of the CR escape time slows down around  eV for a weak turbulent field. In the energy range between  eV and  eV, the neutrino intensity scales as .

In summary, we found that the neutrino intensity below  eV reflects the slope of protons and agrees therefore in all three models. In contrast, the exact position of the “neutrino knee” and the slope of neutrino intensity above this break depends on the nuclear composition and is therefore model dependent: A comparison of the neutrino intensity in the three models is shown in Fig. 7.

Finally, one may yet ask which of the three composition models considered is the more realistic one in the light of recent CR data. Comparing e.g. the CR spectra predicted by the Polygonato model to the intensities of individual groups of CR nuclei up to  eV, measured by the KASCADE and KASCADE-Grande experiments (12); (15), already an inspection by eye indicates that this model predicts a too heavy composition above the knee. The Hillas model of Ref. (39) describes well the average composition (represented e.g. by ) but fails to reproduce the up-turn of the light component around eV observed by KASCADE-Grande (15); (41). By contrast, such an up-turn around  eV is the characteristique feature of the escape model (34). As a result, the intensity of individual groups of CR nuclei measured by KASCADE and KASCADE-Grande (12); (15) is well reproduced in this model (34).

## Iv Comments on the IceCube events

We discuss now briefly how our results impact the Galactic interpretation of the IceCube excess.

### iv.1 Neutrinos from Galactic Sea CRs

The results of the previous section can be converted to the neutrino intensity resulting from collisions of Galactic sea CRs on gas after evaluating the interaction depth

 τν=σppinel(E)∫l.o.s.dsn(ρ,z), (15)

where is the distance from the Galactic center (GC) in the Galactic plane, the distance above the plane, and the distance from the Sun along the chosen line-of-sight (l.o.s.). We neglect for our estimations the contribution from primary nuclei, because we have seen that the neutrino flux is dominated by interactions.

The gas distribution in the Galactic disk can be modeled as with  cm at (increasing to  cm at the GC) and  kpc. Integrating (15) with  mb results in the maximal interaction depth of towards the GC.

At the reference energy  PeV, all three parametrisations predict a neutrino intensity around  GeV m s sr, corresponding to . Thus even in the direction of the largest expected intensity, the predicted neutrino intensity due to diffuse Galactic CR interactions is about two orders of magnitude too small compared to the IceCube excess. Moreover, the neutrino events should be concentrated within  (43), reflecting the very slim Galactic plane, which is much narrower than the latitude distribution of the IceCube events.

### iv.2 Neutrinos from Galactic CR sources

We consider next the neutrino flux produced close to recent CR sources. The propagation of CRs on distances can be approximated by diffusion (44). For  pc as found in Ref. (45) for the Galactic disk, the diffusion approach is marginally justified for the time scales, to  yr, we consider.

Galactic accelerators able to produce CRs with energies  eV have typically only short life-times: For instance, the highest energy particles produced by a supernova remnant (SNR) are thought to escape at the end of the Sedov phase after few 100 yrs. Approximating therefore the accelerator as a bursting source, the number density of CRs at the distance is given by

 nCR(E,r)=Q(E)π3/2r3diffexp[−r2/r2diff], (16)

with assuming no energy losses. We use as diffusion coefficient  cm/s at our reference energy  PeV and assume that the source injects instantaneously CRs with the total energy  erg in protons with an injection spectrum between the minimal energy  GeV and a maximal energy  PeV. We choose suggested by shock acceleration. Then PeV CRs are concentrated within  pc after 3000 yr.

The brightest spots in the Galactic neutrino sky are likely giant molecular clouds (GMC) immersed into the CR overdensities close to recent CR sources. The neutrino flux from a point source at the distance is given by

 ϕν(E)=~εMcσinel4πd2MclmpnCR(E)Yν(E), (17)

where denotes the neutrino yield . Assuming as cloud mass and as distance  kpc results in the neutrino flux

 E2ϕν(E)≃140eVcm−2sr−1. (18)

Sources of this kind would be clearly visible on the neutrino sky as seen by IceCube. Choosing as source rate , which coincides with the Galactic SN rate, implies that the average number of such sources present in the Galaxy equals . Using as volume of the Galactic disk  kpc, there are on average 0.5 sources within one kpc distance to an observer. We note also that the presence of GMCs close to SNRs is not unnatural, since they are born most likely in OB associations.

Any source of high-energy neutrinos produces also -rays. In Ref. (7), the 1–10 TeV -ray flux from known sources in the direction towards the GC was compared with the IceCube excess. Extrapolating the -ray flux of these sources to higher energies implies a neutrino flux which is an order of magnitude smaller than the one required to explain the IceCube excess. This discrepancy could be explained by the slower diffusion of low-energy CRs which have not yet reached the GMC.

Limits on the fraction of photons in the CR flux as e.g. those of CASA-MIA, KASCADE and IceCube can be used to constrain Galactic neutrino sources (5). However, PeV -rays can be absorbed both in sources and during propagation by pair production on star light and CMB photons. Moreover, these gamma-ray limits are biased towards the northern hemisphere. As a result, they do not exclude the case that (a fraction of) the IceCube excess has a Galactic origin.

## V Summary

We have calculated the neutrino yield from collisions of CR nuclei on gas using the event generator QGSJET-II. Our numerical results assuming power-law fluxes for the primary CR nuclei can be used in all applications where the neutrino yield is dominated by the contribution of heavy nuclei. In the case of Galactic CRs, we found that the neutrino flux is well approximated by accounting only for the proton component in the CR flux. Since the proton intensity above  eV varies considerably in different parametrisations of the elemental composition of Galactic CRs, the resulting variations in the predicted neutrino intensity offer the possibility to distinguish between these options. The helium contribution in the interstellar medium to the neutrino flux can be accounted for approximately by employing an enhancement factor, around .

We have compared the spectral shape and the magnitude of the predicted neutrino intensity to the IceCube excess. The slope of the neutrino intensity from interactions of sea CRs is close to both in the poly-gonato and the Hillas model at high energies, while it reflects the proton slope at low energies. The break is less pronounced in the escape model, where the neutrino intensity scales as in the knee region. Thus the expected slopes in the poly-gonato and the Hillas model are compatible with a cutoff (or break) in the energy spectrum of the IceCube events. However, the energy scale of the break, the total number of events expected and their arrival directions make an explanation of these events by Galactic sea CRs very unlikely. By contrast, the energy spectrum of neutrinos produced by proton–gas interactions close to sources has the same slope as the CR injection spectrum. Since the injection spectrum is flatter, this results in a better agreement with the energy spectrum of the IceCube excess.

The required magnitude of the neutrino flux can be achieved, if the sources illuminate a close GMC. Moreover, the sources have to be relatively nearby,  kpc. A small distance to the sources would also explain, why the IceCube neutrino events are not as concentrated towards the Galactic plane as expected. In conclusion, we consider it as a viable option that the IceCube excess is partly connected to nearby Galactic CR sources.

Finally, let us comment on the position of the neutrino knee: While its exact shape is model dependent, its position at  eV reflects simply the fact that the maximal neutrino energy in interactions is approximately 10% of the energy of the proton primary. Our results for the Milky Way can be applied also to the neutrino energy spectrum expected from other normal galaxies, in particular from starburst galaxies. These galaxies have magnetic fields which are two orders of magnitude higher than the Galactic magnetic field (46). As a consequence, in the model of Ref. (34) where the knee is caused by the escape of CRs, its position and thus the one of the neutrino knee is shifted by up too two orders of magnitude. In contrast, such a shift is not expected in models of the Hillas type, because there the maximal rigidity to which the dominant population of supernovae can accelerate CRs is determined not by the ambient magnetic field but by the magnetic field created by CR instabilities.

###### Acknowledgements.
We would like to thank Tom Gaisser for communications about the Hillas model. We are grateful to Gwenael Giacinti and Dima Semikoz for discussions and collaboration on CR propagation on which the “escape model” from Ref. (34); (42) is based. MK would like to thank Nordita for hospitality and support during the program “News in Neutrino Physics.” SO acknowledges support from NASA through the grants NNX13AC47G and NNX13A092G.

### Footnotes

1. While the cascade bound is derived considering CR interactions with the extragalactic -ray background during propagation, a similiar bound applies to many neutrino sources themselves.
2. Note that in contrast to its usual definition, see e.g. (19), accounts only for nuclei in the interstellar medium.
3. Originally, -moments were introduced calculating inclusive fluxes of high energy muons and neutrinos resulting from CR interactions in the atmosphere, without the factor .
4. For the -factors of kaons, the relation holds.
5. Since both and are defined for the same spectra of primary mass groups, the effects of primary abundances and of the energy-dependence of these -factors are largely cancelled in their ratios in the expression (14) for the enhancement factor .

### References

1. M. G. Aartsen et al. [IceCube Collaboration], Science 342, 1242856 (2013) [arXiv:1311.5238 [astro-ph.HE]].
2. M. G. Aartsen et al. [IceCube Collaboration], Phys. Rev. Lett. 111, 021103 (2013) [arXiv:1304.5356 [astro-ph.HE]].
3. L. A. Anchordoqui, H. Goldberg, M. H. Lynch, A. V. Olinto, T. C. Paul and T. J. Weiler, arXiv:1306.5021 [astro-ph.HE].
4. I. Cholis and D. Hooper, JCAP 06, 030 (2013) [arXiv:1211.1974 [astro-ph.HE]]; M. D. Kistler, T. Stanev and H. Yuksel, [arXiv:1301.1703 [astro-ph.HE]]; F. W. Stecker, Phys. Rev. D 88, 047301 (2013) [arXiv:1305.7404 [astro-ph.HE]]; K. Murase and K. Ioka, Phys. Rev. Lett. 111, 121102 (2013) [arXiv:1306.2274 [astro-ph.HE]]; W. Winter, Phys. Rev. D 88, 083007 (2013) [arXiv:1307.2793 [astro-ph.HE]]; J. C. Joshi, W. Winter and N. Gupta, arXiv:1310.5123 [astro-ph.HE].
5. M. Ahlers and K. Murase, arXiv:1309.4077 [astro-ph.HE].
6. D. B. Fox, K. Kashiyama and P. Meszaros, Astrophys. J. 774, 74 (2013) [arXiv:1305.6606 [astro-ph.HE]].
7. A. Neronov, D. V. Semikoz and C. Tchernin, [arXiv:1307.2158 [astro-ph.HE]]; J. C. Joshi, W. Winter and N. Gupta, arXiv:1310.5123 [astro-ph.HE].
8. E. Roulet, G. Sigl, A. van Vliet and S. Mollerach, JCAP 1301, 028 (2013) [arXiv:1209.4033 [astro-ph.HE]].
9. V. Berezinsky, A. Gazizov, M. Kachelrieß and S. Ostapchenko, Phys. Lett. B 695, 13 (2011); [arXiv:1003.1496 [astro-ph.HE]].
10. A. Esmaili and P. D. Serpico, JCAP 1311, 054 (2013) [arXiv:1308.1105 [hep-ph]] and references therein.
11. R. Silberberg and M.M. Shapiro, Proc. 15th ICRC, v. 6, p. 237 (1977); F. W. Stecker, Astrophys. J. 228, 919 (1979).
12. T. Antoni et al. [KASCADE Collaboration], Astropart. Phys. 24, 1 (2005) [arXiv:astro-ph/0505413].
13. H. S. Ahn et al., Astrophys. J. 714, L89 (2010) [arXiv:1004.1123 [astro-ph.HE]].
14. R. Abbasi et al. [IceCube Collaboration], Astropart. Phys. 42, 15 (2013) [arXiv:1207.3455 [astro-ph.HE]].
15. W. D. Apel et al. [KASCADE-Grande Collaboration], Astropart. Phys. 47, 54 (2013) [arXiv:1306.6283 [astro-ph.HE]].
16. For an earlier study discussing the connection of a rigidity dependent knee and neutrino fluxes see J. Candia and E. Roulet, JCAP 0309, 005 (2003).
17. S. Ostapchenko, Phys. Rev. D 83, 014018 (2011).
18. S. Ostapchenko, Progr. Theor. Phys. Suppl. 193, 204 (2012); EPJ Web Conf. 52 (2013) 02001.
19. M. Kachelrieß, I. V. Moskalenko and S. S. Ostapchenko, Astrophys. J. 789, 166 (2014).
20. T. K. Gaisser, Cosmic Rays and Particle Physics (Cambridge University Press, Cambridge, England 1990).
21. A. Białas, M. Bleszynski and W. Czyz, Nucl. Phys. B 111, 461 (1976).
22. N. N. Kalmykov and S. S. Ostapchenko, Phys. At. Nucl. 56, 346 (1993).
23. J. P. Meyer, Astrophys. J. Suppl. 57, 173 (1985).
24. K. Werner, F.-M. Liu and T. Pierog, Phys. Rev. C 74, 044902 (2006).
25. T. Pierog, Yu. Karpenko, J. M. Katzy, E. Yatsenko and K. Werner, arXiv:1306.0121 [hep-ph].
26. J. Blümer, R. Engel and J. R. Hörandel, Prog. Part. Nucl. Phys. 63, 293 (2009) [arXiv:0904.0725 [astro-ph.HE]].
27. M. Aglietta, B. Alessandro, P. Antonioli et al. [EAS-TOP Collaboration], Astropart. Phys. 21, 583 (2004).
28. D. d’Enterria, R. Engel, T. Pierog, S. Ostapchenko and K. Werner, Astropart. Phys. 35, 98 (2011).
29. T. Stanev, P. L. Biermann and T. K. Gaisser, Astron. Astrophys. 274, 902 (1993); K. Kobayakawa, Y. Sato and T. Samura, Phys. Rev. D 66, 083004 (2002).
30. A. M. Hillas, J. Phys. G 31, R95 (2005).
31. Ptuskin, V. S. et al.., Astron. Astrophys. 268, 726 (1993).
32. J. Candia, E. Roulet and L. N. Epele, JHEP 0212, 033 (2002).
33. J. Candia, S. Mollerach, and E. Roulet, JCAP 0305, 003 (2003).
34. G. Giacinti, M. Kachelrieß and D. V. Semikoz, Phys. Rev. D 90, 041302 (2014).
35. B. Peters, Nuovo Cim. 22, 800 (1961); see also G. T. Zatsepin, N. N. Gorunov and L. G. Dedenko, lzv. Akad. Nauk USSR Ser. Fiz. 26 (1962) 685.
36. T. K. Gaisser and M. Honda, Ann. Rev. Nucl. Part. Sci. 52, 153 (2002) [hep-ph/0203272]; for an update see T. K. Gaisser, arXiv:1303.1431 [hep-ph].
37. J. R. Hörandel, Astropart. Phys. 19, 193 (2003) [astro-ph/0210453].
38. T. K. Gaisser, Astropart. Phys. 35, 801 (2012).
39. T. K. Gaisser, T. Stanev and S. Tilav, Front. Phys. China 8, 748 (2013) [arXiv:1303.3565 [astro-ph.HE]].
40. O. Adriani et al. [PAMELA Collaboration], Science 332, 69 (2011) [arXiv:1103.4055 [astro-ph.HE]].
41. W. D. Apel et al. [KASCADE-Grande Collaboration], Phys. Rev. D 87, 081101 (2013).
42. G. Giacinti, M. Kachelrieß and D. V. Semikoz, in preparation.
43. V. S. Berezinsky, T. K. Gaisser, F. Halzen and T. Stanev, Astropart. Phys. 1, 281 (1993); C. Evoli, D. Grasso and L. Maccione, JCAP 0706, 003 (2007).
44. G. Giacinti, M. Kachelrieß and D. V. Semikoz, Phys. Rev. Lett. 108, 261101 (2012); Phys. Rev. D 88, 023010 (2013) [arXiv:1306.3209 [astro-ph.HE]].
45. M. Iacobelli et al., Astron. Astrophys. 558, A72 (2013).
46. B. C. Lacki and R. Beck, Mon. Not. R. Astron. Soc. 430, 3171 (2013) [arXiv:1301.5391 [astro-ph.CO]].
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters