A&A 409, 799-807 (2003)
DOI: 10.1051/0004-6361:20030909

The knee in the Galactic cosmic ray spectrum and variety in Supernovae

L. G. Sveshnikova

Skobeltsyn Institute of Nuclear Physics of Moscow State University, 119992, Leninskie Gory, MSU, Moscow, Russsia

Received 17 March 2003 / Accepted 16 May 2003

Abstract
This paper proposes a qualitative and quantitative solution of a long-standing problem in astrophysics: the origin of the knee in the Galactic cosmic ray (GCR) spectrum. We calculate GCR flux averaged over Supernova explosion energies and types, applying only the formulae of the standard model of CR acceleration in Supernova remnants (SNR) and the latest astronomical data on the variety in Supernovae. For this purpose we estimate the distribution of SNe in explosion energies and show this distribution to be probably a very asymmetric function with large dispersion. In the case under consideration the cosmic ray flux in the whole energy range should be predominantly formed by the most energetic SN explosions. The knee in the GCR spectrum at energy around $E_{\rm knee}=3$ PeV can quantitatively be explained by the dominant contribution of Hypernovae. The model sketches the all-particle cosmic ray spectrum up to 1018 eV.

Key words: ISM: cosmic rays - stars: supernovae: general - acceleration of particles

1 Introduction

Supernovae represent catastrophic explosions that mark the end of the life of some stars. It is well known that the mechanical energy input to the Galaxy from each supernova is about 1051 erg, so with a rate of about 0.01-0.03 year-1 the total power of SNe in our Galaxy is enough to provide the total energy of Galactic cosmic rays (GCR), $\sim$10-12 erg/cm3(Berezinsky et al. 1990). It is shown that there exists a mechanism needed for the channeling of about 10% (or even more Berezhko & Volk 2000a) of the mechanical energy of the explosion into relativistic particles. Considerable collective efforts have been made during recent years to clarify the mechanism of CR acceleration in SRs (Drury et al. 2001). Theoretical progress is connected to the development of a kinetic nonlinear theory of diffusive shock acceleration (Berezhko & Volk 1997, 2000a; Berezhko 2000b; Berezhko & Ellison 1999; Ellison et al. 1997, 2000; Drury et al. 2001; Malkov & Drury 2001) and the main advances have been made due to improved understanding of the nonlinear reaction effects on the shock structure. This theory can explain not only the main characteristics of the observed all-particle cosmic ray spectrum up to an energy of $10^{14}{-}4\times 10^{14}$ eV (Drury et al. 2001), but also heavy element abundances in CR flux relative to the solar system (Ellison et al. 1997).

The standard model predicts (Drury et al. 2001):

1) The power-like and approximately similar spectra of various nuclei of CRs beyond an energy of 100 GeV/n with the slope $\gamma_{\rm sour}\sim 2.0{-}2.1.$ It may be that spectra have a "curvature'', being even harder before the maximal energy of acceleration, if nonlinear reaction effects are strong (Berezhko $\&$ Volk 2000a; Ellison et al. 1997). Spectra of heavy component can be slightly harder than proton ones due to more effective acceleration of dust grains and ions (Ellison et al. 1997).

2) The maximal energy of accelerated particles  $E_{\rm max}$ (cut-off energy) is $\sim$ $ Z \times 10^{14}$ eV for the average SNe exploding into the average interstellar medium (ISM).

3) There is a possibility to move the maximal energy to higher energies assuming an unusual medium for any class of explosions: explosions into the wind of Wolf-Rayet stars or explosions into superbubbles (Bykov & Toptygin 1997), (see also the reviews by Ptuskin 2000 and Biermann 2000). This effect is mainly due to a higher magnetic field in the stellar or superbubble interior.

4) The real source spectrum inferred from observations after propagation corrections is $\gamma_{\rm sour}=\gamma_{\rm obs}-\Delta \gamma $. The value of  $\Delta \gamma$ varies from 0.3 for a model with reacceleration to 0.8 for a model with Galactic wind (Jones et al. 2001). The value of  $\Delta \gamma$ is not well known yet.

The most "nasty'' problem, as it is called by Drury et al. (2001), is the knee problem, i.e. the origin of a pronounced change in spectrum slope from $\gamma_{\rm obs}\sim2.7$ to $\gamma_{\rm obs}\sim 3.1$ at the energy $E_{\rm knee}\sim 3{-}4$ PeV, discovered many years ago (Kulikov & Khristiansen 1959).

The standard picture makes a clear prediction that the GCR spectrum should start to cut off at energy about 1014 eV or less for all species and drop exponentially as one goes to higher energies (Drury et al. 2001). Only some subclasses of SNR can provide the knee particles while most SNRs have spectra cutting off at considerably lower energies (Reynolds $\&$ Keohane 1999).

The upper limit of acceleration is determined essentially by the product of the shock radius $R_{\rm sk}$, shock velocity $V_{\rm sk}$ (usually normalized to 1000 km s-1), ejected mass $M_{\rm ej}$, remnant age $T_{\rm snr}$, explosion energy $E_{\rm snr}$ (usually normalized to a value of 1051 erg denoted by E51). All these values are connected to each other and vary from explosion to explosion. The cut-off energy per particle  $E_{\rm max}$ can be expressed by the simple formula if only the Sedov phase of SNR expansion is considered (Ellison et al. 1997):

                         $\displaystyle E_{\rm max}$ = $\displaystyle 200 \cdot Z \left(\frac{0.3\cdot B}{3~{\rm\mu G}}\right) \left(\frac{n_{\rm H}}{{\rm cm^3}} \right)^{-1/3}$  
    $\displaystyle \times \left(\frac{E_{\rm snr}}{10^{51}~{\rm erg}}\right)^{1/3} \left(\frac{V_{\rm sk}}{10^3~ {\rm km~s^{-1}}}\right)^{1/3}~{\rm TeV}$  
  = $\displaystyle Z\cdot E_{\rm max}^0(B,n_{\rm H}) (E_{51}\cdot V_{\rm sk})^{1/3}~{\rm TeV}.$ (1)

So the cut-off energy depends on three factors: the first factor Z is the charge of the nucleus, the second factor  $E_{\rm max}^0$ (we introduced it formally) strongly depends on the interstellar medium state, where the SN remnant is expanding: density of protons $n_{\rm H}$ and the magnetic field B; the third factor weakly depends on the energy of the explosion and on the velocity of the shock. The value of B is close to 3 $\mu$G. The "warm'' or "hot'' phases of the medium are selected usually as a common place of explosions (Berezinsky et al. 1990). In the first case $n_{\rm H}\sim 0.3$ cm-3 and $E_{\rm max}^0\sim 100$ TeV. In the second case $n_{\rm H}\sim 0.003$ cm-3 and the cut-off energy is about 3 times that: $E_{\rm max}^0 \sim 300$ TeV. So the value of $E_{\rm max}^0$ can be considered as the cut-off energy for average SN explosion with E51=1 and $V_{\rm sk}=1000$ km s-1. An accurate determination of  $E_{\rm max}$ in a more complete kinetic theory (Berezhko & Volk 1997, 2000a) taking into account timescale evolution of the shock in both the Sedov phase and the free expansion phase of SNR evolution gives, within a factor of 2, the same results (Drury et al. 2001).

A usual way to raise the cut-off energy is to increase the magnetic field in the interior of the progenitor, since sensitivity of  $E_{\rm max}^0$to B is very high. But now it is clear that the parameters of SN explosions can also be varied. Detailed observations of a growing number of supernovae show the nature of this phenomena to be complex (Turatto et al. 2002). Many new peculiar events discovered in recent years display a wide range of luminosities, expansion velocities and chemical abundances, that is evidence for large variations in explosion energy and in the properties of their progenitors (Hamuy 2003).

The main idea of this work is an attempt to obtain the cosmic ray particle spectrum averaged over SNe types and explosion energies. In this case the total CR flux $F(E)={\rm d}N/{\rm d}E$ can be expressed by the formula

\begin{displaymath}{\rm d}N/{\rm d}E=\sum_{i=1}^{\rm Nz} \sum_{j=1}^{\rm Ntp} \... ...^{\rm max}} \Psi_j (E_{51}) G(E, E_{\rm max}) {\rm d}E_{51}, \end{displaymath} (2)

where $\sum_i$ is the summation of different cosmic ray nucleus groups Nz, $\sum_j$ is the summation of different types of SNe explosions Ntp, $\Psi_j (E_{51})$ is the SN distribution in explosion energy inside each SN group within the limits of $E_{51}^{\rm min}\div E_{51}^{\rm max}$. In the main variant we used $E_{51}^{\rm min}=0.1$, $E_{51}^{\rm max}=80$.

$G(E,E_{\rm max}(E_{51}, B, Z))$ is the spectrum of comic rays in every explosion approximated by the power law:

\begin{displaymath}G(E,E_{\rm max})=I_{0} E^{-\gamma}, \end{displaymath} (3)


                                     $\displaystyle \gamma$ = $\displaystyle 2.0~ {\rm in~ the~ interval} ~ 10~{\rm GeV} < E < E_{\rm max}/5;$  
$\displaystyle \gamma$ = $\displaystyle 1.70~ {\rm in~ the~ interval}~ E_{\rm max}/5 < E < E_{\rm max};$  
$\displaystyle \gamma$ = $\displaystyle 5~{\rm in~ the~ interval}~ E > E_{\rm max}.$  

This spectrum shape takes into account nonlinear reaction of CRs to shock structure (Berezhko & Volk 1997; Ellison et al. 1997): decreasing $\gamma$ before  $E_{\rm max}$. $E_{\rm max}$ depends on  $E_{51},~ V_{\rm sk},~ B,~ n_{\rm H}$ according to formula (1). For each type of explosion the B and $n_{\rm H}$ can be different.

Intensity of CRs produced in the each SNR (I0) is found from the condition that the fraction of SNR kinetic energy transformed to CRs is fixed:

\begin{displaymath}\int G(E, E_{\rm max})E {\rm d}E = 0.1\cdot E_{51}. \end{displaymath} (4)

As it was shown in Berezhko $\&$ Volk (2000a), cosmic rays can carry away as much as 30% of the kinetic energy of the explosion.

The function $\Psi (E_{51})$ is the most uncertain distribution; the third section will be devoted to the problem of how to estimate this function. In the second section the latest astronomical observational data on variety in SNe explosions will be reviewed. In the fourth section we present the numerical results of calculations. In the fifth part of the paper we discuss the physical interpretation of the results obtained.

It is worth to note in advance that the diversity of SNe by explosion energies being taken into account results in a very important conclusion: the knee region occurs around several PeV, although only the standard model of CR acceleration and the latest astronomical data on supernovae explosions are used in the calculations.

2 Variety in Supernovae

Due to the growing number of SNe observations, the widely accepted conventional classification of SNe by two types (SNII with hydrogen in the observed spectra and SNI without hydrogen in the observed spectra) has been significantly complicated (Turatto 2003; Hamuy 2003).

The thermonuclear explosions of accreting white dwarfs as they approach the Chandrasecar mass ($\sim$ $ 1.4~ M_\odot$) produce type Ia SNe. Due to their high luminosity and accurate calibration, they are successfully used to determine the geometry of the Universe (Leibundgut 2000), as "standard candles''. The energies of the explosions are practically fixed.

Core collapse supernovae (CCSNe: SNIb/c, SNII) are thought to be the gravitational collapse of massive stars ( $M > 8 ~M_\odot$), which makes the neutron star compact remnants. CCSNe prove to comprise the most common general class of exploding stars in the Universe and they come in a great variety of flavors (Hamuy 2003). Even subclasses of "normal'' SNII: plateau SNII-P and linear SNII-L demonstrate a wide range of explosion energy, from 0.6 to $5.5 (\times 10^{51})$ erg among classical SNII (Hamuy 2003). The ejected masses also are in a broad range between 14 and 56 $M_ \odot$ (Hamuy 2003). Despite the great diversity displayed by SNII-P, these objects show a tight luminosity-velocity correlation. This suggests that while the explosion energy increases so to do the kinetic energies (Hamuy 2003). These stars explode as isolated stars.

A distinct class of SNIIdw can be identified which are believed to be strongly interacting with a "dense wind'' produced by SN progenitors prior to explosion. When the narrow line is present, the SN is classified as IIn ("narrow''). A strong degree of individuality is seen in their spectra, but despite the great photometric diversity among SNIIdw, these objects share the property of being generally more luminous than the classical SNII (Hamuy 2003). Among this type of SNe, one event, SN 1997cy, is much more energetic than any other SNII ( $E\sim 30 \times 10^{51}$ erg, Turatto 2000). SN 1997 cy and its twin SN 1999E (Rigon 2003) are associated with GRBs. As in the case of others, these events show strong ejecta-CSM interactions with explosion energies as high as $3\times 10^{52}$ ergs (Turatto 2003).

Hydrogen-deficient supernovae SNIb and SNIc are associated with the gravitational collapse of massive stars (maybe Wolf-Rayet stars), which have lost their hydrogen envelope during the phase of strong wind. In the case of SNIc most of the helium is gone as well. There is as yet no direct observational proof for binary companions in SNIb/c, but this seems likely (Turatto et al. 2002).

In the past few years 3 SNe (1997ef, SN 1998bw, SN 2002 ap) have been found to display very particular spectra: they are extremely smooth and featureless, which can be interpreted as the result of unusual expansion velocities (Hamuy 2003). This suggests that these objects are hyper-energetic so they are called "hypernovae''. The estimated energies of explosions are very high: $7\times 10^{51}$ erg for 2002 ap (Mazzaly 2002), $8\times 10^{51}$ erg for 1997ef (Nomoto et al. 2000), $60\times 10^{51}$ erg for 1998 bw (Nomoto et al. 2000). The estimated expansion velocity of this object is as high as >30 000 km s-1 (Turatto et al. 2002). SNe 1998bw was not only remarkable for its great expansion velocity and luminosity, but also because it exploded at nearly the same location and time as GRB 980425 (Galama et al. 1998).

Hamuy, in a review (Hamuy 2003), made the following conclusions. Despite the great diversity of core-collapse SNe, several regularities emerge which suggest that 1) there is a continuum in the properties of these objects, 2) the mass of the envelope is one of the driving parameters of the explosion, 3) the physics of the core and explosion mechanism of all core-collapse SNe are not fundamentally different, regardless of the external appearance of the supernova.

The great observational diversity of CCSNe has not been fully understood even if it clearly involves the progenitor masses and configurations at the time of explosion. Whereas SNII-P are thought to originate from isolated massive stars, a generalized scenario has been proposed in which common envelope evolution in massive binary systems with varying mass ratios and separations of the components can lead to various degrees of stripping of the envelope (Nomoto et al. 1995). According to this scenario, the sequence of types IIL Ib Ic is ordered according to a decreasing mass of the envelope (Turatto et al. 2002).

3 How to estimate $\vec{\Psi}$ (E51)

The distribution of SNe by explosion energy $\Psi (E_{51})$ is not yet known. But for the calculation of CR flux, it is enough to get an approximate estimation of this function. For this purpose one can use the absolute-magnitude SN distribution N(Mb) and then transform it to  $\Psi (E_{51})$. We use the data of Richardson et al. (2002), where a comparative study of the absolute-magnitude distributions N(Mb) of supernovae has been done. The authors used the Asiago Supernova Catalog (ASC), where the number of events had increased to 1910 by June, 2001, but the number of events suitable for this study is 10 times smaller. For the absolute-magnitude distribution (in blue filter Mb), the authors consider only SNe within 1 Gpc. These distributions for different types of SNe are presented in Fig. 1.

  \begin{figure} \par\includegraphics[width=8cm,clip]{figure1.eps} \end{figure} Figure 1: Absolute-magnitude (in blue filter) distributions of various types of SNe from Richardson et al. (2002).
Open with DEXTER

The analysis shows that (Richardson et al. 2002):

1) At least 7 of 31 SNe in our Galaxy and in galaxies within 10 Mpc appear to have been sub-luminous ( $Mb \ge -15$). Assuming that there is an observational bias, it appears that more (perhaps much more) than 0.2 of all SNe are sub-luminous, but this fraction remains very uncertain.

2) Only 20 of 297 extragalactic maximum-light SNe appear to be over-luminous ($Mb \le-20$), but it has become clear that they do exist. The absolute-magnitude dispersion of SNIb/c has increased in comparison with previous works due to the discovery of some rather luminous events. The SNe IIn are, on average, the most luminous type of core-collapse SNe. Considering the strong observational bias in favor of them, it is safe to conclude that the fraction of all SNe that are over-luminous must be lower than 0.01.

The authors have approximated absolute-magnitude distributions for each type of SNe by Gaussian functions. They consider also "intrinsic'' distributions obtained taking into account not only Galactic extinction, but also calculating extinction distributed for each SN type, averaged over all galaxy inclinations. Moreover, they divided Ib/c into two luminosity groups: "bright'' and "normal'' ones. The II-L group was also divided into two groups. So 5 or 7 groups of SNe can be analyzed. The parameters of the Gaussian distribution <Mb> and $\sigma (Mb)$ are listed in Table 1 for 5 (7) groups of SNe together with the fractional weight of each group. In the present calculation the fraction of SNIa was decreased to 28% in comparison with the fraction of 60% analyzed in Richardson (2002).

Table 1: Parameters of Gaussian distributions for 5 (7) main types of SNe from Richardson (2002), weight W and value of $E_{\rm max}^0$ used in calculations.

As was pointed out in Hamuy (2003), physics of the core and explosion mechanism of all core collapse SNe are not fundamentally different, so one can expect correlations between average absolute magnitude Mb for given SNe group and its average energy of explosion.

For the estimation of the average dependence E51(Mb) we use the calculation of Nadyozhin (2003), performed in the framework of the LN85 model (Litvinova & Nadyozhin 1983). He makes predictions for correlations between three observable parameters of type II plateau supernovae light curves (the plateau duration $\Delta t$, the absolute magnitude MV measured in V-filter and the photospheric velocity  $V_{\rm ph}$ at the middle of the plateau) and three physical parameters (the explosion energy E51, the mass envelope expelled  $M_{\rm ej}$ and presupernova radius R):

\begin{displaymath}\lg E_{51}=0.135 M_V+2.34 \lg\Delta t+3.13 \lg V_{\rm ph}-4.205\end{displaymath}


\begin{displaymath}\lg M_{\rm ej}=0.234 M_V+2.91 \lg\Delta t+1.96 \lg V_{\rm ph}-1.829 \end{displaymath}


\begin{displaymath}\lg R=-0.572 M_V-1.07 \lg\Delta t-2.74 \lg V_{\rm ph}-3.350, \end{displaymath} (5)

where $M_{\rm ej}$ and R are in solar units and $V_{\rm ph}$ in 1000 km s-1. The analysis of 14 real SNe II-P events shows (Nadyozhin 2003) that the expelled mass, explosion energy and presupernova radius remain approximately within the limits $M_{\rm ej} \sim 10\div30 ~M_\odot$, $R\sim 200\div 600 ~R_\odot$, $E_{51} \sim 0.6\div 2.7$.

We rewrite these formulae to exclude the parameter $\Delta t$ and obtain the following simple relations

\begin{displaymath}\lg E_{51}=-0.43 M_V-0.77 \lg R+0.52\lg M_{\rm ej} -5.83 \end{displaymath} (6)


\begin{displaymath}\lg V_{\rm ph}=+0.57 \lg E_{51}-0.06 \lg R-0.48\lg M_{\rm ej}+1.32, \end{displaymath} (7)

which for the appropriate values $M_{\rm ej}=20$ and R=250 can be written as

\begin{displaymath}\lg E_{51}=-0.43 M_V-7 \end{displaymath} (8)


\begin{displaymath}\lg V_{\rm ph}=+0.57 \lg E_{51}+0.57. \end{displaymath} (9)

Expression (8) will be used in further calculations as a zero approximation for the transformation of the N(Mb) distribution to the $\Psi (E_{51})$ distribution. The sensitivity of the results to this key dependence will be discussed in Sect. 5.

In formula (1) the maximal energy of accelerated CR depends weakly on the parameters E51 and  $V_{\rm sk}$ as $E_{\rm max} \sim Z\cdot E_{\rm max}^0 (E_{51}\cdot V_{\rm sk})^{1/3}$ TeV. To reduce this dependence to the dependence $E_{\rm max} (E_{51})$ we take into account that the highest energy CRs are produced at the end of the free expansion phase (Berezhko 2000b), when the ejecta velocity $V_{0}\sim V_{\rm ph}$. Then we replace  $V_{\rm sk}$ by $V_{\rm ph}/(1-1/r)$ with the compression ratio $r\sim 4-7$ (Ellison et al. 1997) (as in a symple model with a mooving piston) and obtain dependence (10) using $V_{\rm ph}(E_{51})$ (9):

\begin{displaymath}E_{\rm max} \sim Z\cdot E_{\rm max}^0\cdot E_{51}^{0.52}~\rm TeV. \end{displaymath} (10)

Here we neglect the factor $\sim$1.6 to have $E_{\rm max}=Z E_{\rm max}^0$ at E51=1. Expression (10) will be used below in the calculations of average CR flux.

Since the distribution N(Mb) from Richardson et al. (2002) (see Fig. 1) can be represented as a sum of Gaussian distributions with average parameters listed in Table 1, and $\lg E_{51}$ depends linearly on Mb (8), the distribution ${\rm d}N/{\rm d}\lg E_{51}$ can also be represented as a sum of Gaussian functions with parameters

\begin{displaymath}<\lg E_{51}>_j~=-0.43<Mb>_j-7\end{displaymath}


\begin{displaymath}\sigma_j(\lg E_{51})=0.43 \sigma _j(Mb), \end{displaymath} (11)

where j is the type of SN. Then ${\rm d}N/{\rm d}\lg E_{51}$ can be transformed into the distribution ${\rm d}N/{\rm d}E_{51}=\Psi(E_{51})$by multiplying by the factor 0.434/E51.
                       $\displaystyle \Psi(E_{51})$ = $\displaystyle \sum_{j=1}^7 W_j \frac{0.434}{E_{51}\sqrt(2\pi\sigma_j^2)}$  
    $\displaystyle \times\exp\left (-\frac{\left(\lg E_{51}-<\lg E_{51}>_j\right)^2}{2 \sigma_j ^2}\right )\cdot$ (12)

The dependence (8) is used only for core collapse SNe, while for thermonuclear explosions SNIa the values $<\lg E_{51}>~=-0.1$, $\sigma(\lg E_{51})=0.2$ are chosen.

The final distribution $\Psi (E_{51})$ for all types of SNe and contributions of various types are presented in Fig. 2, where the high energy tail is shown separately for a better presentation of the details. It can be seen that most of the events have energies from $(0.5\div 3) \times 10^{51}$ erg. A remarkable peculiarity of the distribution is a very long, tiny tail expanding toward high energies up to E51=60, generally provided by Ib/c1 (bright) and IIn types of SNe (see Table 1). The total fraction of events with energy E51>20 is several percent. The reality and the value of this tail will be discussed in detail in the Sect. 5.

  \begin{figure} \par\includegraphics[width=8.7cm,clip]{figure2.eps} \end{figure} Figure 2: The SN distribution in explosion energies $\Psi (E_{51})$ converted from the Mb-distribution presented in Fig. 1. The high energy part of the distribution is presented separately.
Open with DEXTER

The most abundant sites for SN explosions will be examined below and the values of $E_{\rm max}^0$ for calculation by formula (10) will be selected.

SN explosions are not random in the Galaxy, all of them showing strong spatial concentration toward the center of galaxies and toward the arms in spiral galaxies (Bergh 1997). Numerous regions of very hot and rarefied gas with a temperature $T=5\times 10^6$ K and proton density $n_{\rm H}=0.003$ cm-3 occupy 50% of the volume of spiral arms of Galaxy (Kononovich & Moroz 2001). A source of heating is thought to be the activity of young stars, first of all supernovae explosions. So for young massive stars this site of SNe we selected as the most probable for SNIb/c and SNII. This variant of the interstellar medium ("hot phase'') with parameters T=106 K, $n_{\rm H}=0.003$ cm-3, magnetic field $B=3~{\rm\mu G}$ gives $E_{\rm max} \sim 300$ TeV (Berezhko 2000b).

Numerous regions of neutral HI gas can be divided into two parts (Kononovich $\&$ Moroz 2001): the clouds of gas and dust with $n_{\rm H}=10$ cm-3and T=80 K, occupying a relatively small volume 1%, and intercloud regions, that occupy 50% of the volume of spiral arms with $n_{\rm H}=0.1$ cm-3 and T=104 K. The variant with T=104 K, $n_{\rm H}=0.3$ cm-3, $B=5~\mu$G ("warm phase'') gives a maximal energy of acceleration of about $E_{\rm max}^0= 100$ TeV (Berezhko 2000b). This site we choose as the most probable for SNIa.

Besides that (Drury et al 2001), there exist the temporal correlations resulting from concentration of the majority of core collapse SN progenitors into OB associations. An explosion of the first SN in such an association is followed by several tens of others. This results in formation of multiple supernova remnants powered by both SN explosions and the strong winds of Wolf-Rayet stars in the OB associations, which grows as a large bubble of hot, tenuous plasma known as a superbubble SBs (Tomisaka 1998; Korpi 1999).

The SB acceleration model has been developed by Bykov $\&$ Fleishman (1992), Bykov $\&$ Uvarov (1999). Bykov $\&$ Toptygin (1997) estimated the maximal energy of accelerated nuclei as 1018 eV due to reacceleration effects, in the presence of a magnetic field in the bubble interior of the order of 30 $\mu$G. In this model the spectrum beyond the knee is dominated by heavy nuclei.

Since SNIIn explode in the circumstellar medium (in accordance with the definition) we choose for them the much higher value of $E_{\rm max}^0=4500$ TeV.

The values of $E_{\rm max}^0$ used in calculations in (10) are listed in the last column of Table 1 (this formal selection can be considered only as an example of a possible correlation between the type and the site of SN explosion).

In Sect. 4 we present numerical results of the calculations of all-particle cosmic ray spectrum, using the mentioned above dependences and parameters needed for formula (2). It should be noted that propagation effects were disregarded. The presented source spectra might be easily converted to observable spectra in accordance with the standard model $\gamma_{\rm sour}=\gamma_{\rm obs}-\Delta \gamma $. The value of $\Delta \gamma$ equals 0.3-0.8 depending on the propagation model (Jones et al. 2001).

4 All-particle cosmic ray spectrum, numerical results

All nuclei of cosmic rays are divided in 5 rough groups of p, He, (C, N, O), (Mg, Si, Ne), Fe with relative intensities 0.36, 0.25, 0.15, 0.13, 0.15, respectively. This chemical composition takes into account the fact that heavy components have slightly harder spectra beyond 1 TeV than light ones (Shibata 1995; Hoerandel 2003), the contribution of heavy nuclei increases toward higher energies relative to "normal composition'', obtained around 1 TeV. The spectrum shapes are selected for simplicity in form (3) for all nuclei components.

Figure 3 presents the total proton spectrum generated by 7 different types of SNe with parameters from Table 1, calculated by formulae (2) with $\Psi (E_{51})$ (12) with parameters (11), $E_{\rm max} (E_{51})$ dependence (10). The contribution of each SN type is presented in Fig. 3 separately. The CR intensity dN/dE is multiplied by E2 and presented in relative units.

  \begin{figure} \par\includegraphics[width=8.7cm,clip]{figure3.eps} \end{figure} Figure 3: Total (SUM) proton spectrum (in relative units) generated by 7 different types of SNe with parameters from Table 1. The contribution of each SNe type is presented separately.
Open with DEXTER

It can be seen from Fig. 3 that:

1) The contributions of most energetic explosions are stressed significantly due to expression (4) in our calculations: the total number and total energy of accelerated cosmic rays is proportional to the total kinetic energy of the explosion. For example, only 5% of cosmic rays are generated by SNIa, while they comprise about 30% of all SNe.

2) The total intensity of CR is practically formed by contributions of SNIb/c1 (bright branch, see Table 1) and SNIIn.

3) The location of the knee is determined by maximal energy of accelerated CR protons in the most energetic explosions.

All obtained features can be understood if analytical expression for the average value of  $E_{\rm max}$ is written using formula (10). The statistical weight of  $E_{\rm max}$ should be proportional to the total number of accelerated CRs, i.e. $\sim$E51 (as can be seen from (4)).

                        $\displaystyle <E_{\rm max}>$ = $\displaystyle \frac{\int\limits_{E_{51}^{\rm min}}^{E_{51}^{\rm max}} \Psi (E_{... ...{E_{51}^{\rm min}}^{E_{51}^{\rm max}} \Psi (E_{51}) \cdot E_{51} {\rm d}E_{51}}$  
  = $\displaystyle {E_{\rm max}^0} \frac{<E_{51}^{1.52}>}{<E_{51}^1>}~ {\rm TeV}.$ (13)

For the symmetric distribution function $\Psi_j (E_{51})$ with a small dispersion, as for the case of SNIa, the factor <E511.52>/<E511> is close to unity. For a very asymmetric function, as in Fig. 2, this factor can be many times larger. For example, if we choose a power-like function $\Psi(E_{51})=0.44\cdot E_{51}^{-1.7}$ in the interval $E_{51}=1\div 100$ and $\Psi(E_{51})=0.44$ in the interval $E_{51}=0.1\div 1$, the value of <E511.52>/<E511> is $\sim$5. (This power-like shape of $\Psi (E_{51})$ will be discussed in Sect. 5.)

As can be seen from Fig. 3, the first knee in the proton spectrum is located around 3 PeV, while for the most probable energy of explosion E51=1 maximal energy $E_{\rm max} \sim 300$ TeV (10).

The second step in the proton spectrum is formed (as it is seen Fig. 3) by the contribution of SNIIn explosions, because they are also very energetic (see Fig. 2) and they have a much larger value of $E_{\rm max}^0$ (see Table 1).

Figure 4 presents the spectra of different cosmic ray nuclei (the Mg-Si-Ne group was omitted from Fig. 4). In Fig. 5, the contributions of various types of SN explosions to the all-particle spectrum are shown.

  \begin{figure} \par\includegraphics[width=8.7cm,clip]{figure4.eps} \end{figure} Figure 4: Spectra of different cosmic ray nuclei.
Open with DEXTER

Every nuclear component also has two steps shifted to higher energy by factor Z in comparison with protons. The all-particle spectrum beyond the knee is formed by the sum of coupled steps.

The change of the spectrum slope $\delta \gamma$ beyond the knee (in the interval 3 PeV to 26 $\cdot$ 3 PeV) is determined by a fraction of Fe nuclei w(Fe) in the chemical composition of CRs before the knee in the case when the contribution of SNIIn is small enough:

\begin{displaymath}\delta \gamma= \frac {\lg (1/w({\rm Fe}))}{\lg 26}\cdot \end{displaymath} (14)

For $w{\rm (Fe)}\sim 0.15{-}0.20$ the change of the spectrum slope is close to $\delta \gamma \sim 0.5$.

If the contribution of SNIIn is large, $\delta \gamma$should be less. Besides, as it can be seen from Fig. 5, the more the diversity in explosion energies (as in the case of SNIIP), the smoother the behaviour of the spectrum beyond the knee.

  \begin{figure} \par\includegraphics[width=8.7cm,clip]{figure5.eps} \end{figure} Figure 5: The contributions of various types of SN explosions into all-particle spectrum.
Open with DEXTER

The maximal energy of accelerated Galactic CRs is determined by the Fe nuclei generated in SNIIn explosions and located around 1018 eV. The chemical composition of CRs in the region beyond the knee and up to 1018 eV should be heavier than one in the region before the knee.

In Fig. 6 the mean logarithmic mass $<\ln~A>$ of CR (usually used to characterize the mass composition) is presented compared to the data obtained in the KASCADE experiment (Kampert et al. 2002). The main variant (when the chemical composition of CRs generated in SNIIn is the same as for others) predicts fewer heavy component in the range 1016-1017 eV than in the KASCADE experiment. But in accordance with Bykov & Toptygin (1997), the CRs originated in superbubbles can be enriched by heavy nuclei. Figure 6 presents a variant, when proton and helium components are absent in CRs originating in superbubbles. The experimental dependence of  $<\ln A>(E)$ lies between these two variants.

  \begin{figure} \par\includegraphics[width=8.7cm,clip]{figure6.eps} \end{figure} Figure 6: The average mass number $<\ln~A>$ of CRs in the KASCADE experiment (full squares) and in the calculations: full line - main variant, dashed line - variant when CRs from SNIIn are enriched by heavy nuclei.
Open with DEXTER

5 The physical interpretation of results

The most significant problem in the present model is the reality of the long tail in $\Psi (E_{51})$ and the sensitivity of the knee location to this tail. To analyze the second question, we present in Fig. 7 the proton spectrum of CRs with different upper limits integrating $E_{51}^{\rm max}$ in formula (2). In the main variant $E_{51}^{\rm max}=80$ is used. Figure 7 shows that the point of the knee location is determined by the SNe explosions with $E_{51}\sim 30\div 60$. At $E_{51}^{\rm max}=10$ the knee is around 700 TeV, that is higher than 300 TeV due to Eq. (13), but it is not enough to reproduce experimental data.

  \begin{figure} \par\includegraphics[width=8.7cm,clip]{figure7.eps} \end{figure} Figure 7: The dependence of the knee location on maximal energy of SN explosions.
Open with DEXTER

The discovery of SNe with enormous kinetic energy (Hypernovae) is one of the most interesting recent developments in the study of SNe (Nomoto et al. 2002). They can be directly identified by the explosions determining the cosmic ray knee region not only by the energy of explosions, but also by the type of core collapse SNe. In our calculations, only SNIb/c1 (bright branch from Table 1) and SNIIn give the principal contribution to the formation of the knee region (see Figs. 3, 5). Among 7 possible Hypernovae, 5 have been recognized as type Ic (1998bw, 1997ef,1997dq, 1999as, 2002p) and 2 as type IIn (1997cy, 1999E) (Nomoto et al. 2002). The Hypernova branch might be interpreted as follows. Stars with $M>20{-}25~ M_\odot$ form a black hole as a compact remnant; whether they become hypernovae or faint SNe may depend on the angular momentum in the collapsing core, which, in turn, depends on the stellar winds, metallicity, magnetic field and binarity. Hypernovae might have rapidly rotating cores possibly due to the spiraling-in of a companion star in a binary system (Nomoto et al. 2002).

To test the obtained function $\Psi (E_{51})$, we use the sample of 26 real SNe from Hamuy (2003), where the physical parameters (explosion energy, ejected mass, radius of progenitors) for real supernovae of types II, Ib/c, IIdw are presented. The integral distribution W( $>\!\!E_{51})$ of real SNe was constructed and analyzed. All 26 SNe have an energy of more than 1, so it had to be asumed that the fraction of events with E51>1 is 0.6 among all SNe, that is close to the value in our calculation. (It worth noting that the calculation is sensitive mainly to the slope of the tail in the $\Psi (E_{51})$ function, but not on the absolute value of this normalization factor.) This distribution is denoted as the "real SN'', while in reality it depends on the basic theoretical premises and it can be distorted by the selection bias. In Fig. 8 it is presented together with $W(>\!\!E_{51})$ used in the main variant of our calculation with $\Psi (E_{51})$ shown in Fig. 2.

The "real SN'' distribution can be fitted by power-like function with a slope of -0.78. This means that the differential distribution is $\Psi(E_{51})\sim E_{51}^{-1.78}$ at E51>1. In the region of E51=0.1-1, $N(E_{51})=\rm const.$ was chosen. An all-particle spectrum with "power-like'' functions  $\Psi (E_{51})$ identical for all core collapse SN groups is presented in Fig. 9 also.

Since the physical picture of the Hypernovae explosion can differ from other core collapse SNe, we considered a variant, when in E51(Mb) (6) the parameters of ejected mass and radius of progenitors were chosen differently for SNIb/c1, SNIIn (see Table 1) and other SNe: $M_{\rm ej}=4 ~M_\odot$, $R=80 ~R_\odot$ for SNIb/c1 and SNIIn, but $M_{\rm ej}=10 ~M_\odot$, $R=600 ~R_\odot$ for other types. The weights of SNIb/c1 and SNIIn groups were decreased by a factor of two in comparison with Table 1. This variant is denoted as "Hypernovae'' in Figs. 8 and 9. It fits very well the form of "real SN'' distribution.

  \begin{figure} \par\includegraphics[width=8.7cm,clip]{figure8.eps} \end{figure} Figure 8: The integral distributions $W(>\!\!E_{51})$ for different variants of calculations in comparison with the "real SNe'' distribution from Hamuy (2003).
Open with DEXTER

Figure 9 shows the all-particle CR spectra calculated for all three variants of $W(>\!\!E_{51})$. Examination of Figs. 8 and 9 shows that the input of events with energy $E_{51}\sim 30{-}50$ is $\sim$$ 2\pm 1$% in all cases. One can see that the difference in fraction of events with $E_{51}\sim 10$ among all SNe in the "main'' and other variants of calculations is not crucial for the all-particle CR spectra presented in Fig. 9, because the knee is mostly sensitive to the energy of explosion $E_{51}\sim 30{-}50$, as is shown in Fig. 7. But the sharpness of the knee depends noticeably on the shape of $\Psi (E_{51})$. For the case of a power-like $\Psi$ distribution for all types of SNe the knee is smooth, but for the case when hypernovae are singled out by energy ("hypernova'' variant) the knee looks more sharp.

The experimental all-particle spectrum from Hoerandel et al. (2003) is also presented in Fig. 9 in comparison with calculations. This experimental spectrum was multiplied by E0.65 to take into account the propagation corrections and to reduce the observable spectrum to source spectrum. Here we increase the value of $E_{\rm max}^0$ for SNIIn to $1.5\times 10^{16}$ eV (in comparison with Table 1) to get the better coincidence with the experimental spectrum in the interval 1017-1018 eV.

The present calculations reproduce on the whole the all-particle spectrum measured in EAS experiments: the knee around 3 PeV, the change of slope by $\delta \gamma \sim 0.3{-}0.5$ beyond the knee, start of dip around 1018 eV and may be the knee shape. In the considered model it is possible to obtain the sharp form of the knee, if one selects the narrow distribution in explosion energy for Hypernovae. As it has been pointed out in Erlykin $\&$ Wolfendale (1997), the sharpness of the knee, measured in some individual EAS experiments is quite noticeable.

In general the contribution of high energy SNe to the total CR flux, needed for the explanation of the knee region, can be estimated from basic considerations. If the total power of 100% SNe with average energy E51=1in our Galaxy is enough to provide the total energy of Galactic cosmic rays $\sim$10-12 erg/cm3 (Berezinsky et al. 1990), then 2% of SNe with E51=50 or 3% of SNe with E51=30 can also provide the total energy of Galactic cosmic rays.

  \begin{figure} \par\includegraphics[width=8.7cm,clip]{figure9.eps} \end{figure} Figure 9: All-particle spectra calculated with different variants of the $W(>\!\!E_{51})$ functions presented in Fig. 8 (lines) in comparison with the experimental all-particle spectrum (stars) averaged over many EAS experiments by Hoerandel (2003).
Open with DEXTER

We can draw a conclusion that the fraction of events responsible for formation of the knee in the CR spectrum comes to $\sim$$ 2\pm 1$%. It means that with a usual SN-rate of about 10-2 year-1, the Hypernova rate should be about $2\times 10^{-4}$ year-1. Then, taking into account that the lifetime of CR life in our Galaxy is about $3\times 10^7$ year, one can get that about $6\times 10^3$ explosions provide the intensity of low energy CRs in our Galaxy. But for high energy CRs (around 1 PeV), the lifetime of CR is much shorter due to decreasing of escape length  $\lambda_{\rm esc}$ as $E^{-\alpha}$ ( $\alpha=0.54$) beyond the energy 5 GeV (Jones et al. 2001). The number of explosions giving the dominating input to CRs around $E \sim 1$ PeV might be rather small ($\sim$10-15) in the whole Galaxy and maybe only few explosions provide the CR flux in the Solar system.

The latter conclusion may coincide with the idea proposed by Erlykin $\&$ Wolfendale (1997) that a single nearby local SNR accelerates the particles and gives the dominating input (mainly by O and Fe nuclei) to the knee region. But CRs from this SNR reach the Earth directly without distortion by propagation effects. In our model the most energetic explosions give the dominant contribution to the whole energy range of CRs and their propagation in the Galaxy should be taken into account, while these effects are not considered in the present paper.

6 Conclusions

1.
We calculated the Galactic cosmic ray spectrum averaged over Supernova explosion energies and types, based only on the formulae of the standard model of CR acceleration in Supernova remnants and on the latest astronomical data on the variety in Supernovae. For this purpose we estimated the distribution of SNe in explosion energies and show that this distribution is probably a very asymmetric function with large dispersion. In the case under consideration the cosmic ray flux in the whole energy range should be predominantly formed by the most energetic SNe explosions.

2.
The knee in GCR spectrum at energy around $E_{\rm knee}=3$ PeV can quantitatively be explained by the dominant contribution of SNe with an energy of ($\sim$ $ 30{-}50)\times 10^{51}$ erg, that might be identified with Hypernovae. The estimated rate of these energetic explosions is about 1-3 per 104 years, enough to provide the total power of Galactic cosmic rays.

3.
In the proposed model the location of knee is determined by an abrupt fall of protons generated in the most energetic SNe; the change of a spectrum slope beyond the knee (in the interval $E_{\rm knee}\div 26\times E_{\rm knee}$) is mainly determined by the subsequent abrupt fall of other nuclei; the contribution of CRs generated by SNe exploding into a specific circumstellar medium with a large magnetic field becomes predominant in the range $E> 26\times E_{\rm knee}$. As a result, each nuclear component of CRs should have two steps, and the all-particle spectrum beyond the knee is formed by a sum of double steps. The chemical composition of CRs in this region is much heavier than in the region before the knee.

Acknowledgements
I acknowledge numerous discussions on this problem in private conversations and at seminars with many physicists: K. A. Postnov, V. S. Ptuskin, A. D. Erlykin, G. T. Zatsepin, O. G. Ryazhkaya, M. I. Panasyuk, N. N. Kalmykov, T. M. Roganova, N. V. Sokolskaya, V. Prosin, L. A. Kuzmichev, A. K. Managadze, I am very grateful to all members of the RUNJOB and NUCLEON collaborations, where I am engaged, for discussions and support. This work is supported by the RFBR grant N 03-02-16272.

References


Copyright ESO 2003