© ESO, 2016

1. Introduction

Until a decade ago, the cosmic ray spectrum from ~ 10 GeV to ~ 10¹¹ GeV was seen as a power law with two main features: a steepening from a spectral index γ ≈ −2.7 to γ ≈ −3.1 at about 3 × 10⁶ GeV, commonly called the “knee”, and a flattening back to γ ≈ −2.7 at about 4 × 10⁹ 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 10⁵⁻⁶ GeV (Lagage & Cesarsky 1983; Axford 1994). This may extend to ~ 10⁸ 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 10⁹  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 ~ 10¹¹ 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.

Recent measurements of cosmic rays by several new generation experiments have severely challenged this simple view. At low energies, below ~ 10⁶ GeV, satellite and balloon-borne experiments such as ATIC-2 (Panov et al. 2007), CREAM (Yoon et al. 2011), TRACER (Obermeier et al. 2011), PAMELA (Adriani et al. 2014), AMS-02 (Aguilar et al. 2014, 2015a,b), and Fermi-LAT (Abdo et al. 2009) have measured the energy spectra of various elements of cosmic rays ranging from protons to heavier nuclei such as iron as well as the leptonic component of cosmic rays, and anti-particles such as positrons and anti-protons. Some of their results, for example the rise of the positron fraction above ~ 10 GeV (Aguilar et al. 2013), the harder energy spectrum of helium nuclei with respect to the proton spectrum (Adriani et al. 2011), and the spectral hardening of both the proton and helium nuclei at TeV energies (Yoon et al. 2011), are difficult to explain using standard models of cosmic-ray acceleration in supernova remnants and their subsequent propagation in the Galaxy. At high energies, that is above ~ 10⁶ GeV, ground-based experiments such as KASCADE-Grande (Apel et al. 2013), the Tibet III array (Amenomori et al. 2008), IceTop (Aartsen et al. 2013), the Pierre Auger Observatory (Ghia et al. 2015) and the Telescope Array (Abu-Zayyad et al. 2013) have carried out detailed measurements of the all-particle energy spectrum and the composition of cosmic rays. First, they confirm a third major break in the spectrum, a steepening to γ ≈ −3.3 above about 10⁸ GeV, which has been suggested before both by the Fly’s Eye stereo energy spectrum (Bird et al. 1994) and theoretical arguments about the structure of the ankle (Berezinsky & Grigorieva 1988; Rachen et al. 1993). It has anatomically been named the “second knee” (Hörandel 2006). While this still fits with the original view, the cosmic-ray composition measurements at these energies pose a severe challenge: instead of gradually becoming heavier as expected, the data show that the composition reaches a maximum mean mass at energies around 6 × 10⁷ GeV, and then becomes gradually lighter again up to the ankle. Finally, above the ankle the composition becomes heavier again. It has been shown that the observed spectrum and composition at the highest energies can be explained by a mixed-composition extra-galactic source spectrum with progressive cutoffs at ~ Z × 5 × 10⁹ GeV, where Z is the nuclear charge (Aloisio et al. 2014). This would imply that there is no significant impact of the GZK effect in cosmic ray propagation except through photo-disintegration of nuclei. In addition, the measurement of an ankle-like feature in the light component of cosmic rays at ~ 10⁸ GeV by the KASCADE-Grande experiment (Apel et al. 2013), and the new revelation of a strong light component and a very small iron component by the LOFAR measurements between ~ (1−4) × 10⁸ GeV (Buitink et al. 2016), and by the Pierre Auger Observatory above ~ 7 × 10⁸ GeV (Aab et al. 2014) add further challenges to the standard model.

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 3 × 10⁷ 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 ~ 40% 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 80% 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 Sect. 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 ~ 10⁷ GeV. We then revisit two possibilities for a second Galactic component in Sect. 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 Sect. 4. In Sect. 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 Sect. 6, we present a discussion of our results and their implications, and other views on the cosmic rays below 10⁹ GeV, followed by our conclusions in Sect. 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 2 (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 ~ 10⁴² ergs s^-1 injected by supernova explosions into the Galaxy, considering a supernova explosion energy of ~ 10⁵¹ ergs and an explosion frequency of ~ 1/30 yr^-1, 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 100 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 ~ 1 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, (1)where we have adopted a cylindrical geometry for the propagation region described by the radial r and vertical z coordinates with z = 0 representing the Galactic plane. We assume the region to have a constant halo boundary at z = ± L, 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 ~ L 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. N(r,z,p) represents the differential number density of the cosmic-ray nuclei with momentum/nucleon p, and Q(r,p)δ(z) 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 Eq. (1). The diffusion coefficient D(ρ) is assumed to be a function of the particle rigidity ρ as, D(ρ) = D₀β(ρ/ρ₀)^a, where D₀ is the diffusion constant, β = v/c with v(p) and c representing the velocity of the particle and the velocity of light respectively, ρ₀ = 3 GV is a constant, and a is the diffusion index. The rigidity is defined as ρ = Apc/Ze, where A and Z represent the mass number and the charge number of the nuclei respectively, and e is the charge of an electron. The second term in Eq. (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, σ(p) is the inelastic interaction cross-section, and ξ corresponds to the rate of re-acceleration. We take , where V = 4πℜ³/ 3 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 Eq. (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 ~ 1 GeV/nucleon, are not included explicitly in the transport equation. Instead, we introduce a low-momentum cut-off, p₀ ~ 100 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 s. The source term Q(r,p) can be expressed as , where H(m) = 1(0) for m> 0 ( < 0) represents a Heaviside step function, and the source spectrum Q(p) 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 (2)where Q₀ is a normalisation constant which is proportional to the amount of energy f channelled into cosmic rays by a single supernova event, q is the spectral index, and p_c is the cut-off momentum for protons. The exponential cut-off in Eq. (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 Z 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 Eq. (1), we assume that both cosmic-ray sources and interstellar matter are distributed in the disk (i.e. at z = 0). The distributions are assumed to be uniform, and extended up to a radius R.

Recalling the analytical solution of Eq. (1) derived in Thoudam & Hörandel (2014), the cosmic-ray density at the position r = 0 for p>p₀ follows, (3)where J₁ is a Bessel function of order 1, and the functions B and are given by, (4)From Eq. (3), the cosmic-ray density at the Earth can be obtained by taking z = 0 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 (D₀,a) and the re-acceleration parameters (η,s) have been obtained to be, D₀ = 9 × 10²⁸ cm² s^-1, a = 0.33, η = 1.02, and s = 4.5 (Thoudam & Hörandel 2014). We adopt these values in our present study. The supernova remnant radius is taken to be ℜ = 100 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 L. We choose L = 5 kpc, which gives an averaged surface density of atomic hydrogen of atoms cm^-2 (Thoudam & Hörandel 2013). An extra 10% is further added to to account for the helium abundance in the interstellar medium. The radial extent of the source distribution is taken as R = 20 kpc. Each supernova explosion is assumed to release a total kinetic energy of 10⁵¹ ergs, and the supernova explosion frequency is taken as SNe Myr^-1 kpc^-2. The latter corresponds to a rate of ~ 3 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 Fig. 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 (q,f) 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 2.21−2.29, and out of the total of 8% of the supernova explosion energy channelled into SNR-CRs, the largest fraction goes into protons at the level of 6.95%, followed by helium nuclei with 0.79%. 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 Fig. 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 ~ 10⁶ GeV. The calculation assumes an exponential cut-off for the proton source spectrum at E_c = 4.5 × 10⁶ GeV, and for the heavier nuclei at ZE_c. This value of E_c, which is obtained by comparing the predicted all-particle spectrum with the observed all-particle spectrum as shown in Fig. 2, represents the maximum E_c value permitted by the measurements. While obtaining the all-particle spectrum shown in Fig. 2, we also include contributions from the sub-dominant primary cosmic-ray elements (Z< 26), calculated using elemental abundances at 10³ GeV given in Hörandel (2003a) and a source index of 2.25. Their total contribution amounts up to ~ 8% of the all-particle spectrum. The predicted all-particle spectrum agrees with the data up to ~ 2 × 10⁷ GeV, and reproduces the observed knee at the right position. Choosing E_c values larger than 4.5 × 10⁶ 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 E_c value does not rely on the proton measurements at high energies, it can be noticed from Fig. 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 Fig. 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 ~ 10⁵ GeV. Such a trend is also consistent with the KASCADE measurements above ~ 10⁶ GeV (see Fig. 1). Based on our prediction, helium nuclei dominate the all-particle spectrum up to ~ 1.5 × 10⁷ GeV, while above, iron nuclei dominate. The maximum energy of SNR-CRs, which corresponds to the fall-off energy of iron nuclei, is 26 × E_c = 1.2 × 10⁸ 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 ~ 2 × 10⁷ GeV. At 10⁸ GeV, they contribute only ~ 30% 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 3μG, the Larmor radii for cosmic rays of energy Z × 10⁸ GeV is 36 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 Sect. 6, after investigating the effect of different extra-galactic contributions below the ankle in Sect. 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^-1 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, F_esc, can be calculated as, (5)where N(z,p) is given by Eq. (3). Equation (5) assumes that cosmic rays escape only through the diffusion boundaries located at z = ± L. Under this assumption, the total escape rate of SNR-CRs is given by, (6)where A_esc = πR² 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 2 is to account for the two boundaries at z = ± L. The propagation of the escaped SNR-CRs in the Galactic wind region is governed by the following transport equation: (7)where we have assumed a spherically symmetric geometry characterised by the radial variable r, D_w represents the diffusion coefficient of cosmic rays in the wind region which is taken to be spatially constant, N_w(r,p) is the cosmic-ray number density, is the wind velocity which is assumed to increase linearly with r and directed radially outwards, is a constant that denotes the velocity gradient, and Q_esc(p) is given by Eq. (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 100 kpc (Zirakashvili et al. 1996). The second term on the left-hand side of Eq. (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 Eq. (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 Q_esc to be a point source located at r = 0. By solving Eq. (7) analytically, the density of cosmic rays at distance r is given by (see Appendix A), (8)From Eq. (8), the cosmic-ray flux with momentum/nucleon p at the termination shock is obtained as, (9)where R_sh represents the radius of the termination shock. The total rate of cosmic rays injected into the termination shock is given by, (10)where is the surface area of the termination shock. Assuming that only a certain fraction, k_sh, 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), (11)where we have introduced an exponential cut-off in the spectrum at momentum Zp_sh with p_sh representing the maximum momentum for protons, and γ is the spectral index. In our calculation, p_sh and k_sh 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 Eq. (7). In the absence of adiabatic losses, the density of re-accelerated cosmic rays at the Earth (taken to be at r = 0) is given by, (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 D_w to follow the same rigidity dependence as D, and take D_w = 10D. For the wind velocity, we take the velocity gradient km s^-1 kpc^-1. 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 D_w 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 Eq. (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, P_IGM, 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, (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 Sect. 2.2, the total power invested in SNR-CRs (which dominates the overall cosmic-ray energy density in our model) is ~ 8% of the mechanical power injected by supernova explosions in the Galaxy. This corresponds to a total power of ~ 8 × 10⁴⁰ erg s^-1 injected into SNR-CRs. Using this, and taking an intergalactic pressure of P_IGM = 10^-15 ergs cm^-3 (Breitschwerdt et al. 1991), we obtain R_sh = 96 kpc from Eq. (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 Eq. (12) are shown in Fig. 3. Spectra for the individual elements and also the total contribution are shown. The same particle injection fraction of k_sh = 14.5% is applied to all the elements, and the maximum proton energy corresponding to p_sh is taken as E_sh = 9.5 × 10⁷ GeV. These values are chosen so that the total GW-CR spectrum reasonably agrees with the observed all-particle spectrum between ~ 10⁸ and 10⁹ 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 Eq. (12) will lead to further suppression of the flux at low energies. But, at energies of our interest, that is above ~ 10⁷ 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 100 G, and it has been argued that a strong supernova shock in such a field can lead to particle acceleration of energies up to ~ 3 × 10⁹ 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 Eq. (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 Q(p) follows Eq. (2). We assume that each Wolf-Rayet supernova explosion releases a kinetic energy of 10⁵¹ ergs, same as the normal supernova explosion in the interstellar medium. From the estimated total number of Wolf-Rayet stars of ~ 1200 in the Galaxy and an average lifetime of ~ 0.25 Myr for these stars (Rosslowe & Crowther 2015), we estimate a frequency of ~ 1 Wolf-Rayet explosion in every 210 years. This corresponds to ~ 1 Wolf-Rayet explosion in every 7 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 Fig. 4. The results are for two different compositions of the Wolf-Rayet winds available in the literature: carbon-to-helium (C/He) ratio of 0.1 (top panel) and 0.4 (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 ~ 10⁸ and 10⁹ GeV. For C/He = 0.1, we obtain an injection energy of 1.3 × 10⁴⁹ ergs into helium nuclei from a single supernova explosion and a proton source spectrum cut-off of 1.8 × 10⁸ GeV, while for C/He = 0.4, we obtain 9.4 × 10⁴⁸ erg and 1.3 × 10⁸ GeV respectively. For both the progenitor wind compositions, the total amount of energy injected into cosmic rays by a single supernova explosion is approximately 5 times less than the total energy injected into SNR-CRs by a supernova explosion in the Galaxy. The total WR-CR spectrum for the C/He = 0.1 case is dominated by helium nuclei up to ~ 10⁹ GeV, while for the C/He = 0.4 case, helium nuclei dominate up to ~ 2 × 10⁸ GeV. At higher energies, carbon nuclei dominate. One major difference of the WR-CR spectra from the GW-CR spectrum (Fig. 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 ~ 10⁵ GeV. Below the knee, the total WR-CR spectrum is an order of magnitude less than the total SNR-CRs spectrum (Fig. 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 Fig. 5. For the SNR-CRs shown in the figure, we have slightly reduced the value of E_c from 4.5 × 10⁶ GeV (as used in Fig. 2) to 3 × 10⁶ 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 E^-2 and an exponential cut-off at 10¹¹ 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 Fig. 6. For the SNR-CRs, here we take E_c = 4.1 × 10⁶ GeV, and a slightly lower value of ν which corresponds to 6 out of every 7 supernova explosions in the Galaxy (assuming a fraction 1/7 going into Wolf-Rayet supernova explosions as deduced in the previous section). The injection energy f 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 f values are listed in Table 3. The cosmic-ray propagation parameters are the same as in Fig. 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 Fig. 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 (C/He = 0.1) scenario, where the fraction reaches up to ~ 63% at energy ~ 2 × 10⁸ GeV. In contrast to common perceptions, the WR-CR scenarios predict a composition of Galactic cosmic rays dominated mainly by helium (in the C/He = 0.1 case) or carbon nuclei (in the C/He = 0.4) 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 ~ 3 × 10⁵ GeV is not quite as well-measured as at lower energies. Above ~ 10⁶ 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 (X_max). Heavier nuclei interact higher in the atmosphere, resulting in smaller values of X_max as compared to lighter nuclei. For comparison with theoretical predictions, we often use the mean logarithmic mass, ⟨ lnA ⟩, of the measured cosmic rays which can be obtained from the measured X_max values using the relation (Hörandel 2003b), (14)where and represent the average depths of the shower maximum for protons and iron nuclei respectively given by Monte-Carlo simulations, and A_Fe is the mass number of iron nuclei.

In Fig. 8, the ⟨ lnA ⟩ 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 ~ 10⁷ GeV, the GW-CR model deviates from the general trend of the observed composition which reaches a maximum mean mass at ~ 6 × 10⁷ GeV, and becomes gradually lighter up to the ankle. However, in the narrow energy range of ~ (1−5) × 10⁸ 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 (3−5) × 10⁷ 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 X_max 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, 2005; Berezinsky & Grigorieva 1988; Berezinsky et al. 2006; 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 10⁸ 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), (15)where K₀ is a normalisation constant, F_j is the injection fraction which depends on the type of the nuclei j, E₀ = 10⁹ GeV, γ is the source spectral index which is assumed to be the same for the different nuclei, and R_c 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, X_max and variance of X_max above the ankle observed at the Pierre Auger Observatory. We adopt the CTG¹ 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 γ = 0.73, R_c = 3.8 × 10⁹ GV, F_H = 0%, F_He = 0%, F_N = 98.69% and F_Fe = 1.31%. In this model, the EG-CR spectrum below ~ 10¹⁰ 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 ~ 6 × 10¹⁰ GeV, the spectrum is dominated by the CNO group. Above ~ 3 × 10¹⁰ 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 E^-2 source spectrum with a cosmological evolution ∝ (1 + z)^3.5, 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 composition², which is closer to our predicted composition for the WR-CR model (C/He = 0.4) 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, E = ZeBR, if B is the typical magnetic field, and R 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 Fig. 9, we present the all-particle spectrum above 10⁶ 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 (C/He = 0.4). For the SNR-CRs, all the model parameters are the same as in Fig. 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 1.7 × 10⁸ GeV. This is approximately a factor 1.3 larger than the value used in Fig. 6. For the PCS and the UFA models, the proton cut-off energies are almost the same at 1.1 × 10⁸ 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 ~ 7 × 10⁸ GeV. The variation in the injection energy of WR-CRs remain within 6% between the three models. In Fig. 9, spectra of five different mass groups are also shown. The elemental fraction of these mass groups are shown in Fig. 10.

In Fig. 11, we show ⟨ lnA ⟩ predicted by the three EG-CRs model after adding the Galactic contribution. At energies between ~ 3 × 10⁸ GeV and 3 × 10⁹ 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 ~ 10⁹ GeV. The UFA model over predicts the data above the ankle as the model is also tuned to the variance of ⟨ lnA ⟩, 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 10⁹ GeV in the PCS model is due to the dip in the proton spectrum (Fig. 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 ~ 10⁹ 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 Sect. 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 Fig. 10, in which they all have their characteristic “fingerprint”: for example, around 10⁹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 (C/He = 0.1) scenario are given in Appendix B. The main difference from the results of the C/He = 0.4 scenario is the significant dominance of helium up to the transition energy region from Galactic to extra-galactic cosmic rays (see Figs. B.1 and B.2). 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 ~ 10⁹ GeV can be predominantly of Galactic origin. Above 10⁹ 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 ~ 10⁷ GeV, and (ii) GW-CRs or preferably WR-CRs which dominates at higher energies up to ~ 10⁹ 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 (6−8) × 10⁸ 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 10⁹ 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 ~ 4.5 × 10⁶ GeV (see Fig. 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 ~ 10−100 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 ~ 10⁷ GeV in the SNR-CRs, and at ~ 6 × 10⁶ 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 ~ 2 × 10⁷ GeV. SNR-CRs alone cannot explain the observed all-particle spectrum above ~ 2 × 10⁷ GeV. They contribute only ~ 30% of the observed cosmic rays at ~ 10⁸ GeV. On the other hand, in the poly-gonato model, the total contribution from elements with 1 ≤ Z ≤ 28 can explain the observed spectrum up to energies close to 10⁸ 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 ~ 10⁷ 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 X_max 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, E_m ~ 3ZeB(V_s/c)R_s, where B is the magnetic field, V_s is the shock velocity and R_s is the shock radius. From the GW-CR parameters obtained in our study, we can take R_s = R_sh = 96 kpc, km s^-1 which is the terminal wind velocity, and E_m = 9.5 × 10⁷ GeV which is the proton cut-off energy. Using these values, the magnetic field strength in the Galactic halo is estimated to be ~ 73 nG. This is approximately a factor 3 less than the value obtained assuming Parker’s magnetic field topology for the solar wind (Eq. (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 Sect. 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 Fig. 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 ~ 50−100 kpc, can accelerate SNR-CRs up to ~ Z × 10⁸ 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 ~ 10⁸−10⁹ 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 ~ 10⁷−10⁹ 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 ~ 10⁷ 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 10⁸ and 10⁹ GeV by the WR-CR (C/He = 0.1) 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 ~ 10⁸ GeV. In Fig. 12, the elemental fraction predicted by the WR-CR (C/He = 0.1) 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: 1 ≤ A ≤ 6, 7 ≤ A ≤ 19, 20 ≤ A ≤ 39, and 40 ≤ A ≤ 56. 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 R_w from the star follows the relation, (16)where B₀ is the magnetic field at the surface of the star, ω is the angular rotation velocity, R_⋆ is the radius of the star, and V_w is the wind velocity. Using the relation E_m ~ 3ZeB(V_s/c)R_s for the maximum energy as in the case of the GW-CRs, and the proton cut-off energy of E_m = 1.1 × 10⁸ GeV for the WR-CRs (C/He = 0.4) obtained using the PCS/UFA model, we get BR_s ~ 1.2 × 10¹⁵ G cm, where we take the shock velocity V_s = 0.1 c (Soderberg et al. 2012). Using this value of BR_s in Parker’s magnetic field configuration (Eq. (16)) by taking R_w = R_s and other Wolf-Rayet star parameters as R_⋆ = 3 × 10¹² cm, ω = 10^-6 s^-1, and V_w = 2000 km s^-1 (Berezhko & Völk 2000), we obtain the magnetic field at the surface of the star as B₀ ~ 1.5 × 10⁴ 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 100 G in the observable parts of Wolf-Rayet winds, de la Chevrotière et al. (2013) estimated an upper limit for the surface magnetic field of ~ 5400 G. An even stronger field in the wind, up to ~ 2000 G, has been reported (de la Chevrotière et al. 2014), which indicates that the surface magnetic field of these stars can go well above the order of 10⁴ G.

From the total energy of 1.4 × 10⁴⁹ ergs injected into WR-CRs by a single supernova explosion, and the explosion rate of Wolf-Rayet stars in the Galaxy of 1/210 yr^-1, we estimate the total power injected into WR-CRs as 2.1 × 10³⁹ erg s^-1. This is approximately a factor 40 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 1.2−1.6 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 (Eq. (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 Fig. 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 Z × 10⁸ 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 ~ 10⁷ 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 ~ 10⁷ 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 ~ 10⁸ and 10⁹ GeV, the predicted all-particle spectrum in Hillas (2005) consists of a significant iron fraction which may be in agreement with the ⟨ lnA ⟩ data when mixed with a strong extra-galactic proton component, but is in tension with the small iron fraction (~ 2−10%) 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 model that does not require the introduction of an additional Galactic component is to assume that EG-CRs have a significant contribution down to energies below the second knee. Such a scenario would require a steep spectrum of ~ E^-3 and a strong flux suppression below ~ 10⁸ GeV (see Hillas 2005, for a brief discussion; and also Muraishi et al. 2005, in the context of the origin of the knee). To explore this scenario, we inject an additional extra-galactic component of pure protons at the position of the Galactic wind termination shocks, and allow them to propagate diffusively towards the Galactic disk in the presence of the Galactic wind outflow. The injection spectrum is assumed to follow E^{− γ}exp(−E/E_c). The propagation is treated exactly the same as the propagation of GW-CRs from the termination shock towards the Galactic disk. All propagation parameters are kept the same, except for the wind velocity constant which is treated as a free parameter. The best-fit all-particle spectrum obtained after adding the contribution of SNR-CRs and EG-CRs from the minimal model is shown in Fig. 13. The best-fit parameters are γ = 3.3, E_c = 4.1 × 10⁶ GeV for the SNR-CRs protons, E_c = 1.5 × 10⁹ GeV for the additional EG-CRs, and km s^-1 kpc^-1. This value of gives a wind velocity which is about a factor 13 larger than the wind velocity used in the study of GW-CRs. Such a fast wind is required in order to generate a strong modulation for particles below the second knee so that the predicted flux does not exceed the observed data at low energies. For γ< 3.3, the required wind velocity is lower, but the model prediction does not fit the observed data very well (see e.g. the case of γ = 3 in Fig. 13). Replacing the additional extra-galactic protons with heavier elements only slightly reduces the required wind velocity. Having a strong Galactic wind can have serious effects on the spectrum and distribution of low-energy cosmic rays in the Galaxy (see e.g. Bloemen et al. 1993). In the presence of a strong wind, cosmic-ray transport will be dominated by advection rather than diffusion, and will produce a cosmic-ray distribution that resembles the distribution of the sources. But, the cosmic-ray distribution inferred from the observations of diffuse gamma-ray emission from the Galaxy indicates a radial gradient weaker than the distribution of supernova remnants or pulsars in the Galaxy. These observations suggest that if supernova remnants are the main sources of cosmic rays in the Galaxy, the propagation of cosmic rays should be dominated by diffusion, not by advection. In addition, if the transport is dominated by advection, the cosmic-ray spectrum is expected to exhibit a break (steepening) at an energy where the advection boundary, , equals the halo boundary L. Such a break is not observed below the knee, except at ~ 10 GeV which is due to Solar modulation. Attributing the knee to such a break raises issues regarding the cosmic-ray injection index. Below the break, cosmic-ray transport is advection dominated and the spectrum is expected to follow E^{− (γ + a/ 2)}, where γ is the source index and a is the diffusion index. For the observed spectral index of ~ 2.7 and a = 0.33 used in our study, we get γ = 2.53. This is incompatible with the prediction of diffusive shock acceleration theory which predicts an index close to 2 for the strong shocks present in supernova remnants (Ptuskin et al. 2010; Caprioli et al. 2011). Choosing a = 0.6, as in pure diffusion propagation models, gives γ = 2.4. This relaxes the tension a bit, but such a high value of a is not favoured by the observed small level of cosmic-ray anisotropy. Another strong constraint on the Galactic wind velocity is provided by the abundance ratio of radioactive secondary to stable secondary. Measurement of ¹⁰Be/⁹Be ratio puts a constraint at km s^-1 kpc^-1 (Bloemen et al. 1993). All these arguments pose a serious problem to the alternative scenario of a strong EG-CR component with a steep spectrum extending below the second knee, and modulating by Galactic wind. One possibility, but rather unrealistic, for this scenario to work is if the additional EG-CR component has a spectrum and composition almost similar to that of the GW-CRs produced at the Galactic wind termination shocks.

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 (≳ 1 nG) extra-galactic field with a constant coherence length extending over the entire Universe, this effect could set in at around 10⁹ 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 γ = 2.7 as indicated by the KASCADE-Grande measurements above ~ 10⁸ 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 10⁸ GeV.

An additional problem for EG-CRs with an overall spectrum steeper than E^-2.7 is that, if one assumes that they fill the extra-galactic space homogeneously with energies from ~ 1 GeV to 10⁹ 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 10⁸ 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 ~ 10⁶ 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 ~ 5 × 10⁷ GeV and the second component dominating in the range of ~ 5 × 10⁷−10⁹ GeV, can explain almost all observed features in the all-particle spectrum and composition when combined with an extra-galactic component dominating above ~ 10⁹ 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 10⁸ 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 ~ 10⁸ and 10⁹ 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.

---

Acknowledgments

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. S.T. wishes to thank Onno Pols for discussions on Wolf-Rayet stars. A.v.V. 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

Appendix A: Derivation of Eq. (8)

The Green’s function, G(r,r′,p,p′), of Eq. (7) satisfies, (A.1)In rectangular coordinates, the above equation can be written as, (A.2)where we have written with î, ĵ and representing the unit vectors along the x, y and z directions. Following a similar procedure adopted in Lerche & Schlickeiser (1982b), we express, (A.3)and, (A.4)Inserting Eqs. (A.3) and (A.4) into Eq. (A.2), we get, (A.5)We now introduce variables ψ_x, ψ_y and ψ_z such that k_x = ψ_xF(p), k_y = ψ_yF(p) and k_z = ψ_zF(p), where (A.6)This reduces Eq. (A.5) to (A.7)where, (A.8)The solution of Eq. (A.7) is given by, (A.9)where the Heaviside step function H[p−p′] = 1(0) for p>p′( <p′). Taking inverse Fourier transform of , we obtain the required Green’s function as, (A.10)where, (A.11)and, (A.12)Then, for a given cosmic-ray source characterised by q(x′,y′,z′,p′), the differential number density of particles with momentum p at a distance (x,y,z) is given by, (A.13)For any point source located at (0,0,0) and emitting q(p) spectrum of particles, that is q(x′,y′,z′,p′) = δ(x′)δ(y′)δ(z′)q(p′), the solution becomes, (A.14)From Eq. (A.6), since F(p) reduces to p, and so also F(u) to u, by writing (x² + y² + z²) = r² in spherical coordinates and replacing q(p) by Q_esc(p) as given by Eqs. (6), (A.14) can be reduced in the form of Eq. (8): (A.15)

Appendix B: All-particle spectrum and composition of cosmic rays obtained using different EG-CR models and WR-CRs (C/He = 0.1)

The predicted all-particle spectrum, elemental fraction and ⟨ lnA ⟩ obtained for the three different models of EG-CRs (the minimal, PCS and UFA), combined with the WR-CR (C/He = 0.1) scenario for the additional Galactic component, are shown in Figs. B.1−B.3, respectively. The proton cut-off energies for the WR-CRs required to produce a good-fit to the measured spectrum are 2.4 × 10⁸ GeV for the minimal model, 1.5 × 10⁸ GeV for the PCS model, and 1.6 × 10⁸ GeV for the UFA model.These values are about a factor 1.4 larger than the cut-off energies obtained in the case of C/He = 0.4. The variation in the injection energy of WR-CRs between the three cases remain within 6% as in the C/He = 0.4 scenario.

The predicted composition is dominated by helium nuclei up to around the second knee for the minimal and the UFA models, while for the PCS model, helium dominates up to around 10¹⁰ GeV. The Galactic component at the transition energy region from Galactic to extra-galactic cosmic rays is dominated by helium, unlike in the case of C/He = 0.4, where it is dominated by a mixture of helium and CNO group.

The ⟨ lnA ⟩ predicted by the minimal model shows some deviation from the general trend of the measurements between ~ 10⁸ and 5 × 10⁹ GeV, although the discrepancy is less than that observed in the C/He = 0.4 scenario. The predictions of both the PCS and the UFA models show better agreement with the data below ~ 10⁹ GeV. Between around 10⁷ and 10⁹ GeV, they predict a mean mass lighter than the prediction of the C/He = 0.4 case, and show a better agreement with the data (EPOS-LHC) from the Pierre Auger Observatory in the 10⁸−10⁹ GeV energy range, but slightly under predict the available measurements at around ~ 10⁸ GeV. The two WR-CR scenarios should be possible to differentiate by accurate measurements of the elemental composition between 10⁷ and 10⁹ GeV.

All Tables

All Figures