On the cosmic ray spectrum from type II Supernovae expanding in their red giant presupernova wind

While from the energetic point of view SNRs are viable sources of Galactic CRs, the issue of whether they can accelerate protons up to PeV remains unsolved. Here we discuss particle acceleration at the forward shock of SN and discuss the possibility that the escaping particle current may excite a non-resonant instability that in turn leads to the formation of resonant modes confining particles close to the shock and increasing the maximum energy. This mechanism works throughout the expansion of the SN explosion, from the ejecta dominated (ED) to the Sedov-Taylor (ST) phase. Because of their higher explosion rate,we focus on type II SNae expanding in the slow, dense red supergiant wind. When the explosion occurs in such winds, the transition between the ED and the ST phase is likely to take place within a few tens of years. As a result, the spectrum of accelerated particles shows a break in the slope, at the maximum energy (Em) achieved at the beginning of the ST phase. Above this energy, the spectrum becomes steeper but remains a power law than developing an exponential cutoff. We show that for type II SNae typical parameters, proton Em can easily reach PeV energies, confirming that type II SNRs are the best candidate sources for CRs at the knee. We have tried to fit KASCADE-Grande, ARGO -YBJ and YAC1-Tibet Array data with our model but we could not find any parameter combination that could explain all data sets. Indeed the recent measurement of the proton and helium spectra in the knee region, with the ARGO-YBJ and YAC1-Tibet Array, has made the situation very confused. These measurements suggest that the knee in the light component is at 650 TeV, appreciably below the overall spectrum knee. This finding would resolve the problem of reaching very high energies in SNae, but, on the other hand, it would open a critical issue in the transition region between Galactic and extragalactic CRs.


Introduction
Supernova remnants (SNRs) can provide the correct order of magnitude of energy injection in the form of accelerated particles to explain the density of CRs observed in the Galaxy (for reviews see, e.g., [18,8]). However, the issue of whether particles can be accelerated diffusively at the SN shock to sufficiently high energies is still open for debate. One thing that is clear is that effective magnetic field amplification is necessary, since in the absence of this phenomenon the maximum energy is bound to be in the ∼ GeV region. This need was already recognized by [39,40], who found that the resonant instability induced by CR streaming could amplify the magnetic field to δB/B ∼ 1 and allow the maximum energy to reach ∼ 10 3 − 10 4 GeV. This latter upper bound on the maximum energy, still well below the knee, derives from the limit δB/B 1, which is in turn related to the resonant nature of the instability considered [38,63,54,55,56]. Non resonant instabilities may in fact potentially lead to larger amplification factors, but on scales either much larger or much smaller than the particle gyroradius, so that larger fields would not necessarily reflect into increased particle scattering, which is what is needed to reduce the acceleration time and increase the maximum achievable energy.
The pioneering papers of [41] and [12] showed how non-resonant modes can grow very fast ahead of a SNR shock. The fastest growing mode is characterized by wavenumber k much larger than the Larmor radius of the particles contributing the CR current, but the non-linear evolution of these modes is shown to form fluctuations on larger scales, up to the Larmor radius of the particles constituting the CR current. At that point, resonant scattering may in principle facilitate the return of particles to the shock surface and possibly increase the maximum energy. It has been argued that this non-resonant branch of the streaming instability may account for very high energy CRs in the Galaxy [13,52] and, in particular, may allow to reach the knee if applied to the case of supernovae expanding in the wind of their red super-giant (RSG) progenitor star. The role of winds of evolved pre-supernova stars for CR acceleration has been discussed several times in the literature, with different spins. Since massive stars often go through a Wolf-Rayet (WR) evolutionary phase and they are usually clustered within OB associations (superbubbles), several authors (see [15] for a recent review) have proposed that the collective effect of the WR fast rarefied winds may help accelerating CRs to very high energies and at the same time provide an explanation for some anomalies in the CR composition [33]. The most notable of such anomalies, and the one that is most solidly characterized from the observational point of view, is the ratio of abundances of two Ne isotopes, 22 N e/ 20 N e [34], as observed in CRs, in comparison with its value in the solar environment. Additional constraints come from observations of 59 N i in CRs, which impose limits on the rate of SN explosions in superbubbles, or, to state it more accurately, on the time when CR acceleration starts as compared with the time between two SN explosions in the same OB association. While this latter constraint is not a serious challenge for current models [15], it was recently pointed out [44] that the average 22 N e/ 20 N e in superbubbles must approach the solar value, which would kill one of the main arguments in support of super bubbles as the main GCR sources.
In this paper we focus on the explosion of SNe in the wind of the RSG progenitor star and we argue that these explosions are the best sites to reach PeV energies. A number of SNe expanding in their RSG wind have been detected, among which the best studied case is that of SN 1993J [51]. This source is especially interesting for our purposes because radio measurements provide us with direct estimates of the magnetic field strength at different times, through the detection of the synchrotron self-absorption feature [32]. This suggests a scenario in which not only the field is amplified at the SN shock, but the amplification turns out to be consistent with the expectations of the non-resonant instability [58]. Because of the high density of such winds, the highest energies are reached within a few decades after the SN explosion, namely before the beginning of the ST phase of the explosion. In this scenario, particles are accelerated from a pool that is mainly made of wind particles, rather than ejecta of other SN explosions. Eventually one may expect that the shock reaches the edge of the bubble excavated by the stellar wind, and at that point particle acceleration proceeds in the ordinary way, although the details will depend on whether the star was originally in an OB association or in the ISM. At such late times, the particles escaping the SN are less energetic and we will not be concerned much with this stage of the acceleration process. On other hand, it is worth recalling that most constraints on composition and composition anomalies are limited to such lower energies.
The origin of the knee in the all-particle spectrum of CRs is inextricably connected with the issue of the end of the Galactic CR spectrum and the transition to extragalactic CRs: in Ref. [6], the authors try to model the transition region by summing the Galactic SNR contribution and the flux of nuclei of extragalactic origin, required to fit the Auger data [4]. The authors reach the conclusion that in order to satisfy the observational requirement of a predominantly light composition in the EeV region [4,2,3,60], an extra component of extragalactic protons is required. Such an extra component appears to be in good agreement with the proton spectrum as measured by KASCADE-Grande [10]. However, both the proton and iron spectra measured by KASCADE-Grande in the energy region 10 16 −10 18 eV suggest that either the injection spectra are not cut off exponentially at the maximum energy (namely the number of particles decreases more slowly than an exponential at energies E ≫ EM ), or there is some, as yet unknown, class of sources with maximum energy much in excess of the knee energy.
The observational situation has recently become even more confusing after the publication of the spectra of light CRs (protons and helium) by the ARGO-YBJ [26] and YAC1-Tibet Array [36] experiments. Both measurements suggest a knee in the light component at ∼ 650 TeV, appreciably below the knee. On one hand, it is clear that this finding would make the requirements in terms of particle acceleration less severe. On the other, however, it makes the requirements on the transition region much harder to accommodate. In other words, even in the scenario arising from ARGO data a population of sources able to reach maximum energy much above the knee is required to exist, thereby leaving the problem of the maximum energy little affected.
In this paper we consider the viability of the two hypotheses above: we investigate whether non-exponential tails at E > EM are to be expected in the spectra of particles escaping supernova shocks, and we study the possibility of having supernovae that contribute maximum (proton) energies much higher than the knee. The paper is organised as follows: in § 2 we discuss the role of escaping particles for the estimate of the maximum energy achievable at a supernova shock. In § 3 we illustrate our calculations of the spectra of nuclei and in § 4 we compare our predictions with the observed spectra, with particular emphasis on the knee and the transition region. In § 4.1 we discuss the implications of the recent measurements of the ARGO-YBJ and YAC1-Tibet Array, comparing them with KASCADE-Grande data. We conclude in § 5.

Non-resonant instability excited by escaping CRs
CR streaming leads to the excitation of both resonant Alfvén waves and non-resonant, almost-purely-growing, modes [12]. The former have a growth rate that, for typical parameters of a supernova remnant, leads to maximum energy ≪ P eV , mainly as a consequence of the saturation at δB/B ∼ 1. On the other hand, the non-resonant instability is very fast growing, but the fastest growing mode has a wavelength much shorter than the Larmor radius of the particles generating the current, and also the wrong polarisation for resonant interaction with positively charged particles moving away from the shock surface [7]. As a result, to first order, the non-resonant waves do not lead to efficient particle scattering, and would therefore be useless to the goal of increasing the maximum achievable energy.
However, it was already suggested by [12], and then confirmed by the numerical work of [50] and more recently [22], that the non-linear evolution of these modes leads to the generation of power on larger spatial scales. Here it is useful to think in physical terms about the process that may be responsible for particle return to the shock from the upstream region. Let us consider a first generation of particles of energy equal to the current achievable maximum, E, moving from downstream to upstream. By definition, upstream there are no waves able to scatter resonantly such particles, hence they are going to move quasi-ballistically and escape the system at the speed of light (or fractions of it, depending on the level of anisotropy of the distribution of escaping particles). The current of particles at distance R from the explosion center as due to particles escaped when the shock location was at r can be written as: In the above expressions we have used the relation between the number density of accelerated particles with energy E, nCR(E), at some location R upstream of the shock, and the energy density in accelerated particles at the shock ξCRρ(r)v 2 s (r), where ξCR is the CR acceleration efficiency. This relation is easily found through the transport equation for accelerated particles at a plane shock. The function Ψ(EM ) accounts for normalisation of the particle distribution function which is taken to be fs(E) ∝ E −(2+β) and extending between a minimum energy E0, that does not depend on time and a maximum energy EM which does depend on time, as we will see. In Eq. 1, the factor (r/R) 2 should account for the dilution due to spherical propagation of particles in ballistic motion.
The underlying assumption is that the large scale magnetic field is radial or absent. Whether this assumption is a realistic one in the two environments we consider in the following, namely the ISM or the wind of a massive progenitor star, may be questionable. In the latter case the zeroth order expectation is that the field should be toroidal, although one can envision that instabilities induced by the interaction between the slow and dense wind produced during the red supergiant phase of the progenitor and the faster and more dilute wind later blown as a Wolf Rayet star, could stretch the field lines and make them radial in at least a portion of the shock surface. In the case of a SNR expanding in the ISM, more caution should be taken, in that there is no reason why the upstream magnetic field should display a radial geometry.
The fastest growing mode has a growth rate where vA is the Alfvén speed in the unperturbed magnetic field B0. The wavenumber where the growth is the fastest can also be easily estimated using the condition which corresponds to balance between current and magnetic tension. The above expression implies that: where rL is the particle gyroradius, ξCR is the CR efficiency and VA = B 0 √ 4πρ the Alfvén velocity. Upon saturation of the instability, the energy density in the form of amplified magnetic field in units of the ram pressure of the plasma can be estimated as (v sh /c)(ξCR/ ln(EM /E0)), that for typical values of the parameters is ∼ 10 −4 (for simplicity here we assumed that the spectrum of accelerated particles is ∼ E −2 ). This reflects in about one order of magnitude larger pressure downstream of the shock (as due to compression), an estimate which is still a factor of 10 below that of [48], where the magnetic pressure downstream was assumed to be a fixed fraction (∼ 3.5%) of the ram pressure, independent of the shock speed.
Using the expression for kM derived in Eq. 4, we can rewrite this condition, at a point R upstream of the shock, as [13,52,53,19]: where we have used the shock velocity to change the integration variable from time to shock coordinate. In Eq. 7, the integral in r is assumed to extend over all previous radii of the shock. In fact, since the values of the growth rate that are found for typical values of SN parameters are rather large (γM ≫ v sh /R), in case of ballistic motion of the escaping particles, one would expect that only a narrow range of radii close to R would give contribution to the particle current. Here we keep the form of Eq. 7 as given in Refs. [13,52,53], since the difference that this correction would introduce is only a factor of ∼ 2.
In order to find an expression for the evolution of the maximum energy with time, we need to assume a density profile of the medium in which the supernova remnant is expanding. For a supernova expanding in the uniform ISM, one can assume that the ISM density is constant ρ(R) = ρISM . For a type II supernova, instead, it is often the case that the explosion takes place in the wind produced by a red giant pre-supernova progenitor star. The density of the wind can be written as whereṀ is the rate of mass loss of the red giant and Vw is the wind velocity. In the following we will simply write ρ(R) ∝ R −m , having in mind the two cases of m = 0, corresponding to a uniform medium, and m = 2, appropriate for expansion in the progenitor wind. Differentiating Eq. 7 with respect to R, we obtain an implicit expression for the maximum energy (see Fig. 1): 3. From escape of accelerated particles to CRs

Spectrum of the escaping particles
For a supernova shock moving with velocity v sh (t) in the surrounding medium (ISM or wind of the progenitor star) the shock radius is R sh (t) = t 0 dt ′ v sh (t ′ ). At each time t we assume that particles with given energy EM (t) ≡ E can escape the accelerator. The number of particles that escape the shock at each give time is related to the current at the upstream radius R through: where Nesc (E) is the number of particles per unit energy, with energy E = EM (t), that are released in the ISM. Using the expression in Eq. 1 for the current one can rewrite: with dR/dE = (dR/dΨ)(dΨ/dE) to be computed using Eq. 2 and 9. In order to derive the first term on the right in the latter equation, we need to make explicit assumptions on the shock dynamics. For purposes of illustration we assume in the following a general power-law dependence of the shock radius on time: R ∝ t λ , which also implies vs ∝ R (λ−1)/λ . Using the latter dependencies in Eq. 9 it is straightforward to find for dR/dΨ: and rewrite Eq. 11 in a particularly useful way: where ) is a function of E alone, to be computed from Eq. 2. The advantage of Eq. 13 is that it makes it immediately clear that in the Sedov-Taylor phase, during which ρv 2 s R 3 is constant with time, the CR spectrum released in the ISM is bound to be just given by χ(E), and hence only depend on the spectrum in the remnant, unless ξCR or E0 depend on time. It should also be noticed, however, that this dependence is not just a one to one identity and the spectrum of escaped CRs will not be the same as the spectrum of accelerated particles in the remnant, as already noted by [13], [22]. In fact, when computing χ(E) explicitly from Eq. 2 we find: where the sign ≈ refers to the fact that we have used the assumption E >> E0. It is apparent that the power-law dependence of χ(E) can only be −2 or steeper. Therefore the spectrum of CRs released during the Sedov-Taylor phase is the same as the spectrum of accelerated particles in the source only if the latter is E −2 or steeper, while for flatter source spectra Nesc(E) ∝ E −2 . While in the Sedov-Taylor phase, only the environmental details are important while the ejecta density profile is irrelevant to the spectral slope, during the free expansion phase both density profiles are important in order to determine the position of the breaks and the maximum energy to which the released spectrum extends. This kind of information enters the calculation of the dependence of the shock radius on time, which we expressed as In the following, we take the value of λ appropriate to describe the different evolutionary stages of the SNR from Ref. [24,45,59]. This work gives the remnant expansion law during the different stages for generic density profiles of both the SN ejecta ρej ∝ R −k and the ambient medium ρ ∝ R −m . The general expression for λ is: where the ejecta profile appropriate to describe the two types of SN has k = 7 and k = 9 for type I and II explosions respectively. Once λ is given, from Eq. 13 and 14 we can derive the spectrum released by the SN during the ED and ST phases as: Inverting the time dependence of E as can be found from Eq. 9, and expressing the dynamical quantities as functions of E, one finally finds: This result is very interesting because it shows that the spectrum of escaping CRs integrated over time during the SN shock expansion is a broken power law, with an index close to 2 at low energies (belowEM ) and steeper at higher energies. In more general terms, the spectrum of CRs released into the ISM is different from the postshock spectrum of accelerated particles, and is rather sensitive to how magnetic field amplification evolves in time [20,21]. A similar conclusion was previously obtained in the context of simple box models of CR acceleration at shocks [30,31]. At some energy ≫ EM , the spectrum will eventually suffer an exponential cutoff. However such cutoff does not have an immediate phenomenological impact, because it would appear at energies at which the spectrum has already dropped appreciably due to the steeper power law. The transition between the two power laws at the maximum energy EM corresponds to the transition between the end of the ED phase, where most mass has already been processed, and the beginning of the ST (adiabatic) phase. Notice that in principle the highest value reached by particles at the shock is larger than EM at earlier times, but the flux of CRs with such energies is suppressed (though not exponentially but as a power law) because of the small amount of mass that is processed during the ED phase. In order to provide an estimate of the maximum achievable energy as a function of time, we proceed with modeling the transition from the ED to the ST phase of the SN expansion by using a parametric form for the time dependence of the shock radius: where a is a smoothing parameter (in the following we use a = −5 unless indicated otherwise), λED and λST are the indices describing the asymptotic time dependencies during the ED and ST phases, respectively, and R0 and t0 are the values of radius and time at the transition between the two phases. The value of R0 can be estimated by using the condition that at such distance the swept up mass equals that of the ejected material Mej : where ρ(R) is the ISM density in the case of type I SNe and the density in the wind in the case of type II explosions. For the former type we then find, for the transition radius: where ρ is the constant ISM density, which we take to correspond to 1 particle cm −3 . For type II events, instead: The assumption of self-similarity for the time evolution of the SN shock leads to a generic form for the density evolution of the ejecta [59] ρej( where and ESN is the SN explosion energy. The condition that the density of the ejecta and the density of the wind be equal at the forward shock leads to the following expression for the time t0: where B = ρISM for type I and B =Ṁ /(4πVw) for type II respectively. The shock velocity can be easily written as which can be used to calculate the maximum energy as a function of time, and hence the spectrum of escaping particles from an individual supernova. Before proceeding with detailed calculations, let us briefly discuss our choice of focusing on type II SNe as the best candidate producers of CRs at the knee. Let us consider the maximum energy that can be reached by each type of SN events, as given by Eq. 9 evaluated at the transition between the ED and the ST phase. To the goal of a simple estimate, let us specialise this equation to the case of a E −2 source spectrum.
For type I events we can estimate: while for type II: It is apparent that for standard parameters pertaining the two types of explosions, type II SNe can reach about one order of magnitude larger maximum energies than type Ia, and in particular that if particles are accelerated and escape as described here, they easily seem to be able to provide proton acceleration up to the knee.
In order to illustrate how the maximum energy depends on the assumed spectral index within this framework, in Fig. 1 we plot the evolution with time of the instantaneous maximum energy a type II SNR can achieve for three different values of the source spectral index (all the other parameters have standard values, as specified in the caption). The figure shows that EM as derived from Eq. 9 and 2 changes by about 1 order of magnitude at the transition between ejecta dominated and Sedov phase (indicated by the vertical line in the plot) depending on the source spectral index: spectral indices steeper than 2 lead to lower values of EM compared to the estimate in Eq. 28, while the opposite is true for flatter spectra.

CR Spectra observed at the Earth
Here we adopt a simple model for the description of the transport of nuclei in the Galaxy. Assuming that the transport has reached a stationary regime and that it is dominated by diffusion, the equation describing the transport of a nuclear specie α is: where Qα is the rate of injection per unit volume of CR of type α and σsp,α is the cross section for spallation. In principle one should include an injection term of particles of type α due to the spallation of heavier nuclei, but as long as we focus on primary species and we are at sufficiently high energies, this contribution can be neglected [14]. Eq. 29 can be easily solved if we assume that both injection and spallation occur in a thin disc of size 2h: Here the injection spectra Nα(p) are calculated as discussed in Sec. 3.1, n d is the gas density of the Galactic disc, assumed of thickness 2h and radius RD. If we introduce the grammage and we impose a free escape boundary condition at |z| = H, the solution of Eq. 29 in the Galactic disc can be easily written as follows: where µ = 2hn d mp is the surface density of the Galactic disc and Xα = mp/σsp,α. Written in this form, the solution is very clear. At particle momenta for which the grammage is small compared with Xα the standard solution is recovered: while in the range of momenta at which the solution is spallation-dominated, the observed spectrum reproduces the shape of the injection spectrum. In order to calculate the spectra of CRs at the Earth, an evaluation of the grammage is needed. In Ref. [47], the authors develop a leaky-box model in order to reproduce the GALPROP results and, for a source spectrum E −2 find a grammage in the form: where R is the rigidity of particles and β = v/c ∼ 1 for 1 TeV particles. They calculate also the grammage considering the reacceleration and find: Given that the CR spectrum observed at Earth has an energy dependence E −2.65 as taken from TRACER and CREAM data [11,61], in our calculations we used a slope of the diffusion coefficient δ = 2.65 − pinj , where for pinj we considered the values 2 and 2.31 as we will discuss in Sec 4. In the first case, the diffusion coefficient has a slightly stronger energy dependence than the above mentioned E 0.6 ; in the other case, instead, we use the exact form of Eq. 35. For the spallation cross section we adopt the simple formulation of Ref. [35]: where A is the atomic mass number, and the two energy dependent coefficients are It is worth noticing that one of the values of δ we adopt here, and in fact the one that leads to most interesting results, is rather large, δ = 0.65 (corresponding to pinj = 2), and this is known to lead to exceedingly large anisotropy in the CR signal at Earth [46,17].

Results
Having shown that type II SNe are the most likely sources to accelerate protons up to the knee, in this section we focus on this type of explosions and attempt at fitting the all particle spectrum. To this aim we vary the different parameters, including the ejecta density profile.
As discussed in § 3, the maximum energy up to which CRs can be accelerated in a SN is regulated by the escape of particles from the shock region and in the very early stages of the SN evolution, in principle, this energy can exceed the knee energy. However, the amount of mass processed at such times is relatively small and the net flux of particles is correspondingly small. The spectrum of accelerated particles that escape the shock region during the Sedov-Taylor phase and become CRs reflects the spectrum in the source if fs ∝ E −(2+β) with β ≥ 0, and is ∝ E −2 if β ≤ 0. This spectral slope extends up to what we can name the effective maximum energy, EM , reached at the end of the ejecta dominated phase of expansion. We find that, at energies larger than EM , the spectrum is not falling exponentially, but rather as a power-law with a steeper index as compared to the lower energy trend. After escape from individual sources, CRs propagate diffusively through the Galaxy. At the high energies we are interested in here, the only relevant process during propagation is spallation, while other processes, such as reacceleration or Galactic winds, induce negligible effects and are therefore ignored in the present calculation. It is also worth stressing that the particles that are actually accelerated in the early phases of the SN explosion (say, within the first 100 years) are the highest energy particles; consequently, no direct constraint on these accelerators can be inferred from 59 N i, which is often used to impose a lower limit on the time lag between the explosion and the acceleration phase. In any case, in this scenario, the material from which particles are accelerated to the highest energies is not polluted with SN ejecta, since it is basically the material expelled during the RSG phase of the pre-supernova star.
As discussed in § 3, the value of EM is determined by the combination of the CR acceleration efficiency, ξCR, and the energetics of the SN, ESN . On the other hand, the flux of CRs at the Earth derives from a combination of ξCR, ESN and the rate of SNe, ℜ. The shaded (yellow) region in Fig. 2 shows the region of the ξCR − ESN plane for which EM ≥ 1 PeV. The two panels refer to two different spectra at the source: E −2 is on the left and E −2.31 on the right.
The need to fit the proton spectrum in the knee region imposes an additional constraint on the combination ξCR ESN ℜ. In Fig. 2 the lines indicate such a constraint projected on the ξCR − ESN plane for SN rates as specified in each panel.
One can see that for the case of a "flat" source spectrum E −(2+β) with β ≤ 0, the standard SN energetics of 10 51 erg is compatible with the need to have maximum energy in the PeV range with totally reasonable values of the CR acceleration efficiency (ξCR ∼ 10%) and SN rate (ℜ ∼ 1/30 yr −1 ). If the rate is decreased to 1/100 yr −1 , higher energetics and higher efficiencies are required to reproduce the CR flux at Earth and still accelerate particles up to 1 PeV, whereas for a SN rate ℜ ≥ 1/30 yr −1 , the knee simply cannot be reached: this is because for large SN rates the allowed CR energetics is low in order to avoid overproducing the flux at Earth. The conclusion is even more severe for SNRs with a steeper internal spectrum: in this case, the sources that can reach the knee and fit the spectrum are extremely rare, with a frequency less than 1/800 yr −1 ; in addition the source energetics is required to exceed a few ×10 51 erg and the acceleration efficiency several tens of percent.
The conclusion is that the most common type II SNe, expanding in their red supergiant wind, with standard energetics and an efficiency of order 10% are the most likely contributors of CRs up to the knee. In Fig. 3   show the spectra of the different nuclei, as well as the all-particle spectrum, for a population of type II SNe with a density profile of the ejecta described by k = 9, an energy release corresponding to ESN = 10 51 erg, exploding at a rate ℜ = 1/30 yr −1 . The spectra of nuclei are normalised to the data of CREAM for protons [61], PAMELA for Helium [5] and TRACER for CNO and Iron [11]. The diffusion coefficient is taken as described in Sec. 3.2. Clearly, in a simple approach such as the one discussed here, it is not possible to account for subtleties such as the harder spectrum of helium compared to hydrogen. We simply normalize the abundances in such a way that the CR composition at Earth is reproduced. We consider this assumption as a weak point of all models trying to explain the origin of CRs: there is still no comprehensive theory of acceleration of nuclei in SN shocks. In fact, it is known, and can be easily understood, that there is a preferential injection of heavier nuclei at a collisionless shock [27], but a quantitative theory of this phenomenon is not available yet. The situation is made even more complex by the fact that most refractory elements are embedded in dust grains, so that dust sputtering at shocks is bound to have en important role in determining chemical composition of CRs observed at the Earth. At present, the last and most comprehensive work on nuclei injection from dust sputtering is the one of Refs. [43,29], but a physical understanding, from first principles, of injection of particles with different masses and charges, is still lacking, so that one is forced to parametrize the injection of nuclei: this is what is commonly done and what we have done in this work, fixing the abundances so as to fit observations. Fig. 3 shows that for the values of the parameters that are typical of type II SNe the maximum energy of protons is in the PeV range, close to the knee at 3 × 10 15 eV. The maximum energy of heavier nuclei is Z times larger. The spectrum of protons (as well as of other nuclei) is not cut off sharply at the highest achievable energy, but rather shows a severe steepening to a different power law. The slope of the high energy part is determined by the profile of the ejecta. For k = 9 the spectrum at the source is steepened from E −2 to E −4 . The low energy part of the spectra of nuclei heavier than He shows a hardening that is due to spallation during propagation. The superposition of the spectra of all nuclei returns the all particle spectrum (black curve) that shows clear evidence for a knee, at energy of ∼ 1 PeV. The high energy hardening of the protons spectrum is due to the introduction of an ad hoc extragalactic proton component with spectrum ∼ E −2.7 , as found by the KASCADE-Grande experiment.
In Fig. 4 we plot the maximum energy of accelerated particles as a function of time for an acceleration efficiency of 10% and different energetics of the parent SN, as labeled. One can see that more energetic SNe accelerate particles to higher energies and the onset of the Sedov phase in the wind occurs at earlier times. In any case, one should appreciate how in the framework of particle acceleration in SNe that explode in the wind of the presupernova star, the most important phase in the acceleration process occurs a few tens of years after the explosion, thereby implying a change of paradigm in which SNe, rather than SNRs, play a central role. The fact that in principle increasingly larger energies are reached at earlier times raises the issue of what is the minimum time that we need to take into account in the calculation of the spectra of escaping particles. One obvious limitation comes from the fact that, after escape, particles need to cross the wind thereby suffering different kinds of energy losses. For protons, inelastic pp scattering occurs with a loss time scale τpp = 1/nσc, and it is unimportant when where reference values of the parameters have been adopted. For a shock velocity of 10 4 km/s, this corresponds to times t ≫ 10 3 s. Given the large photon background in the early phases of the SN explosion, it is worth considering  also inverse Compton scattering losses of protons, with rate of energy change: where we approximate the photon energy density as U ph = ESN /( 4π 3 R 3 ). By requiring that E/(dE/dt) ≫ r/c one gets that t ≫ 6.3 × 10 −6 E 1/2 GeV yrs. For protons with PeV energy this constrains times useful for escape to later than the first ∼ 2 × 10 5 s. For more energetic SNe (namely larger values of the shock velocity) this bound becomes less constraining.
Photodisintegration of nuclei in the dense photon background may also become important. However, this is a threshold process, which is turned on when ǫγA > E bind ∼ 10 MeV, E bind being the binding energy of a nucleon inside the nucleus. This leads to a lower bound on the nucleus Lorentz factor γA > 10 7 . Above this threshold, the time scale for the process to become irrelevant, assuming a typical photon energy in the SN of ǫ ∼ 1 eV, is 1 The bound becomes less severe for more energetic SNe, namely for larger shock velocities. We stress once more that this bound is of some concern only for nuclei with rigidity quite above the rigidity of the knee.

Comparison with high energy data
Any understanding of the origin of CRs at the knee is bound to have important implications for higher energy CRs, and more specifically for the transition from Galactic to extragalactic CRs. In this section we discuss this issue in some detail. The standard lore about the origin of the knee in the all-particle spectrum is based on the premise that this feature coincides with a change of mass composition of CRs from light to heavy. In this perspective, the knee is associated with the end of the spectrum of the light component of CRs, say H and He. The natural implication of this line of thought is that the sources of Galactic CRs must be able to accelerate particles to at least the energy of the knee in the all-particle spectrum, say ∼ 3 PeV. The important conclusion has stimulated endless discussion for the last thirty years on the sources and on the type of acceleration mechanism: in 1983, two seminal papers [39,40] reached the conclusion that even in the presence of resonant growth of Alfvén waves the maximum energy reachable in SNRs could not exceed ∼ 10 3 − 10 5 GeV, thereby falling short of the knee by more than one order of magnitude. Measurement of the proton and He spectrum by KASCADE [9] seemed to confirm that the spectrum of the light CR component is steepening in the PeV region. A note of caution should be added, in that the information about mass composition is accessible to KASCADE only through a complex Monte Carlo dependent procedure. For instance using Sybill and QGSJet to simulate the showers, would lead to different spectra inferred for H and He.
More recently, the KASCADE-Grande experiment [10] has measured the spectra and mass composition in the energy region E ≥ 10 16.6 eV: the spectrum of protons was claimed to show an-ankle like feature at energies E = 10 17.08±0.08 eV, where there is a hardening to a slope of about 2.7. This harder, high energy component was associated with the beginning of the extragalactic CR component. At the same time, the spectrum of Fe nuclei was observed to have a knee-like feature at E = 10 16.92±0.04 eV, which also happens to be at ∼ 26 times the energy of the knee. This feature was readily interpreted as the sign of a rigidity dependent knee in the spectra of the different chemicals. Notice however that the high energy part of the spectra (above the knee) does not appear to fall as an exponential, but rather as a steep power law.
Qualitatively, this behavior seems similar to the one discussed above and connected with the transition from ejecta dominated to Sedov expansion of a SN shell. In Fig. 5 we show the spectra of different elements as calculated using the model discussed above, and compared with the KASCADE-Grande data on protons and Fe. An extragalactic proton component with power law shape E −2.7 has been added to fit the highest energy data points of KASCADE-Grande. For the case k = 9, most appropriate for type II SNe (left panel in Fig. 5), a maximum energy EM ∼ = 3.7 PeV is required. Reaching this high energy requires a SN energy ESN = 2 × 10 51 erg, a factor of a few larger than the standard ESN = 10 51 erg, an efficiency ξCR = 20% and a rate of explosion of ℜ = 1/110 yr −1 . The value of EM is constrained by the need to fit the ankle-like feature on the proton spectrum as measured by KASCADE-Grande. For k = 7, less justified for this type of SNe, the high energy part of the injected spectra are somewhat harder and match the observations better. This is illustrated in the right panel of Fig. 5, where EM ∼ = 781 TeV and ℜ = 1/25 yr −1 .
Recent measurements of the spectrum of the light nuclei and the all-particle spectrum carried out with the ARGO experiment have raised serious doubts about the very foundations of the idea of light nuclei being accelerated to ∼ P eV energies [26,28]. ARGO finds a knee in the light component (H+He) at ∼ 650 TeV, while the   all-particle spectrum is shown to have a knee at an energy compatible with all other measurements. This result appears to be compatible with some previous results obtained by experiments built at altitudes close to the shower maximum. The discrepancy between the ARGO and KASCADE results is an intrinsically observational issue and may be suggesting a rather poor knowledge of the development of showers and question our way to infer mass composition from shower-related observables.
Since it is hard to imagine a way to make ARGO and KASCADE data compatible with each other, we decided to take ARGO data at face value and tried to infer the consequences that would follow. This is illustrated in Fig. 6, where we plot the results of our calculations and compare them with the ARGO data, shown as a shadowed area. The relatively large area derives from the fact that a few different analysis techniques have been adopted so as to establish the independence on the Monte Carlo procedure adopted of the knee-like feature in the light component.
The best fit parameters we find in order to reproduce the ARGO data with our model are ESN = 10 51 erg, ℜ = 1/15 yr −1 , EM ∼ = 507 TeV and ξCR ∼ = 5.2%. While it is clear that the spectrum of light nuclei (H+He) is easily accounted for, it is equally clear that the all-particle spectrum knee cannot be reproduced within a simple approach such as the one discussed here.    [26,28], with k = 9, ESN = 10 51 erg, ℜ = 1/15 yr −1 , EM ∼ = 507 TeV and ξCR ∼ = 5.2%.

Summary
Most of the models for the origin of Galactic CRs still rely on SNRs as being the most likely sources. Yet, the issue of whether SNRs are able to accelerate CR protons up to the energy of the knee, a few PeV, remains open. From the point of view of energetics, it seems that SNRs, in one flavour or another, are realistic candidates. From the point of view of the chemical composition, the anomalous ratio of 22 N e/ 20 N e observed in CRs has led to suggest that the bulk of Galactic CRs could be accelerated in OB associations [33,34,15], although this argument has been recently revised by [44] who reached the conclusion that OB associations cannot contribute the majority of CRs.
In the present paper we focused on the issue of reaching the energy of the knee in particle acceleration during the expansion of a SN shock in the dense region occupied by the wind of a RSG, and stressing the physical aspects of the acceleration process, namely the nature of the instabilities that are required to bring the magnetic field strength and topology to the level that is suitable for particle acceleration. In particular, we investigated in detail the implications of the so called Non-Resonant Hybrid (NRH) instability described by [12,52,13] by computing the maximum energy and the overall particle spectrum produced during the whole SN expansion, both in the ED and ST phases. The general idea is that at any given time, the particles that escape the system with the highest available energy E produce the turbulence needed to scatter the next generation of particles with the same energy, thereby leading to acceleration to higher energies. We discussed the role of the NRH instability in both type Ia and type II SNe, in the context of three models of particle acceleration at the SN shock: 1) a flat spectrum ∝ E −2 ; 2) a hard spectrum ∝ E −(2+β) with β < 0, and 3) a steep spectrum ∝ E −(2+β) with β > 0. The maximum energy and the spectrum of escaping particles was calculated in each case.
The effective maximum energy, which is reached at the beginning of the ST phase in all cases, is found to be at most of order a few hundred TeV for type Ia SNe exploding in the ISM (with the typical values of the relevant parameters). For type II SNe, the situation is quite more complex because these explosions often occur in the wind environment created by the pre-supernova red-giant star. Particle acceleration during the expansion of the SN shock in the dense slow wind of the RSG progenitor is found particularly promising from the point of view of exciting the NRH instability. In this case, the instantaneous maximum energy is shown to always decrease as a function of time after the explosion. However, the ST phase typically starts a few tens of years after the SN event and most of the mass is processed by the end of the ejecta dominated phase. These facts have important consequences. The early beginning of the ST phase allows for a maximum energy which is still high; the effective maximum energy does not coincide with a cut-off, but rather with a steepening of the spectrum, since higher energy particles are indeed produced during the ED phase. Our calculations show that a type II SN with standard energetics (ESN = 10 51 erg, ξCR = 10%) can accelerate CRs up to knee energies, EM ∼ 1 PeV.
It is interesting to appreciate how particle acceleration at PeV energies typically occurs ∼ 10 − 30 years after the SN explosion, which represents a change of paradigm with respect to the standard SNR paradigm where the highest energies are reached several hundred years after explosion. In this scenario, the probability of catching a PeVatron in action in our own Galaxy (for instance through its gamma ray emission) is very low.
For given parameters of the SN explosion, the maximum energy reached at the beginning of the ST phase is a function of the CR acceleration efficiency while the spectrum observed at the Earth also depends on the rate of SN explosions. We calculated the required efficiency and rates necessary to reach the knee and fit the overall CR flux: for steep spectra at the shock, the escape current is lower and this forces one to have larger SN energetics and larger CR acceleration efficiencies, which is counterintuitive for a steep spectrum of accelerated particles. For flat E −2 spectra at the shock, the spectrum of escaping particles is also flat and PeV energies can be reached for ordinary parameters of the SN explosion. Reaching PeV energies is even easier for hard spectra at the shock, however in both the flat and the hard case, one has to face the well known severe problem with CR anisotropy [46,17,49].
The comparison of the expected spectra of light nuclei with the recent data of KASCADE-Grande and ARGO appears very problematic: while individually both sets of data can be fitted with reasonable values of the SN parameters and CR acceleration efficiencies, there is no model that can fit both at the same time. KASCADE-Grande data require a SN energy ESN = 2 × 10 51 erg, an efficiency ξCR = 20% and a rate of explosion of ℜ = 1/110 yr −1 that lead to a maximum energy EM ∼ = 3.7 PeV, whereas we can fit ARGO data with ESN = 10 51 erg, ℜ = 1/15 yr −1 and ξCR ∼ = 5.2% that lead to EM ∼ = 507 TeV. The disagreement between these two experiments is a purely experimental problem, that requires a serious and careful assessment of the systematic uncertainties involved in the adopted experimental techniques.