Cosmic-ray energy spectrum and composition up to the ankle - the case for a second Galactic component

# Cosmic-ray energy spectrum and composition up to the ankle − the case for a second Galactic component

S. Thoudam Department of Astrophysics/IMAPP, Radboud University, P.O. Box 9010, 6500 GL Nijmegen, The Netherlands Now at: Department of Physics and Electrical Engineering, Linnéuniversitetet, 35195 Växjö, Sweden E-mail: satyendra.thoudam@lnu.se    J.P. Rachen Department of Astrophysics/IMAPP, Radboud University, P.O. Box 9010, 6500 GL Nijmegen, The Netherlands    A. van Vliet Department of Astrophysics/IMAPP, Radboud University, P.O. Box 9010, 6500 GL Nijmegen, The Netherlands    A. Achterberg Department of Astrophysics/IMAPP, Radboud University, P.O. Box 9010, 6500 GL Nijmegen, The Netherlands    S. Buitink Astronomical Institute, Vrije Universiteit Brussel, Pleinlaan 2, 1050 Brussels, Belgium    H. Falcke Department of Astrophysics/IMAPP, Radboud University, P.O. Box 9010, 6500 GL Nijmegen, The Netherlands NIKHEF, Science Park Amsterdam, 1098 XG Amsterdam, The Netherlands ASTRON, Postbus 2, 7990 AA Dwingeloo, The Netherlands    J.R. Hörandel Department of Astrophysics/IMAPP, Radboud University, P.O. Box 9010, 6500 GL Nijmegen, The Netherlands NIKHEF, Science Park Amsterdam, 1098 XG Amsterdam, The Netherlands
August 4, 2019
###### Key Words.:
Galaxy — cosmic rays — diffusion — ISM: supernova remnants — Stars: winds — Stars: Wolf-Rayet

Motivated by the recent high-precision measurements of cosmic rays by several new-generation experiments, we have carried out a detailed study to understand the observed energy spectrum and composition of cosmic rays with energies up to about  eV. Our study shows that a single Galactic component with subsequent energy cut-offs in the individual spectra of different elements, optimised to explain the observed elemental spectra below  eV and the ‘knee’ in the all-particle spectrum, cannot explain the observed all-particle spectrum above  eV. We discuss two approaches for a second component of Galactic cosmic rays – re-acceleration at a Galactic wind termination shock, and supernova explosions of Wolf-Rayet stars, and show that the latter scenario can explain almost all observed features in the all-particle spectrum and the composition up to  eV, when combined with a canonical extra-galactic spectrum expected from strong radio galaxies or a source population with similar cosmological evolution. In this two-component Galactic model, the knee at  eV and the ‘second knee’ at  eV in the all-particle spectrum are due to the cut-offs in the first and second components, respectively. We also discuss several variations of the extra-galactic component, from a minimal contribution to scenarios with a significant component below the ‘ankle’ (at  eV), and find that extra-galactic contributions in excess of regular source evolution are neither indicated nor in conflict with the existing data. We also provide arguments that an extra-galactic contribution is unlikely to dominate at or below the second knee. Our main result is that the second Galactic component predicts a composition of Galactic cosmic rays at and above the second knee that largely consists of helium or a mixture of helium and CNO nuclei, with a weak or essentially vanishing iron fraction, in contrast to most common assumptions. This prediction is in agreement with new measurements from LOFAR and the Pierre Auger Observatory which indicate a strong light component and a rather low iron fraction between and  eV.

## 1 Introduction

Until a decade ago, the cosmic ray spectrum from  GeV to  GeV was seen as a power law with two main features: a steepening from a spectral index to at about  GeV, commonly called the ‘knee’, and a flattening back to at about  GeV, consequently denoted as the ‘ankle’. Phenomenological explanations for the knee have been given due to propagation effects in the Galaxy (Ptuskin et al., 1993), progressive cutoffs in the spectra of nuclear components from hydrogen to lead (Hörandel, 2003a), or re-acceleration at shocks in a Galactic wind (Völk & Zirakashvili, 2004), but left open the question of the primary Galactic accelerators producing these particles. Explanations based on source physics have been mostly built on the assumption that supernova remnants, on grounds of energetics known as one of the most promising sources for cosmic rays (Baade & Zwicky, 1934), accelerate cosmic rays at shocks ploughing into the interstellar medium to energies up to about  GeV (Lagage & Cesarsky, 1983; Axford, 1994). This may extend to  GeV if they are propagating in fast and highly magnetised stellar winds (Völk & Biermann, 1988; Biermann & Cassinelli, 1993), or if non-linear effects in the acceleration process are considered (Bell & Lucek, 2001). The combination of such components could eventually explain cosmic rays below and above the knee as a superposition of components of different nuclei, as shown, for example by Stanev et al. (1993). At energies above   GeV this steep component was assumed to merge into a flatter extra-galactic component (Rachen et al., 1993; Berezinsky et al., 2004), explaining the ankle in the spectrum. For this extra-galactic component, sources on all scales have been proposed: From clusters of galaxies (Kang et al., 1996) through radio galaxies (Rachen & Biermann, 1993), compact AGN jets (Mannheim et al., 2001) to gamma-ray bursts (Waxman, 1995). It was commonly assumed to be dominated by protons. Eventually, at  GeV the cosmic ray spectrum was believed to terminate in the so-called GZK cutoff (Greisen, 1966; Zatsepin & Kuzmin, 1996) due to interaction with cosmic microwave background (CMB) photons.

The new data have led to a number of theoretical modifications of the standard model. The spectral hardening at TeV energies has been explained as due to the hardening in the source spectrum of cosmic rays (Biermann et al., 2010a; Ohira et al., 2011; Yuan et al., 2011; Ptuskin et al., 2013), as a propagation effect (Tomassetti, 2012; Blasi et al., 2012), the effect of re-acceleration by weak shocks (Thoudam & Hörandel, 2014) or the effect of nearby sources (Thoudam & Hörandel, 2012, 2013; Erlykin & Wolfendale, 2012). At high energies, the increasing mean mass around the knee still fits well the idea of progressive cut-offs (Hörandel, 2003a), if the nuclear species are constrained to masses up to iron and thus limited to energies below about  GeV. The light composition around the ankle revived interest in the so-called ‘proton dip model’, which explains the ankle feature as due to an extra-galactic propagation effect of protons producing electron-positron pairs at the CMB (Berezinsky & Grigorieva, 1988; Berezinsky et al., 2006). This would imply that the cosmic ray spectrum below the ankle is, at least in part, of extra-galactic origin. While the recent measurement of proton fraction at the ankle by the Pierre Auger Collaboration (Aab et al., 2014) has raised problems with this approach, as the model is compatible only with more than protons (Berezinsky et al., 2006), a number of new models have been suggested, involving compact sources with significant photo-disintegration of nuclei during acceleration (Globus et al., 2015a; Unger et al., 2015), or as a component with primordial element composition accelerated at clusters of galaxies and limited by pair production losses in the CMB (Rachen, 2016). However, with all these new ideas, big questions remain open: How does the cosmic ray component at the knee connect to the one at the second knee to ankle regime, and where is the transition from Galactic to extra-galactic cosmic rays?

In this work, we revisit the basic models of Galactic cosmic-ray production in view of the currently available data. We start by developing a detailed model description for low-energy cosmic rays assuming them to be primarily produced inside supernova remnants (SNRs) present in the interstellar medium (hereafter, these cosmic rays will be referred to as the ‘SNR-CRs’). This model, described in Section 2, has been demonstrated to explain the observed spectral hardening of protons and helium nuclei in the TeV region and, at the same time, explains the observed composition of cosmic rays at low energies (Thoudam & Hörandel, 2014). The model prediction will be extended to high energies, and compared with the observed all-particle energy spectrum. It will be shown that SNR-CRs cannot explain the observed energy spectrum above  GeV. We then revisit two possibilities for a second Galactic component in Section 3: (a) The re-acceleration of SNR-CRs escaped into the Galactic halo by the Galactic wind termination shocks (Jokipii & Morfill, 1987; Zirakashvili & Völk, 2006), and (b) the contribution of cosmic rays from the explosions of Wolf-Rayet stars in the Galaxy (Biermann & Cassinelli, 1993). The possibility of a second Galactic component has also been discussed in Hillas (2005) who considered Type II SNRs expanding into a dense wind of their precursor stars. For both the scenarios considered in the present work, we assume the extra-galactic proton component used by Rachen et al. (1993) to obtain proper results for total spectrum and composition at energies just below the ankle in Section 4. In Section 5 we then check the effect of other hypotheses for the extra-galactic component, using (1) a phenomenological ‘minimal model’ derived from composition results measured at the Pierre Auger Observatory (di Matteo et al., 2015), (2) the minimal model plus the ‘primordial cluster component’ introduced by Rachen (2016), and (3) the ‘extra-galactic ankle’ model by Unger et al. (2015). In Section 6, we present a discussion of our results and their implications, and other views on the cosmic rays below  GeV, followed by our conclusions in Section 7.

## 2 Cosmic rays from supernova remnants (SNR-CRs)

Although the exact nature of cosmic-ray sources in the Galaxy is not yet firmly established, supernova remnants are considered to be the most plausible candidates both from the theoretical and the observational points of view. It has been theoretically established that shock waves associated with supernova remnants can accelerate particles from the thermal pool to a non-thermal distribution of energetic particles. The underlying acceleration process, commonly referred to as the diffusive shock acceleration process, has been studied quite extensively, and it produces a power-law spectrum of particles with a spectral index close to (Krymskii, 1977; Bell, 1978; Blandford & Ostriker, 1978; Drury, 1983; Ptuskin et al., 2010; Caprioli et al., 2011), which is in good agreement with the values inferred from radio observation of supernova remnants (Green, 2009). Moreover, the total power of ergs s injected by supernova explosions into the Galaxy, considering a supernova explosion energy of ergs and an explosion frequency of yr, is more than sufficient to maintain the cosmic-ray energy content of the Galaxy. In addition to the radio measurements, observational evidence for the presence of high-energy particles inside supernova remnants is provided by the detection of non-thermal X-rays (Vink & Laming, 2003; Parizot et al., 2006) and TeV gamma rays from a number of supernova remnants (Aharonian et al., 2006, 2008; Albert et al., 2007). For instance, the detection of TeV gamma rays up to energies close to  TeV from the supernova remnant RX J1713.7-3946 by the H.E.S.S. Cherenkov telescope array indicates that particles with energies up to PeV can be accelerated inside supernova remnants (Aharonian et al., 2007).

### 2.1 Transport of SNR-CRs in the Galaxy

After acceleration by strong supernova remnant shock waves, cosmic rays escape from the remnants and undergo diffusive propagation through the Galaxy. During the propagation, some fraction of cosmic rays may further get re-accelerated due to repeated encounters with expanding supernova remnant shock waves in the interstellar medium (Wandel, 1988; Berezhko et al., 2003). This re-acceleration is expected to be produced mainly by older remnants, with weaker shocks, because of their bigger sizes. Therefore, the re-acceleration is expected to generate a particle spectrum which is steeper than the initial source spectrum of cosmic rays produced by strong shocks. This model has been described in detail in Thoudam & Hörandel (2014), and it has been shown that the re-accelerated cosmic rays can dominate the GeV energy region while the non-re-accelerated cosmic rays dominate at TeV energies, thereby explaining the observed spectral hardening in the TeV region. Below, we briefly summarise some key features of the model which are important for the present study.

The steady-state transport equation for cosmic-ray nuclei in the Galaxy in the re-acceleration model is described by,

 ∇⋅(D ∇N)−[¯nvσ+ξ]δ(z)N +[ξsp−s∫pp0duN(u)us−1]δ(z)=−Qδ(z), (1)

where we have adopted a cylindrical geometry for the propagation region described by the radial and vertical coordinates with representing the Galactic plane. We assume the region to have a constant halo boundary at , and no boundary in the radial direction. This is a reasonable assumption for cosmic rays at the galacto-centric radius of the Sun as the majority of them are produced within a radial distance from the Sun (Thoudam, 2008). Choosing a different (smaller) halo height for the Galactic centre region, as indicated by the observed WMAP haze (Biermann et al., 2010b), will not produce significant effects in our present study. represents the differential number density of the cosmic-ray nuclei with momentum/nucleon , and is the injection rate of cosmic rays per unit volume by supernova remnants in the Galaxy. The diffusive nature of the propagation is represented by the first term in Equation 2.1. The diffusion coefficient is assumed to be a function of the particle rigidity as, , where is the diffusion constant, with and representing the velocity of the particle and the velocity of light respectively,  GV is a constant, and is the diffusion index. The rigidity is defined as , where and represent the mass number and the charge number of the nuclei respectively, and is the charge of an electron. The second term in Equation 2.1 represents the loss of particles during the propagation due to inelastic interaction with the interstellar matter, and also due to re-acceleration to higher energies, where represents the surface density of matter in the Galactic disk, is the inelastic interaction cross-section, and corresponds to the rate of re-acceleration. We take , where is the volume occupied by a supernova remnant of radius re-accelerating the cosmic rays, is a correction factor that is introduced to account for the actual unknown size of the remnants, and is the frequency of supernova explosions per unit surface area in the Galactic disk. The term containing the integral in Equation 2.1 represents the gain in the number of particles due to re-acceleration from lower energies. The effect of Galactic wind and ionisation losses which are important mostly at low energies, below  GeV/nucleon, are not included explicitly in the transport equation. Instead, we introduce a low-momentum cut-off, MeV/nucleon, in the particle distribution to account for the effect on the number of low-energy particles available for re-acceleration in the presence of these processes (Wandel et al., 1987). We assume that re-acceleration instantaneously produces a power-law spectrum of particles with spectral index . The source term can be expressed as , where for represents a Heaviside step function, and the source spectrum is assumed to follow a power-law in total momentum with an exponential cut-off which, in terms of momentum/nucleon, can be written as

 Q(p)=AQ0(Ap)−qexp(−ApZpc), (2)

where is a normalisation constant which is proportional to the amount of energy channelled into cosmic rays by a single supernova event, is the spectral index, and is the cut-off momentum for protons. The exponential cut-off in Equation 2 represents a good approximation for particles at the shock produced by the diffusive shock acceleration mechanism (see e.g. Malkov & Drury, 2001). We assume that the maximum energy for cosmic-ray nuclei produced by the supernova shock is times the maximum energy for protons. Based on the observed high concentration of supernova remnants and atomic and molecular hydrogen near the Galactic disk, in Equation 2.1, we assume that both cosmic-ray sources and interstellar matter are distributed in the disk (i.e. at ). The distributions are assumed to be uniform, and extended up to a radius .

Recalling the analytical solution of Equation 2.1 derived in Thoudam & Hörandel (2014), the cosmic-ray density at the position for follows,

 N (z,p)=¯νR∫∞0dksinh[k(L−z)]sinh(kL)×J1(kR)B(p){Q(p) +ξsp−s∫pp0dp′p′sQ(p′)A(p′)exp(ξs∫pp′A(u)du)}, (3)

where is a Bessel function of order 1, and the functions and are given by,

 B(p)=2D(p)kcoth(kL)+¯nv(p)σ(p)+ξ A(u)=1uB(u). (4)

From Equation 2.1, the cosmic-ray density at the Earth can be obtained by taking considering that our Solar system lies close to the Galactic plane.

### 2.2 Model prediction for the low-energy measurements

By comparing the abundance ratio of boron-to-carbon nuclei predicted by the model with the measurements, the cosmic-ray propagation parameters () and the re-acceleration parameters () have been obtained to be,  cm s, , , and (Thoudam & Hörandel, 2014). We adopt these values in our present study. The supernova remnant radius is taken to be  pc. The inelastic interaction cross-section for protons is taken from Kelner et al. (2006), and for heavier nuclei, the cross-sections are taken from Letaw et al. (1983). The surface matter density is taken as the averaged density in the Galactic disk within a radius equal to the size of the diffusion boundary . We choose  kpc, which gives an averaged surface density of atomic hydrogen of  atoms cm (Thoudam & Hörandel, 2013). An extra is further added to to account for the helium abundance in the interstellar medium. The radial extent of the source distribution is taken as  kpc. Each supernova explosion is assumed to release a total kinetic energy of  ergs, and the supernova explosion frequency is taken as  SNe Myr kpc. The latter corresponds to a rate of supernova explosions per century in the Galaxy.

Using the values of various parameters mentioned above, the energy spectra of SNR-CRs for different elements are calculated. In Figure 1, results for eight elements (proton, helium, carbon, oxygen, neon, magnesium, silicon and iron, which represent the dominant species at low energies) are compared with the measured data at low energies. The source parameters for the individual elements are kept free in the calculation, and they are optimised based on the observed individual spectra at low energies. The parameter values that best reproduce the measured data are listed in Table 1. The source spectral indices are in the range of , and out of the total of of the supernova explosion energy channelled into SNR-CRs, the largest fraction goes into protons at the level of , followed by helium nuclei with . The calculated spectra reproduce the measured data quite well including the behaviour of spectral hardening at TeV energies observed for protons and helium nuclei. In our model, the absence of such a spectral hardening for heavier nuclei is explained as due to the increasing effect of inelastic collision over re-acceleration with the increase in mass (Thoudam & Hörandel, 2014).

### 2.3 Extrapolation of the SNR-CR spectrum to high energies

In Figure 1, we also show an extrapolation of the model prediction to high energies. For protons, helium, carbon, silicon and iron nuclei, the predictions are compared with the available measurements from the KASCADE experiment above  GeV. The calculation assumes an exponential cut-off for the proton source spectrum at  GeV, and for the heavier nuclei at . This value of , which is obtained by comparing the predicted all-particle spectrum with the observed all-particle spectrum as shown in Figure 2, represents the maximum value permitted by the measurements. While obtaining the all-particle spectrum shown in Figure 2, we also include contributions from the sub-dominant primary cosmic-ray elements , calculated using elemental abundances at  GeV given in Hörandel (2003a) and a source index of 2.25. Their total contribution amounts up to of the all-particle spectrum. The predicted all-particle spectrum agrees with the data up to  GeV, and reproduces the observed knee at the right position. Choosing values larger than  GeV will produce an all-particle spectrum which is inconsistent both with the observed knee position and the intensity above the knee. Although our estimate for the best-fit value does not rely on the proton measurements at high energies, it can be noticed from Figure 1 that both the predicted proton and helium spectra are in good agreement (within systematic uncertainties) with the KASCADE data. For carbon, silicon and iron nuclei, the agreement with the data is less convincing, which may be related to the larger systematic uncertainties in the shapes of the measured spectra.

From Figure 2, it can be observed that, at energies around the knee, the all-particle spectrum is predicted to be dominated by helium nuclei, not by protons. The CREAM measurements have shown that helium nuclei become more abundant than protons at energies  GeV. Such a trend is also consistent with the KASCADE measurements above  GeV (see Figure 1). Based on our prediction, helium nuclei dominate the all-particle spectrum up to GeV, while above, iron nuclei dominate. The maximum energy of SNR-CRs, which corresponds to the fall-off energy of iron nuclei, is GeV. Although this energy is close to the position of the second knee, the predicted intensity is not enough to explain the observed intensity around the second knee. Our result shows that SNR-CRs alone cannot account for the observed cosmic rays above  GeV. At GeV, they contribute only of the observed data.

## 3 Additional component of Galactic cosmic rays

Despite numerous studies, it is not clearly understood at what energy the transition from Galactic to extra-galactic cosmic rays (EG-CRs) occurs. Although it was pointed out soon after the discovery of the CMB and the related GZK effect that it is possible to construct an all-extra-galactic spectrum of cosmic rays containing both the knee and the ankle as features of cosmological propagation (Hillas, 1967), the most natural explanation was assumed to be that the transition occurs at the ankle, where a steep Galactic component is taken over by a flatter extra-galactic one. To obtain a sharp feature like the ankle in such a construction, it is necessary to assume a cut-off in the Galactic component to occur immediately below it (Rachen et al., 1993; Axford, 1994), thus this scenario is naturally expecting a second knee feature. For a typical Galactic magnetic field strength of  G, the Larmor radii for cosmic rays of energy  GeV is  pc, much smaller than the size of the diffusion halo of the Galaxy, which is typically considered to be a few kpc in cosmic-ray propagation studies, keeping comic rays around the second knee well confined in the Galaxy. This suggests that the Galactic cut-off at this energy must be intrinsic to a source population or acceleration mechanism different from the standard supernova remnants we have discussed above. In an earlier work, Hillas (2005) considered an additional Galactic component resulting from Type II supernova remnants in the Galaxy expanding into a dense slow wind of the precursor stars. In the following, we discuss two other possible scenarios. The first is the re-acceleration of SNR-CRs by Galactic wind termination shocks in the Galactic halo (Jokipii & Morfill, 1987; Zirakashvili & Völk, 2006), and the second is the contribution of cosmic rays from the explosions of Wolf-Rayet stars in the Galaxy (Biermann & Cassinelli, 1993). Both these ideas have been explored in the past when detailed measurements of the cosmic-ray spectrum and composition at low and high energies were not available. Using new measurements of cosmic rays and astronomical data (like the Wolf-Rayet wind composition), our study can provide a more realistic estimate of the cosmic-ray contribution from these two possible mechanisms. In the following, the re-accelerated cosmic rays from Galactic wind termination shocks will be referred to as ‘GW-CRs’, and cosmic rays from Wolf-Rayet stars as ‘WR-CRs’. Some ramifications of these basic scenarios will be discussed in Section 6, after investigating the effect of different extra-galactic contributions below the ankle in Section 5.

### 3.1 Re-acceleration of SNR-CRs by Galactic wind termination shocks (GW-CRs)

The effect of Galactic winds on the transport of cosmic rays in the Galaxy has been discussed quite extensively (Lerche & Schlickeiser, 1982a; Bloemen et al., 1993; Strong & Moskalenko, 1998; Jones et al., 2001; Breitschwerdt et al., 2002). For cosmic rays produced by sources in the Galactic disk such as the SNR-CRs, the effect of winds on their transport is expected to be negligible above a few GeV as the transport is expected to be dominated mainly by the diffusion process. However, Galactic winds can lead to the production of an additional component of cosmic rays which can dominate at high energies. Galactic winds, which start at a typical velocity of about few km/s near the disk, reach supersonic speeds at distances of a few tens of kpc away from the disk. At about a hundred kpc distance or so, the wind flow terminates resulting into the formation of termination shocks. These shocks can catch the SNR-CRs escaping from the disk into the Galactic halo, and re-accelerate them via the diffusive shock acceleration process. The reaccelerated cosmic rays can return to the disk through diffusive propagation against the Galactic wind outflow. For an energy dependent diffusion process, only the high-energy particles may be effectively able to reach the disk.

To obtain the contribution of GW-CRs, we will first calculate the escape rate of SNR-CRs from the inner diffusion boundary, then propagate the escaped cosmic rays through the Galactic wind region, and calculate the cosmic-ray flux injected into the Galactic wind termination shocks. The escaped flux of SNR-CRs from the diffusion boundary, , can be calculated as,

 Fesc=[D∇N]z=±L=[DdNdz]z=±L, (5)

where is given by Equation 2.1. Equation 5 assumes that cosmic rays escape only through the diffusion boundaries located at . Under this assumption, the total escape rate of SNR-CRs is given by,

 Qesc=Fesc×2Aesc, (6)

where is the surface area of one side of the cylindrical diffusion boundary which is assumed to have the same radius as the Galactic disk, and the factor is to account for the two boundaries at . The propagation of the escaped SNR-CRs in the Galactic wind region is governed by the following transport equation:

 ∇.(Dw∇Nw−VNw)+∂∂p{∇.V3pNw}=−Qescδ(r), (7)

where we have assumed a spherically symmetric geometry characterised by the radial variable , represents the diffusion coefficient of cosmic rays in the wind region which is taken to be spatially constant, is the cosmic-ray number density, is the wind velocity which is assumed to increase linearly with and directed radially outwards, is a constant that denotes the velocity gradient, and is given by Equation 6. The exact nature of the Galactic wind is not known. The spatial dependence of the wind velocity considered here is based on the model of magnetohydrodynamic wind driven by cosmic rays, which shows that the wind velocity increases linearly with distance from the Galactic disk until it reaches an asymptotic value at a distance of around  kpc (Zirakashvili et al., 1996). The second term on the left-hand side of Equation 7 represents the loss of particles due to advection by the Galactic wind, and the third term represents momentum loss due to the adiabatic expansion of the wind flow which is assumed to be spherically symmetric. In writing Equation 7, considering that the size of the wind region is much larger than the size of the escaping region of the SNR-CRs, we neglect the size of the escaping region and consider to be a point source located at . By solving Equation 7 analytically, the density of cosmic rays at distance is given by (see Appendix A),

 Nw(r,p)=√~Vp28π3/2∫∞0 dp′Qesc(p′)[∫p′puDw(u)du]3/2 ×exp⎛⎜⎝−r2~Vp24∫p′puDw(u)du⎞⎟⎠. (8)

From Equation 3.1, the cosmic-ray flux with momentum/nucleon at the termination shock is obtained as,

 Fw(p)=[−Dw∂Nw∂r+VNw]r=Rsh, (9)

where represents the radius of the termination shock. The total rate of cosmic rays injected into the termination shock is given by,

 Qinj(p)=Fw(p)×Ash, (10)

where is the surface area of the termination shock. Assuming that only a certain fraction, , participates in the re-acceleration process, the cosmic-ray spectrum produced by the termination shock under the test particle approximation can be written as (Drury, 1983),

 Qsh(p)=γp−γexp(−ApZpsh)∫pp0kshQinj(u)uγ−1du, (11)

where we have introduced an exponential cut-off in the spectrum at momentum with representing the maximum momentum for protons, and is the spectral index. In our calculation, and will be kept as model parameters, and their values will be determined based on the measured all-particle spectrum.

After re-acceleration, the transport of cosmic-rays from the termination shock towards the Galactic disk also follows Equation 7. In the absence of adiabatic losses, the density of re-accelerated cosmic rays at the Earth (taken to be at ) is given by,

 NGW−CRs(p)=Qsh4πDwRshexp[−~VR2sh2Dw] (12)

The diffusion in the wind region is assumed to be much faster than near the Galactic disk as the level of magnetic turbulence responsible for particle scattering is expected to decrease with the distance away from the Galactic disk. We assume to follow the same rigidity dependence as , and take . For the wind velocity, we take the velocity gradient  km/s/kpc. This value of is within the range predicted in an earlier study using an advection-diffusion propagation model (Bloemen et al., 1993), but slightly larger than the constraint given in Strong & Moskalenko (1998). It may be noted that as long as both and are within a reasonable range, it is not their individual values that is important in determining the flux of GW-CRs, but their ratio , as can be seen from Equation 12. The larger this ratio, the more the flux will be suppressed, and vice-versa.

The distance to the termination shock can be estimated by balancing the Galactic wind ram pressure, , against the intergalactic pressure, , at the position of the termination shock, where is the mass density of the wind and represents the terminal velocity of the wind. The ram pressure is related to the total mechanical luminosity of the wind at the termination shock as, . Using this, we obtain,

 Rsh=(Lw2πPIGM~V)1/3. (13)

For Galactic wind driven by cosmic rays (Zirakashvili et al., 1996), the total mechanical luminosity of the wind cannot be larger than the total power of the cosmic rays. From Section 2.2, the total power invested in SNR-CRs (which dominates the overall cosmic-ray energy density in our model) is of the mechanical power injected by supernova explosions in the Galaxy. This corresponds to a total power of ergs s injected into SNR-CRs. Using this, and taking an intergalactic pressure of ergs cm (Breitschwerdt et al., 1991), we obtain  kpc from Equation 13. The spectral indices are taken to be the same as the source indices of the SNR-CRs listed in Table 1. Having fixed these parameter values, the spectra of the GW-CRs calculated using Equation 12 are shown in Figure 3. Spectra for the individual elements and also the total contribution are shown. The same particle injection fraction of is applied to all the elements, and the maximum proton energy corresponding to is taken as GeV. These values are chosen so that the total GW-CR spectrum reasonably agrees with the observed all-particle spectrum between and  GeV.

The GW-CRs produce a negligible contribution at low energies. This is due to the increasing effect of advection over diffusion at these energies, preventing particles from reaching the Galactic disk. Higher energy particles, which diffuse relatively faster, can overcome the advection and reach the disk more effectively. The flux suppression at low energies is more significant for heavier nuclei like iron which is due to their slower diffusion relative to lighter nuclei at the same total energy. Adding adiabatic losses to Equation 12 will lead to further suppression of the flux at low energies. But, at energies of our interest, that is above  GeV, the result will not be significantly affected as the particle diffusion time, ), is significantly less than the adiabatic energy loss time,  yr. The steep spectral cut-offs at high energies are due to the exponential cut-offs introduced in the source spectra.

### 3.2 Cosmic rays from Wolf-Rayet star explosions (WR-CRs)

While the majority of the supernova explosions in the Galaxy occur in the interstellar medium, a small fraction is expected to occur in the winds of massive progenitors like Wolf-Rayet stars (Gal-Yam et al., 2014). Magnetic fields in the winds of Wolf-Rayet stars can reach of the order of  G, and it has been argued that a strong supernova shock in such a field can lead to particle acceleration of energies up to  GeV (Biermann & Cassinelli, 1993; Stanev et al., 1993).

Since the distribution of Wolf-Rayet stars in the Galaxy is concentrated close to the Galactic disk (see e.g. Rosslowe & Crowther (2015)), the propagation of WR-CRs can also be described by Equation 2.1 with the source term replaced by , where represents the frequency of Wolf-Rayet supernova explosions per unit surface area in the Galactic disk, and the source spectrum follows Equation 2. We assume that each Wolf-Rayet supernova explosion releases a kinetic energy of  ergs, same as the normal supernova explosion in the interstellar medium. From the estimated total number of Wolf-Rayet stars of in the Galaxy and an average lifetime of  Myr for these stars (Rosslowe & Crowther, 2015), we estimate a frequency of Wolf-Rayet explosion in every years. This corresponds to Wolf-Rayet explosion in every supernova explosions occurring in the Galaxy. The source indices of the different cosmic-ray species and the propagation parameters for the WR-CRs are taken to be the same as for the SNR-CRs.

The contribution of WR-CRs to the all-particle spectrum is shown in Figure 4. The results are for two different compositions of the Wolf-Rayet winds available in the literatures: Carbon-to-helium (C/He) ratio of (top panel) and (bottom panel), given in Pollock et al. (2005). The abundance ratios of different elements with respect to helium for the two different wind compositions are listed in Table 2. In our calculation, these ratios are assumed to be proportional to the relative amount of supernova explosion energy injected into different elements. The overall normalisation of the total WR-CR spectrum and the maximum energy of the proton source spectrum are taken as free parameters. Their values are determined based on the observed all-particle spectrum between and  GeV. For , we obtain an injection energy of  ergs into helium nuclei from a single supernova explosion and a proton source spectrum cut-off of  GeV, while for , we obtain  ergs and  GeV respectively. For both the progenitor wind compositions, the total amount of energy injected into cosmic rays by a single supernova explosion is approximately times less than the total energy injected into SNR-CRs by a supernova explosion in the Galaxy. The total WR-CR spectrum for the case is dominated by helium nuclei up to  GeV, while for the case, helium nuclei dominate up to  GeV. At higher energies, carbon nuclei dominate. One major difference of the WR-CR spectra from the GW-CR spectrum (Figure 3) is the absence of the proton component, and a very small contribution of the heavy elements like magnesium, silicon and iron. Another major difference is the much larger flux of WR-CRs than the GW-CRs below  GeV. Below the knee, the total WR-CR spectrum is an order of magnitude less than the total SNR-CRs spectrum (Figure 2).

## 4 All-particle spectrum and composition of cosmic rays at high energies

The all-particle spectrum obtained by combining the contributions of SNR-CRs, GW-CRs and EG-CRs is compared with the measured data in Figure 5. For the SNR-CRs shown in the figure, we have slightly reduced the value of from  GeV (as used in Figure 2) to GeV in order to reproduce the measurements better around the knee. The extra-galactic contribution, denoted by EG-RSB93 in the figure, is taken from Rachen et al. (1993), which represents a pure proton population with a source spectrum of and an exponential cut-off at  GeV as expected from strong radio galaxies or sources with a similar cosmological evolution. Also shown in the figure are the spectra of the individual elements. The model prediction reproduces the observed elemental spectra as well as the observed features in the all-particle spectrum.

The total spectra for the two WR-CR scenarios are shown in Figure 6. For the SNR-CRs, here we take  GeV, and a slightly lower value of which corresponds to out of every supernova explosions in the Galaxy (assuming a fraction going into Wolf-Rayet supernova explosions as deduced in the previous section). The injection energy for the different elements of the SNR-CRs has been re-adjusted accordingly, so that the sum of SNR-CRs and WR-CRs for the individual elements agree with the measured elemental spectra at low energies. The values are listed in Table 3. The cosmic-ray propagation parameters are the same as in Figure 2. The predicted all-particle spectra are in good agreement with the measurements. The WR-CR scenarios are found to reproduce the second knee and the ankle better than the GW-CR model.

In Figure 7, we show the elemental fraction at high energies predicted by the GW-CR and WR-CR models. In all the models, the composition consists of a large fraction of helium nuclei over a wide energy range. The maximum helium fraction is found in the case of WR-CR (=0.1) scenario, where the fraction reaches up to at energy GeV. In contrast to common perceptions, the WR-CR scenarios predict a composition of Galactic cosmic rays dominated mainly by helium (in the case) or carbon nuclei (in the ) near the transition energy region from Galactic to extra-galactic cosmic rays. The GW-CR model predicts an almost equal contribution of helium and iron nuclei at the transition region.

The cosmic-ray composition at energies above  GeV is not quite as well-measured as at lower energies. Above  GeV, KASCADE has provided spectral measurements for groups of elements by measuring the electron and muon numbers of extensive air showers induced by cosmic rays in the Earth’s atmosphere. Several other experiments such as LOFAR, TUNKA, and the Pierre Auger Observatory have also provide composition measurements at high energies by measuring the depth of the shower maximum . Heavier nuclei interact higher in the atmosphere, resulting in smaller values of as compared to lighter nuclei. For comparison with theoretical predictions, we often use the mean logarithmic mass, , of the measured cosmic rays which can be obtained from the measured values using the relation (Hörandel, 2003b),

 ⟨lnA⟩=(Xmax−XpmaxXFemax−Xpmax)×lnAFe, (14)

where and represent the average depths of the shower maximum for protons and iron nuclei respectively given by Monte-Carlo simulations, and is the mass number of iron nuclei.

In Figure 8, the values predicted by the different models are compared with the measurements from different experiments. Although all our model predictions are within the large systematic uncertainties of the measurements, at energies above  GeV, the GW-CR model deviates from the general trend of the observed composition which reaches a maximum mean mass at  GeV, and becomes gradually lighter up to the ankle. However, in the narrow energy range of  GeV, the behaviour of the GW-CR model is in good agreement with the measurements from TUNKA, LOFAR and Yakutsk experiments which show a nearly constant composition that is different from the behaviour observed by the Pierre Auger Observatory at these energies. Understanding the systematic differences between the different measurements at these energies will be important for further testing of the GW-CR model. Up to around the ankle, the WR-CR models show an overall better agreement with the measurements than the GW-CR model. At around  GeV, the WR-CR models seem to slightly under predict the KASCADE measurements, and they are more in agreement with the TUNKA measurements. Cosmic-ray composition measured by experiments like KASCADE, which measures the particle content of air showers on the ground, is known to have a large systematic difference from the composition measured with fluorescence and Cherenkov light detectors using measurements (Hörandel, 2003b). The large discrepancy between the model predictions and the data above the ankle is due to the absence of heavy elements in the EG-CR model considered in our calculation. The effect of choosing other models of EG-CRs will be discussed in the next section.

## 5 Test with different models of extra-galactic cosmic rays

Despite of the dominance of the ankle-transition model in the general discussion, it has often been pointed out that the essential high-energy features of the cosmic ray spectrum, that is the ankle and, in part, even the second knee, can be explained by propagation effects of extra-galactic protons in the cosmologically evolving microwave background (Hillas, 1967; Berezinsky & Grigorieva, 1988; Berezinsky et al., 2006; Hillas, 2005; Aloisio et al., 2012, 2014). While the most elegant and also most radical formulation of this hypothesis, the so-called ‘proton dip model’, is meanwhile considered disfavoured by the proton fraction at the ankle measured by the Pierre Auger Observatory (Aab et al., 2014), the light composition below the ankle recently reported by the LOFAR measurement (Buitink et al., 2016) and a potential ‘light ankle’ at about  GeV found by the KASCADE-Grande experiment (Apel et al., 2013) have reinstated the interest in such models, and led to a number of ramifications, all predicting a more or less significant contribution of extra-galactic cosmic rays below the ankle. As such a component can greatly modify the model parameters, in particular the maximum energy, for the additional Galactic component – if not removing its necessity altogether – we study this effect using the WR-CR models, which show an overall best agreement with the data below the ankle, as a Galactic paradigm.

Before, however, discussing a stronger extra-galactic component below the ankle, we want to think about the minimal extra-galactic contribution we can have, if we assume the largely heavy spectrum above the ankle is all extra-galactic and consider their propagation over extra-galactic distances. To construct this ‘minimal model’, we follow di Matteo et al. (2015) and use the Monte-Carlo simulation code CRPropa 3.0 (Batista et al., 2016), which takes into account all important interaction processes undergone by EG-CRs while propagating through the CMB and the extra-galactic background light, and also the energy loss associated with the cosmological expansion. The effects of uncertainties in the simulations are discussed in Batista et al. (2015). We assume the sources to be uniformly distributed in a comoving volume, and they produce cosmic rays with a spectrum given by (di Matteo et al., 2015),

 QEG =K0Fj(EE0)−γ, EZRc (15)

where is a normalisation constant, is the injection fraction which depends on the type of the nuclei ,  GeV, is the source spectral index which is assumed to be the same for the different nuclei, and is the rigidity at which the spectrum deviates from a power law. The model parameters are determined by simultaneously fitting the cosmic-ray energy spectrum, and variance of above the ankle observed at the Pierre Auger Observatory. We adopt the CTG111CRPropa with the default TALYS photo-disintegration cross sections and the EBL model of Gilmore et al. (2012) model for our calculation (di Matteo et al., 2015), and consider that the sources inject protons, helium, nitrogen and iron nuclei. The best-fit model parameters values are ,  GV, , , and . In this model, the EG-CR spectrum below  GeV is dominated by protons and helium nuclei which are secondary products from the photo-disintegration of heavier nuclei during the propagation. At higher energies up to  GeV, the spectrum is dominated by the CNO group. Above  GeV, the spectrum exhibits a steep cut-off which is mostly due to the intrinsic cut-off in the injection spectrum, and not due to the GZK absorption during the propagation. This gives an overall best agreement with the measured data (di Matteo et al., 2015).

The first assumption we consider for an additional component of light particles below the ankle is based on the same physics, that is photo-disintegration of energetic nuclei in photon backgrounds, but considering this effect already in potentially densely photon loaded sources during acceleration. The physical motivation for this scenario is the acceleration of heavy nuclei at external/internal shocks in gamma ray bursts (Murase et al., 2008; Globus et al., 2015b), or in tidal disruption events (Farrar & Gruzinov, 2009). Two variants of this assumptions have been recently suggested: the first, by Globus et al. (2015a), assumes that diffusion losses in the source are faster than the photo-disintegration time scale over a large range of energies, leading to a significantly steeper spectrum of the secondary protons than for the escaping residual nuclei, while in the second model by Unger et al. (2015) only the highest energy particles have an escape time which is smaller than the photo-disintegration time. While the predictions of the former model for secondary protons below the ankle are phenomenologically quite similar to the extra-galactic component of Rachen et al. (1993) at these energies, that is an approximate source spectrum with a cosmological evolution , the second model Unger et al. (2015, hereafter the ‘UFA model’) predicts a strong pure-proton component concentrated only about one order of magnitude in energy below the ankle. Within their fiducial model, they consider a mix with a pure iron Galactic cosmic-ray component in Unger et al. (2015). For our study, we use results which are optimised for a pure nitrogen Galactic composition222Michael Unger, private communication., which is closer to our predicted composition for the WR-CR model () around the second knee.

A second assumption for an additional extra-galactic component is based on a universal scaling argument, which links the energetics of extra-galactic cosmic-ray sources on various scales and predicts that a dominant contribution to extra-galactic cosmic rays is expected from clusters of galaxies, accelerating a primordial proton-helium mix at their accretion shocks during cosmological structure formation (Rachen, 2016). As it has been shown already by Kang et al. (1997) that, for canonical assumptions on the diffusion coefficient around shocks (e.g. Bohm diffusion), the particle acceleration in this scenario is limited by pair-production losses in the CMB, this extra-galactic component is rather expected not to reach ultra-high energies, except for very optimistic assumptions on the acceleration process, but to be confined to energies below the ankle. As so far no detailed Monte-Carlo propagation for this model has been calculated, we use here the analytical approximation developed in Rachen (2016). Assuming that both injection and acceleration of primordial protons and helium nuclei are only dependent on particle rigidity, the model predicts a succession of a proton and helium component with increasing energy, which are fixed in relative normalisation by the know primordial abundances. The more energetic helium component sharply cuts off at the ankle, merging into the cosmic-ray spectrum produced by extra-galactic sources at smaller scales, for which acceleration even in the conservative case is not limited by CMB or other photon interactions, and thus reaches the so-called Hillas limit, , if is the typical magnetic field, and the typical size of the accelerator (Hillas, 1984). In our treatment, we hereby keep the exact cut-off energy and the total normalisation of this primordial cluster shock component as free parameters and determine them from fitting the all-particle spectrum, where we use the minimal model derived above as the second extra-galactic component extending into ultra-high energies. This model is henceforth denoted as ‘PCS model’.

In Figure 9, we present the all-particle spectrum above  GeV obtained using the three different EG-CR models – minimal model only, UFA and PCS model. The galactic contributions are from SNR-CRs and WR-CRs (). For the SNR-CRs, all the model parameters are the same as in Figure 6 (bottom). For the WR-CRs, the cut-off energy and the normalisation of the source spectrum are re-adjusted in order to produce an overall good fit to the measured spectrum and composition. They are allowed to vary in the three different cases. For the minimal model, the best-fit proton cut-off energy of the WR-CRs is found to be  GeV. This is approximately a factor larger than the value used in Figure 6. For the PCS and the UFA models, the proton cut-off energies are almost the same at  GeV, which are about a factor 1.5 less than that of the minimal model. This relaxation in the cut-off energy is due to the strong contribution of EG-CRs below the ankle in the two models. In the minimal model, the transition from Galactic to extra-galactic components occurs around the ankle, while in the PCS and UFA models, it occurs at  GeV. The variation in the injection energy of WR-CRs remain within between the three models. In Figure 9, spectra of five different mass groups are also shown. The elemental fraction of these mass groups are shown in Figure 10.

In Figure 11, we show predicted by the three EG-CRs model after adding the Galactic contribution. At energies between  GeV and  GeV, the minimal model shows a bump that follows the trend of LOFAR and the data from other experiments, but contradicts the composition data from the Pierre Auger Observatory at  GeV. The UFA model over predicts the data above the ankle as the model is also tuned to the variance of , but it is well within the systematic uncertainties (experimental as well as theoretical) as discussed in Unger et al. (2015). The sharp feature present just above  GeV in the PCS model is due to the dip in the proton spectrum (Figure 9, middle panel, black-thin-solid line) that results from the intersection of the components from galaxy clusters and the minimal model, and is partially an artefact of the simplified propagation approach applied to this model. We expect it to be much smoother for realistic propagation. At energies below  GeV, both the PCS and the UFA models produce similar results which are in better agreement with the observed trend of the composition, but do not introduce a significant improvement over the canonical extra-galactic component used in Section 4. In all the three cases for the EG-CR model, the CNO group dominates the composition of Galactic cosmic rays at the transition region from Galactic to extra-galactic cosmic rays. A clear distinction between the models would be possible from a detailed measurement of the five major mass groups shown in Figure 10, in which they all have their characteristic ‘fingerprint’: for example, around GeV the minimal model is dominated by the CNO group, the PCS model by helium, and the UFA model by protons.

Results obtained using the WR-CR () scenario are given in Appendix B. The main difference from the results of the scenario is the significant dominance of helium up to the transition energy region from Galactic to extra-galactic cosmic rays (see Figures 14 and 15). The main results and the parameter values of the different models discussed in the present work are summarised in Table 4.

## 6 Discussions

Our study has demonstrated that cosmic rays below  GeV can be predominantly of Galactic origin. Above  GeV, they are most likely to have an extra-galactic origin. We show that both the observed all-particle spectrum and the composition at high energies can be explained if the Galactic contribution consists of two components: (i) SNR-CRs which dominates the spectrum up to  GeV, and (ii) GW-CRs or preferably WR-CRs which dominates at higher energies up to  GeV. When combined with an extra-galactic component expected from strong radio galaxies or a source population with similar cosmological evolution, the WR-CR scenarios predict a transition from Galactic to extra-galactic cosmic rays at around  GeV, with a Galactic composition mainly dominated by helium or the CNO group, in contrast to most common assumptions. In the following, we discuss our results for the SNR-CRs, GW-CRs, and WR-CRs in the context of other views on the Galactic cosmic rays below  GeV, the implication of our results on the strength of magnetic fields in the Galactic halo and Wolf-Rayet stars, and also the case of a steep extra-galactic component extending below the second knee.

### 6.1 SNR-CRs

The maximum contribution of the SNR-CRs to the all-particle spectrum is obtained at a proton cut-off energy of  GeV (see Figure 2). Such a high energy is not readily achievable under the standard model of diffusive shock acceleration theory in supernova remnants for magnetic field values typical of that in the interstellar medium (see e.g. Lagage & Cesarsky, 1983). However, numerical simulations have shown that the magnetic field near supernova shocks can be amplified considerably up to times the mean interstellar value (Lucek & Bell, 2000; Reville & Bell, 2012). This is also supported by observations of thin X-ray filaments in supernova remnants which can be explained as due to rapid synchrotron losses of energetic electrons in the presence of strong magnetic fields (Vink & Laming, 2003; Parizot et al., 2006). Such strong fields may lead to proton acceleration up to energies close to the cut-off energy obtain in our study (Bell, 2004).

The main composition of cosmic rays predicted by the SNR-CRs alone looks similar to the prediction of the poly-gonato model (Hörandel, 2003a). Both show a helium dominance over proton around the knee, and iron taking over at higher energies at  GeV in the SNR-CRs, and at  GeV in the poly-gonato model.

The helium dominance is more significant in the SNR-CRs than in the poly-gonato model which is due to the flatter spectral index required to reproduce the recent measurements from CREAM and ATIC experiments with the SNR-CRs. The main difference, however, is in the total contribution above  GeV. SNR-CRs alone cannot explain the observed all-particle spectrum above  GeV. They contribute only of the observed cosmic rays at  GeV. On the other hand, in the poly-gonato model, the total contribution from elements with can explain the observed spectrum up to energies close to  GeV. This difference is mainly due to the difference in the shapes of the spectral cut-offs of particles between the two models. For the SNR-CRs, we consider a power-law with an exponential cut-off, while the poly-gonato model assumes a broken power-law with a smooth break around the cut-off (break) energy. This leads to a higher flux around the cut-off energy in the poly-gonato model. On adding GW-CRs or WR-CRs as an additional Galactic component, the composition above  GeV in our model has a large fraction of helium or a mixture of helium and CNO group, which is quite different from the prediction of the poly-gonato model where the composition is mainly dominated by iron nuclei. Our prediction (in particular, that of the WR-CR scenario) is more in agreement with the measurements from fluorescence and Cherenkov light detectors, while the poly-gonato model is in agreement with data from the measurements of air shower particles on the ground.

Recently, Globus et al. (2015a) claimed that a single Galactic component with rigidity dependent cut-off is sufficient to explain the observed all-particle spectrum when combined with an extra-galactic component. Their claim that an additional Galactic component is not needed does not contradict our claim of having one. It is simply that they assume the particle spectrum as a broken power law with an exponential cut-off which leads to an increased flux above the break energy (knee) as in the poly-gonato model. However, we have demonstrated that if one considers a power-law spectrum with an exponential cut-off which is expected for particles produced by diffusive shock acceleration process in supernova remnants (Malkov & Drury, 2001), a single component cannot explain the observed spectrum beyond the knee, and a second Galactic component is inevitable. Their single component, which they had not assigned to any specific source class, would correspond to the superposition of multiple components similar to the ones proposed in our model. Based on the physical models of the most plausible sources and the propagation of cosmic rays in the Galaxy, we show that two Galactic components are sufficient to explain the measured spectrum, but do not exclude the existence of more than two components.

### 6.2 GW-CRs

Assuming that the maximum energy of particles produced by the Galactic wind termination shock is limited by the condition that the particle diffusion length must be less than the size of the shock, the maximum energy under Bohm diffusion can be written as, , where is the magnetic field, is the shock velocity and is the shock radius. From the GW-CR parameters obtained in our study, we can take  kpc,  km s which is the terminal wind velocity, and  GeV which is the proton cut-off energy. Using these values, the magnetic field strength in the Galactic halo is estimated to be  nG. This is approximately a factor less than the value obtained assuming Parker’s magnetic field topology for the solar wind (Equation 16).

An intrinsic issue in the case of re-acceleration by Galactic wind termination shock is the difficulty to observe the re-accelerated particles in the Galactic disk because of advection by the wind flow, except for the highest energy particles, as discussed in Section 3.1. As a consequence, the spectrum in the disk may not show a continuous transition between the SNR-CRs and GW-CRs (see e.g. Zirakashvili & Völk, 2006). This effect is actually visible in the predicted spectra of the individual elements shown in Figure 5. However, we notice that the superposition of the individual spectra smears out this effect in the all-particle spectrum. Nevertheless, in order to avoid this effect, Zirakashvili & Völk (2006) considered termination shocks which are stronger near the Galactic poles and weaker towards the Galactic equator, unlike in our study where the shocks are considered to have equal strengths in all the directions. In their configuration, the maximum energy of particles decreases from the poles towards the equator, and therefore, the superposition of spectra from different colatitudes produces a continuity in the total spectrum. Another consideration is the particle re-acceleration by spiral shocks in the Galactic wind which are formed by the interaction between fast winds originating from the Galactic spiral arms and slow winds from the interarm regions (Völk & Zirakashvili, 2004). These shocks, which can be formed at distances of  kpc, can accelerate SNR-CRs up to  GeV. An alternative possibility is the re-acceleration by multiple shock waves in the Galactic wind generated by time dependent outflows of gas from the Galactic disk (Dorfi & Breitschwerdt, 2012). These shocks, which are long-lived like the termination shocks, can accelerate particles up to  GeV in the lower Galactic halo. An attractive feature of this model is the advection of particles downstream of the shocks towards the Galactic disk, thereby, resolving the difficulty of observing the re-accelerated particles in the disk. Despite having different features, the cosmic-ray composition predicted by all these different models in the energy range of  GeV are expected to be similar to the result presented here since they consider the same seed particles (cosmic rays from the Galactic disk) for re-acceleration as in our study. Below  GeV where the GW-CRs are significantly suppressed in our case, the other wind models discussed above will give a different result.

### 6.3 WR-CRs

The prediction of a large helium fraction and a small iron fraction between around and  GeV by the WR-CR () model seems to be in agreement with new measurements from the LOFAR radio telescope (Buitink et al., 2016), and the Pierre Auger Observatory (Aab et al., 2014). These measurements have revealed a strong light component, and an almost negligible iron component above  GeV. In Figure 12, the elemental fraction predicted by the WR-CR () model combined with the PCS model for the EG-CRs is compared with the best-fit composition of the LOFAR data for four mass groups: , , , and . The model predictions are found to show a good agreement with the data.

Using the maximum energy of particles for the WR-CRs, it is possible to estimate the strength of the magnetic field at the surface of Wolf-Rayet stars. Assuming that the magnetic field configuration in the Wolf-Rayet winds follows Parker’s model (Parker, 1958), the toroidal magnetic field strength near the equatorial plane of the star at the position from the star follows the relation,

 B=B0ωR2⋆VwRw, (16)

where is the magnetic field at the surface of the star, is the angular rotation velocity, is the radius of the star, and is the wind velocity. Using the relation for the maximum energy as in the case of the GW-CRs, and the proton cut-off energy of  GeV for the WR-CRs () obtained using the PCS/UFA model, we get  G cm, where we take the shock velocity (Soderberg et al., 2012). Using this value of in Parker’s magnetic field configuration (Equation 16) by taking and other Wolf-Rayet star parameters as  cm,  s, and  km s (Berezhko & Völk, 2000), we obtain the magnetic field at the surface of the star as  G. Such a strong magnetic field was also predicted in an earlier study by Biermann & Cassinelli (1993), and is found to be in agreement with recent magnetic field measurements from Wolf-Rayet stars. Based on an upper limit of  G in the observable parts of Wolf-Rayet winds, Chevrotieŕe et al. (2013) estimated an upper limit for the surface magnetic field of  G. An even stronger field in the wind, up to  G, has been reported (Chevrotieŕe et al., 2014), which indicates that the surface magnetic field of these stars can go well above the order of  G.

From the total energy of  ergs injected into WR-CRs by a single supernova explosion, and the explosion rate of Wolf-Rayet stars in the Galaxy of  yr, we estimate the total power injected into WR-CRs as  ergs s. This is approximately a factor less than the power injected into SNR-CRs by supernova explosions in the interstellar medium. The required amount of supernova explosion energy injected into helium nuclei for WR-CRs is about times that of the SNR-CRs. This indicates that the average abundance of helium nuclei swept up by supernova shocks in the Wolf-Rayet winds must be higher than the helium abundance present in the interstellar medium if the particle injection fraction and the acceleration efficiency of the shocks are the same for the SNR-CRs and the WR-CRs.

Our results for the WR-CRs are obtained by assuming that the particle injection fraction into the shocks is the same for all the different elements. The injection fraction may depend on the type of the element, and the nature of this dependence is not quite understood. By taking the ratio of the SNR-CRs source spectra (Equation 2) at a fixed rigidity to the known Solar system elemental abundances (Lodders & Palme, 2009), we estimate the relative injection fraction of particles for the different elements. Applying these relative injection fractions to the WR-CRs, we find that the composition is significantly dominated by carbon nuclei, in contrast to the results shown in Figure 4 where the composition is mainly dominated by helium or a mixture helium and carbon nuclei. Thus, the contribution of WR-CRs in this case is strongly constrained by the measured carbon spectrum at low energies. The all-particle spectrum for this case, after adding the contributions of SNR-CRs and EG-CRs, underpredicts the measured data between the second knee and the ankle. This problem might be resolved if we consider that both GW-CRs and WR-CRs contribute at the same time. In future, we will explore the parameter space of this combined scenario.

### 6.4 Comparison with Hillas’s ‘Component B’

Bell & Lucek (2001) showed that magnetic field upstream of supernova shock fronts can be amplified non-linearly by cosmic rays up to many times the pre-shock magnetic field. They showed that these highly amplified magnetic fields can facilitate cosmic-ray acceleration up to energies  GeV for supernova shocks expanding in the interstellar medium, even higher by an order of magnitude for shocks expanding into pre-existing stellar winds. Based on the Bell-Lucek’s version of diffusive shock acceleration, Hillas (2005) proposed a second Galactic component ‘Component B’, produced by Type II supernova remnants in the Galaxy expanding into dense slow winds of the preceding red supergiants, to accommodate for the observed cosmic rays above  GeV. In the Hillas (2005) model, a Galactic component ‘Component A’, produced by Type Ia supernova remnants in the Galaxy, dominates the all-particle energy spectrum below  GeV. The ‘Component A’ has a similar composition to the SNR-CRs in our model, but the ‘Component B’ has a large iron fraction in contrast to the WR-CR component in our model which is dominated mostly by helium or a mixture of helium and CNO group with a small iron fraction. Between and  GeV, the predicted all-particle spectrum in Hillas (2005) consists of a significant iron fraction which may be in agreement with the data when mixed with a strong extra-galactic proton component, but is in tension with the small iron fraction () preferred by the recent measurements of LOFAR (Buitink et al., 2016) and the Pierre Auger Observatory (Aab et al., 2014). These new measurements disfavour the general view that the Galactic component above the second knee is dominated by heavy (iron) nuclei.

### 6.5 A steep EG-CR component extending below the second knee

An alternative to the modulation of EG-CRs by the Galactic wind is the ‘magnetic horizon effect’ (Stanev et al., 2000; Lemoine, 2005; Aloisio & Berezinsky, 2005), which leads to a flattening of the extra-galactic spectrum below an energy where the diffusive propagation distance in a partly turbulent extra-galactic magnetic field, over the time scale set by energy losses of the cosmic rays through interactions with ambient photon backgrounds, gets below the average distance of cosmic ray sources. Assuming a relatively strong ( nG) extra-galactic field with a constant coherence length extending over the entire universe, this effect could set in at around  GeV, effectively cutting off the extra-galactic component at lower energies slightly below the ankle (Aloisio et al., 2012), or even above (Mollerach & Roulet, 2013). However, more detailed treatments in the context of large scale structure formation (Kotera & Lemoine, 2008), have indicated that this effect is much less efficient due to the large voids in the universe which are essentially free of magnetic fields. As shown recently in detailed simulations, the magnetic horizon effect should play virtually no role above the second knee for any type of nuclei, and for protons in some extra-galactic magnetic field scenarios, not even above the knee (Batista & Sigl, 2014).

We point out that neither the Galactic wind nor the magnetic horizon effects discussed above prevent a hard extra-galactic component, like the light component with as indicated by the KASCADE-Grande measurements above  GeV (Apel et al., 2013), from contributing around the second knee as such a hard component will be already consistent with the measured data at low energies. Even if such a hard extra-galactic component is present, an additional Galactic component will still be required as the extra-galactic component will remain subdominant in the all-particle spectrum below  GeV.

An additional problem for EG-CRs with an overall spectrum steeper than is that, if one assumes that they fill the extra-galactic space homogeneously with energies from  GeV to  GeV, it contains more energy than the gravitational binding energy released in the universe during structure formation (Rachen, 2016). Using realistically low efficiencies for this energy – which is, besides the lower overall nuclear binding energy released in fusion by all primordial baryonic matter going into stars, the only fundamental energy budget present in the late universe – to be converted into cosmic rays, one can conclude that spectral indices as discussed here for a dominant extra-galactic component below the second knee cannot easily be reconciled with this energy budget, no matter which kind of sources one proposes. Mainly on the basis of this argument, together with the difficulties of a sufficient spectral modification at low energies discussed above, we consider a dominantly extra-galactic explanation of cosmic rays below  GeV as implausible.

## 7 Conclusions

We have demonstrated that a single Galactic component with progressive energy cut-offs in the individual spectra of different elements, and describing the low-energy measurements below  GeV from balloon and satellite-borne experiments, cannot explain simultaneously the knee and the second knee observed in the all-particle spectrum. We show that a two-component Galactic model, the first component dominating up to  GeV and the second component dominating in the range of  GeV, can explain almost all observed features in the all-particle spectrum and composition when combined with an extra-galactic component dominating above  GeV. Discussing two different scenarios for the second Galactic component, we find that a contribution of Wolf-Rayet supernovae explain best both the measured energy spectrum and composition. Our main result is that this component predicts a Galactic contribution at and above the second knee which is mainly dominated by helium or a mixture of helium and CNO nuclei, and is consistent with a ‘regular’ extra-galactic contribution from sources with a flat spectral index and a cosmological evolution typical for AGNs or star formation. Using re-acceleration at the Galactic wind termination shock as a second Galactic component also allows to fit the all-particle energy spectrum, but not the observed composition very well. Tests of the two-component Galactic model using different hypotheses for a significant extra-galactic cosmic-ray component below the ankle, do neither significantly improve nor deteriorate this result, mostly because both the Galactic and extra-galactic components have a rather light composition, and contain little or no heavy nuclei like iron, in contrast to common assumptions. In all cases, the transition from Galactic to extra-galactic cosmic rays occurs between the second knee and the ankle, and we see neither the need nor a theoretical case for an extra-galactic component significantly contributing at or below  GeV. Our findings are in agreement with recent measurements from LOFAR and the Pierre Auger Observatory, which have revealed a strong light component and a rather low iron fraction between and  GeV. A clear distinction of the various discussed Galactic and extra-galactic scenarios would be possible if we could separately measure the spectra of at least four major mass groups, that is protons, helium, CNO, and heavier, at energies between the second knee and the ankle.

## Acknowledgements

We wish to thank Michael Unger, Glennys Farrar and Luis Anchordoqui, for their comments on the manuscript, and for providing us results for the UFA model. We also thank Peter Biermann for his insightful comments, and Torsten Enßlin and Christoph Pfrommer for helping to improve the manuscript during a discussion at the meeting of ‘International Team 323’ at the International Space Science Institute in Bern. ST wishes to thank Onno Pols for discussions on Wolf-Rayet stars. AvV acknowledges financial support from the NWO Astroparticle physics grant WARP. We furthermore acknowledge financial support from an Advanced Grant of the European Research Council (grant agreement no. 227610), European Union’s Horizon 2020 research and innovation programme (grant agreement no. 640130), the NWO TOP grant (grant agreement no. 614.001.454), and the Crafoord Foundation (grant no. 20140718).

## References

• Aab et al. (2014) Aab, A., et al. (Pierre Auger Collaboration), 2014, PRD, 90 122006
• Aartsen et al. (2013) Aartsen, M. G., et al. 2013, PRD, 88, 042004
• Abbasi et al. (2009) Abbasi, R. U., et al. 2009, APh, 32, 53
• Abdo et al. (2009) Abdo, A. A., et al. 2009, PRL, 102, 181101
• Abu-Zayyad et al. (2013) Abu-Zayyad, T., et al. 2013, ApJ, 768, L1
• Aharonian et al. (2006) Aharonian, F. A., et al. 2006, ApJ, 636, 777
• Aharonian et al. (2007) Aharonian, F. A., et al. 2007, A&A, 464, 235
• Aharonian et al. (2008) Aharonian, F. A., et al. 2008, A&A, 477, 353
• Adriani et al. (2011) Adriani, O., et al. 2011, Science, 332, 69
• Adriani et al. (2014) Adriani, O., et al. 2014, Physics Reports, 544, 323
• Aguilar et al. (2013) Aguilar, M., et al. 2013, PRL, 110, 141102
• Aguilar et al. (2014) Aguilar, M., et al. 2014, PRL, 113, 121102
• Aguilar et al. (2015a) Aguilar, M., et al. 2015a, PRL, 114, 171103
• Aguilar et al. (2015b) Aguilar, M., et al. 2015b, PRL, 115, 211101
• Ahn et al. (2009) Ahn, H. S., Allison, P. S., Bagliesi, M. G., et al. 2009, ApJ, 707, 593
• Albert et al. (2007) Albert, J., et al. 2007, A&A, 474, 937
• Aloisio & Berezinsky (2005) Aloisio, R., & Berezinsky, V. S., 2005, ApJ, 625, 249
• Aloisio et al. (2012) Aloisio, R., Berezinsky, V. S, & Gazizov, A., 2012, APh, 39, 129
• Aloisio et al. (2014) Aloisio, R., Berezinsky, V. S, & Blasi, P., 2014, JCAP, 10, 020
• Amenomori et al. (2008) Amenomori, M., et al. 2008, ApJ, 678, 1165
• Antoni et al. (2005) Antoni, T., et al. 2005, APh, 24, 1
• Apel et al. (2013) Apel, W. D., et al. 2013, PRD, 87, 081101
• Axford (1994) Axford, W. I., 1994, ApJS, 90, 937
• Baade & Zwicky (1934) Baade, W., & Zwicky, F., 1934, PNAS, 20, 259
• Batista & Sigl (2014) Batista, R. A, & Sigl, G., 2014, JCAP, 1411, 031
• Batista et al. (2015) Batista, R. A., et al. 2015, JCAP, 10, 063
• Batista et al. (2016) Batista, R. A., et al. 2016, JCAP, 05, 038
• Bell (1978) Bell, A. R. 1978, MNRAS 182, 147
• Bell (2004) Bell, A. R., 2004, MNRAS, 353, 550
• Bell & Lucek (2001) Bell, A. R., & Lucek, S. G., 2001, MNRAS, 321, 433
• Berezhko & Völk (2000) Berezhko, E. G., & Völk, H. J., 2000, A&A, 357, 283
• Berezhko et al. (2003) Berezhko, E. G., Ksenofontov, L. T., Ptuskin, V. S., Zirakashvili, V. N., Völk, H. J., 2003, A&A, 410, 189
• Berezhnev et al. (2013) Berezhnev, S. F., et al. 2013, in Proc. 33rd ICRC, Paper ID 326
• Berezinsky & Grigorieva (1988) Berezinsky, V. S., & Grigorieva, S. I.,1988, A&A, 199, 1
• Berezinsky et al. (2004) Berezinsky, V. S., Grigorieva, S. I., & Hnatyk, B. I., 2004, APh, 21, 617
• Berezinsky et al. (2006) Berezinsky, V. S., Gazizov, A., & Grigorieva, S. I., 2006, PRD, 74, 043005
• Biermann & Cassinelli (1993) Biermann, P. L., & Cassinelli, J. P., 1993, A&A, 277, 691
• Biermann et al. (2010a) Biermann, P. L., Becker, J. K., Dreyer, J., Meli, A., Seo, E., & Stanev, T., 2010a, ApJ, 725, 184
• Biermann et al. (2010b) Biermann, P. L., Becker J. K., Caceres, G., Meli, A., Seo, E., & Stanev, T., 2010b, ApJ, 725, 184
• Bird et al. (1994) Bird, D. J., et. al. 1994, ApJ, 424, 491
• Blandford & Ostriker (1978) Blandford, R., & Ostriker, J. P., 1978, ApJ, 221, L29
• Blasi et al. (2012) Blasi, P., Amato, E., & Serpico, P. D., 2012, PRL, 109, 061101
• Bloemen et al. (1993) Bloemen, J. B. G. M., Dogiel, V. A., Dorman, V. L., & Ptuskin, V. S., 1993, A&A, 267, 372
• Buitink et al. (2016) Buitink, S., et al. (LOFAR Collaboration) 2016, Nature, 531, 70
• Breitschwerdt et al. (1991) Breitschwerdt, D., McKenzie, J. F., & Vólk, H. J., 1991, A&A, 245, 79
• Breitschwerdt et al. (2002) Breitschwerdt, D., Dogiel, V. A., & Vólk, H. J., 2002, A&A, 385, 216
• Caprioli et al. (2011) Caprioli, D., Blasi, P., & Amato, E., 2011, APh, 34, 447
• Chevrotieŕe et al. (2013) de la Chevrotieŕe, A., St-Louis, N., & Moffat, A. F. J., (MiMeS Collaboration) 2013, ApJ, 764, 171
• Chevrotieŕe et al. (2014) de la Chevrotieŕe, A., St-Louis, N., & Moffat, A. F. J., (MiMeS Collaboration) 2014, ApJ, 781, 73
• di Matteo et al. (2015) di Matteo, A. for the Pierre Auger Collaboration, PoS (ICRC2015), 249; arXiv:1509.03732
• Dorfi & Breitschwerdt (2012) Dorfi, E. A., & Breitschwerdt, D., 2012, A&A, 540, A77
• Drury (1983) Drury, L. O, 1983, Reports on Progress in Physics, 46, 973
• Engelmann et al. (1990) Engelmann, J. J., Ferrando, P., Soutoul, A., Goret, P., & Juliusson, E., 1990, A&A, 233, 96
• Erlykin & Wolfendale (2012) Erlykin, A. D. Wolfendale, A. W., 2012, APh, 35, 449
• Farrar & Gruzinov (2009) Farrar, G. R., & Gruzinov, A., 2009, ApJ, 693, 329
• Gal-Yam et al. (2014) Gal-Yam, A., et al. 2014, Nature, 509, 471
• Ghia et al. (2015) Ghia, P. L., for the Pierre Auger Collaboration, PoS (ICRC2015), 034; arXiv:1509.03732
• Gilmore et al. (2012) Gilmore, R. C., et al. 2012, MNRAS, 422, 3189
• Globus et al. (2015a) Globus, N., Allard, D., & Parizot, E., 2015, PRD, 92, 021302
• Globus et al. (2015b) Globus, N., Allard, D., Mochkovitch, R., & Parizot, E., 2015, MNRAS, 451, 751
• Green (2009) Green, D. A., 2009, BASI, 37, 45
• Greisen (1966) Greisen, K., 1966, PRL, 16, 748
• Hillas (1967) Hillas, A. M., 1967, Physics Letters A, 24, 677
• Hillas (1984) Hillas, A. M., 1984, ARA&A, 22, 425
• Hillas (2005) Hillas, A. M., 2005, Journal of Physics G: Nuclear and Particle Physics, 31, R95
• Hörandel (2003a) Hörandel, J. R., 2003a, APh, 19, 193
• Hörandel (2003b) Hörandel, J. R., 2003b, Journal of Physics G: Nuclear and Particle Physics, 29, 2439
• Hörandel (2006) Hörandel, J. R., 2006, JPhCS, 47, 41
• Jokipii & Morfill (1987) Jokipii, J. R., & Morfill, G., 1987, 312, 170
• Jones et al. (2001) Jones, F. C, Lukasiak, A., Ptuskin, V., & Webber, W., 2001, ApJ, 547, 264
• Kampert & Unger (2012) Kampert, K. H., & Unger, M., 2012, 35, 660
• Kang et al. (1996) Kang, H., Ryu, D., & Jones, T. W., 1996, ApJ, 456, 422
• Kang et al. (1997) Kang, H., Rachen, J. P., & Biermann, P. L., 1997, MNRAS, 286, 257
• Kelner et al. (2006) Kelner, S. R., Aharonian, F. A., & Bugayov, V. V., 2006, PRD 74, 034018
• Knurenko & Sabourov (2010) Knurenko, S., & Sabourov, A., 2010, Proc. XVI ISVHECRI
• Kotera & Lemoine (2008) Kotera, K., & Lemoine, M., 2008, PRD, 77, 023005
• Krymskii (1977) Krymskii, G. F., 1977, Akad. Nauk SSSR Dokl., 234, 1306
• Lagage & Cesarsky (1983) Lagage P. O., & Cesarsky C. J., 1983, A&A, 125, 249
• Lemoine (2005) Lemoine, M., 2005, PRD, 71, 083007
• Lerche & Schlickeiser (1982a) Lerche, I., & Schlickeiser, R., 1982a, A&A, 116, 10
• Lerche & Schlickeiser (1982b) Lerche, I., & Schlickeiser, R., 1982b, MNRAS, 201, 1041
• Letaw et al. (1983) Letaw, J. R., Silberberg, R., & Tsao, C. H., 1983, ApJS, 51, 271
• Lucek & Bell (2000) Lucek S. G., & Bell A. R., 2000, MNRAS, 314, 65
• Lodders & Palme (2009) Lodders, K., & Palme, H., 2009, Meteoritics and Planetary Science Supplement, 72, 5154
• Malkov & Drury (2001) Malkov M. A., Drury L. O’C, 2001, Reports on Progress in Physics, 64, 429
• Mannheim et al. (2001) Mannheim, K., Protheroe, R. J. and Rachen, J. P., 2001, PRD, 63, 023003
• Mollerach & Roulet (2013) Mollerach, S., & Roulet, E., 2013, JCAP, 1310, 013
• Müller et al. (1991) Müller, D., Swordy, S. P., Meyer, P., L’Heureux, J., & Grunsfeld, J. M., 1991, ApJ, 374, 356
• Muraishi et al. (2005) Muraishi, H., Yanagita, S. & Yoshida, T., 2005, Progress of Theoretical Physics, 113, 721
• Murase et al. (2008) Murase, K., Ioka, K., Nagataki, S. & Nakamura, T., 2008, PRD, 78, 023005
• Obermeier et al. (2011) Obermeier, A., Ave, M., Boyle, P., et al., 2011, ApJ, 742, 14
• Ohira et al. (2011) Ohira, Y., Murase, K., & Yamazaki, R., 2011, MNRAS, 410, 1577
• Panov et al. (2007) Panov, A. D., et al. 2007, Bull. Russ. Acad. Sci., Vol. 71, No. 4, pp. 494
• Parizot et al. (2006) Parizot, E., Marcowith, A., Ballet, J., & Gallant, Y. A. 2006, A&A, 453, 387
• Parker (1958) Parker, E. N., 1958, ApJ, 128, 664
• Pollock et al. (2005) Pollock, A. M. T., Corcoran, M. F., Stevens, I. R., & Williams, P. M., 2005, ApJ, 629, 482
• Porcelli et al. (2015) Porcelli, A. for the Pierre Auger Collaboration, PoS (ICRC2015), 420; arXiv:1509.03732
• Ptuskin et al. (1993) Ptuskin, V. S., Rogovaya, S. I., Zirakashvili, V., Chuvilgin, L. G, et al. 1993, A&A, 268, 726
• Ptuskin et al. (2010) Ptuskin, V. S., Zirakashvili, V., & Seo, E. S., 2010, ApJ, 718, 31
• Ptuskin et al. (2013) Ptuskin, V. S., Zirakashvili, V., & Seo, E. S., 2013, ApJ, 763, 47
• Rachen & Biermann (1993) Rachen, J. P., & Biermann, P. L., 1993, A&A, 272, 161
• Rachen et al. (1993) Rachen, J. P., Stanev, T., & Biermann, P. L., 1993, A&A, 273, 377
• Rachen (2016) Rachen, J. P., 2016, Proc. 28th Texas Symposium on Relativistic Astrophysics, Geneva, Switzerland, 13-18 December 2015, Electronic Proceedings #230
• Reville & Bell (2012) Reville, B., & Bell, A. R. 2012, MNRAS, 419, 2433
• Rosslowe & Crowther (2015) Rosslowe, C. K., & Crowther, P. A., 2015, arXiv:1412.0699
• Schulz et al. (2013) Schulz, A. for the Pierre Auger Collaboration, 2013, in Proc. 33 ICRC, Rio de Janeiro; arXiv:1307.5059
• Soderberg et al. (2012) Soderberg, A. M., et al. 2012, ApJ, 752, 78
• Stanev et al. (1993) Stanev, T., Biermann, P. L., & Gaisser, T. K., 1993, A&A, 274, 902
• Stanev et al. (2000) Stanev, T., Engel, R. Múcke, A. Protheroe, R. J., & Rachen, J. P., 2000, PRD, 62, 093005
• Strong & Moskalenko (1998) Strong, A. W., & Moskalenko, I. V., 1998, ApJ, 509, 212
• Swordy et al. (1990) Swordy, S. P., Müller, D., Meyer, P., L’Heureux, J., & Grunsfeld, J. M., 1990, ApJ, 349, 625
• Thoudam (2008) Thoudam, S., 2008, MNRAS, 388, 335
• Thoudam & Hörandel (2012) Thoudam, S., & Hörandel, J. R. 2012, MNRAS, 421, 1209
• Thoudam & Hörandel (2013) Thoudam, S., & Hörandel, J. R. 2013, MNRAS, 435, 2532
• Thoudam & Hörandel (2014) Thoudam, S., & Hörandel, J. R. 2014, A&A, 567, A33
• Tomassetti (2012) Tomassetti, N. 2012, ApJ, 752, L13
• Unger et al. (2015) Unger, M., Farrar, G. R., & Anchordoqui, L. A. 2015, PRD, 92, 123001
• Vink & Laming (2003) Vink J., & Laming J. M., 2003, ApJ, 584, 758
• Völk & Biermann (1988) Völk, H. J. & Biermann, P. L., 1988, ApJ, 333, L65
• Völk & Zirakashvili (2004) Völk, H. J. & Zirakashvili, V. N, 2004, A&A, 417, 807
• Wandel (1988) Wandel, A., 1988, A&A, 200, 279
• Wandel et al. (1987) Wandel, A., Eichler, D. S., Letaw, J. R., Silberberg, R., & Tsao, C. H., 1987, ApJ, 316, 676
• Waxman (1995) Waxman, E., 1995, PRL, 75, 386
• Yoon et al. (2011) Yoon, Y. S., et al. 2011, ApJ, 728, 122
• Yuan et al. (2011) Yuan, Q., Zhang, B., & Bi, X. -J, 2011, PRD 84, 043002
• Zatsepin & Kuzmin (1996) Zatsepin, G. T., & Kuzmin V. A, 1966, JETPh Lett., 4, 78
• Zirakashvili et al. (1996) Zirakashvili, V., Breitschwerdt, D., Ptuskin, V. S. & Völk, H. J., 1996, A&A, 311, 113
• Zirakashvili & Völk (2006) Zirakashvili, V., & Völk, H. J., 2006, AdSpR, 37, 1923

## Appendix A Derivation of Equation 3.1

The Green’s function, , of Equation 7 satisfies,

 ∇.(Dw∇G−VG)+∂∂p{∇.V3pG}=−δ(r−r′)δ(p−p′). (17)

In rectangular coordinates, the above equation can be written as,

 Dw∂2G∂x2+Dw∂2G∂y2+Dw∂2G∂z2−~V∂∂x(xG)−~V∂∂y(yG) −~V∂∂z(zG)+∂∂p(~VpN)=−δ(x−x′)δ(y−y′)δ(z)δ(p−p′), (18)

where we have written with , and representing the unit vectors along the , and directions. Following a similar procedure adopted in Lerche & Schlickeiser (1982b), we express,

 G(x,x′,y,y′,z,z′,p,p′)= ∫∞−∞dkx∫∞−∞dky∫∞−∞dkz¯G(kx,x′,ky,y′,kz,z′,p,p′) ×eikx(x−x′)eiky(y−y′)eikz(z−z′), (19)

and,

 δ(x