The Origin of the Most Energetic Galactic Cosmic Rays: Supernova Explosions into Massive Star Plasma Winds

: We propose that the high energy Cosmic Ray particles up to the upturn commonly called the ankle , from around the spectral turn-down commonly called the knee , mostly come from Blue Supergiant star explosions. At the upturn, i.e., the ankle , Cosmic Rays probably switch to another source class, most likely extragalactic sources. To show this we recently compiled a set of Radio Supernova data where we compute the magnetic ﬁeld, shock speed and shock radius. This list included both Blue and Red Supergiant star explosions; both data show the same magnetic ﬁeld strength for these two classes of stars despite very different wind densities and velocities. Using particle acceleration theory at shocks, those numbers can be transformed into characteristic ankle and knee energies. Without adjusting any free parameters both of these observed energies are directly indicated by the supernova data. In the next step in the argument, we use the Supernova Remnant data of the starburst galaxy M82. We apply this analysis to Blue Supergiant star explosions: The shock will race to their outer edge with a magnetic ﬁeld that is observed to follow over several orders of magnitude B ( r ) × r ∼ const ., with in fact the same magnetic ﬁeld strength for such stellar explosions in our Galaxy, and other galaxies including M82. The speed is observed to be ∼ 0.1 c out to about 10 16 cm radius in the plasma wind. The Supernova shock can run through the entire magnetic plasma wind region at full speed all the way out to the wind-shell, which is of order parsec scale in M82. We compare and identify the Cosmic Ray spectrum in other galaxies, in the starburst galaxy M82 and in our Galaxy with each other; we suggest how Blue Supergiant star explosions can provide the Cosmic Ray particles across the knee and up to the ankle energy range. The data from the ISS-CREAM (Cosmic Ray Energetics and Mass Experiment at the International Space Station) mission will test this cosmic ray concept which is reasonably well grounded in two independent radio supernova data sets. The next step in developing our understanding will be to obtain future more accurate Cosmic Ray data near to the knee , and to use unstable isotopes of Cosmic Ray nuclei at high energy to probe the “piston” driving the explosion. We plan to incorporate these data with the physics of the budding black hole which is probably forming in each of these stars.


Introduction
Energetic particles far above thermal levels, called Cosmic Rays (CRs), were discovered by Hess (1912) [1] and Kohlhörster (1913) [2]; their energies are now known to range up to about 10 20 eV, and at that energy scale first seen by Linsley (1963) [3]. It was proposed very early (Baade and Zwicky 1934 [4]) that Supernovae (SNe) could easily provide the energies required for Galactic Cosmic Rays (GCRs). Ginzburg and Syrovatskij (1963) [5] suggested that radio galaxies could provide the highest energy component. The mechanism of acceleration was originally proposed by Fermi (1949Fermi ( , 1954 [6,7] as repeated reflections between interstellar clouds with magnetic field bottle configurations which are approaching each other. This was later generalized to today's concept of energetic particles being reflected back and forth across a shock wave, and separating converging flows (Axford et al., 1977;Krymskii 1977;Bell 1978a,b;Blandford and Ostriker 1978; [8][9][10][11][12]), with early general overviews by Ginzburg and Ptuskin (1976) [13] and Drury (1983) [14]. It is now recognized that such particles permeate interstellar space, as well as much of intergalactic space. Particles with energies above order GeV are found throughout the Solar system. The early arguments on the possible origin of CRs were all made on the grounds of sufficient energy supply, and they still stand. However we now know of many other processes and events around stars and black holes that clearly produce energetic particles. Energetic particles are also produced on the Sun, and in the Solar wind, teaching us basic plasma physical concepts.
In addition to the supernova origin of accelerated CRs and outflow (the theme of this paper), is energy outflow and particle acceleration due to central supermassive black holes (SMBHs), and their jets and lobes. They produce CRs, magnetic fields and energy. They have been discussed by, e.g., Ginzburg and Syrovatskij (1963) [5], Lovelace (1976) [15], Biermann and Strittmatter (1987) [16], Kronberg et al. (2001) [17], Kronberg (2002) [18], Kronberg et al., 2004 [19], Colgate (2004) [20], Colgate et al. (2014) [21], Biermann et al. (2016) [22], and in other papers. These phenomena clearly can produce much higher particle energies, and might explain the highest energies detected above the ankle. We note that, since this paper confirms that massive star explosions as well as SN Ia explosions are just two sources of magnetic fields in galaxies and, via galactic winds, in the cosmos. Bigger outflows, such as from super-massive black hole activity via, e.g., radio galaxies, are an additional source.
In Section 2 we develop the theory of massive star plasma winds with the goal of deriving characteristic ankle and knee energies. In Section 3 we compare the energetics of the resolved radio super nova remnants (SNRs) in M82 and demonstrate consistency with our predictions. In Section 4 we discuss a wide variety of concerns involving generation and transport, both theoretical and observational, in both a Galactic and extragalactic context, including alternate acceleration mechanisms found in the literature. We stress that our goal in this paper is to demonstrate one path to understand the CR spectrum, without excluding other mechanisms.

Massive Star Plasma Winds
Massive stars explode as Supernovae; and they have magnetic plasma winds. The shock pushed out by the Supernova (SN) explosion can accelerate CR particles. When this shock races through a stellar magnetic plasma wind, accelerated CR particle energies can go as high a few 10 18 eV , referred to as CR-IV, [39]). Massive stars with magnetized plasma winds come in two varieties, Red Supergiant (RSG) stars above about 25 M Zero Age Main Sequence (ZAMS) mass, and Blue Supergiant (BSG) stars above about 33 M ZAMS mass. All of these stars are in binary systems, but not all of them are in tight orbits. RSG star magnetic plasma winds have a low velocity and high density, whereas BSG star magnetic plasma winds have a very high velocity and low density. During the time preceding the SN explosion, this plasma wind experiences short episodes during which the star may enhance its mass loss and its magnetic fields (as summarized with many references, e.g., in ASR18 [34]). A fair number of both RSG star explosions as well as BSG explosions have now been observed, at many wavelengths and spatial resolutions. Some data are collected in ASR18 [34]. The observation that the strength of the magnetic field in the SN shock is the same for RSG explosions as for BSG explosions is not consistent with all mechanisms that enhance the magnetic field just in the SN shock (Bell and Lucek 2001) [40]. In such a case the magnetic field should be an order of magnitude stronger in RSG star explosions, since the implied magnetic field energy density would scale with the ram pressure in the SN-driven shock. (For a detailed discussion including the possibility of a small contribution from the Bell-Lucek mechanism, see ASR18, [34]). We note that the review by Helder et al. (2012) [41] on SNRs and CRs was written before most of these detailed data became available. A newer corresponding article on CRs in stellar-wind bubbles by Zirakashvili and Ptuskin (2018) [42] does not use any of these data either. The arguments in Jokipii (1982Jokipii ( , 1986) and Zirakashvili and Ptuskin (2018, [42]) all refer quite explicitly to weak particle scattering, whereas Jokipii (1987, [45]) refers to the strong scattering limit, used here as well as in earlier work as cited. The plasma wind terminates in a shell, where it encounters the interstellar medium (ISM) and the molecular cloud (MC) out of which the star formed, and/or the wind-bubble environment of preceding SN explosions in its immediate environment. As the ram pressure of RSG star plasma winds is two orders of magnitude lower than for BSG star plasma winds, the wind shell may have a radial scale about an order of magnitude smaller than BSG star plasma winds. It follows that the life time of activity for RSG star explosions is an order of magnitude shorter than for BSG star explosions. This provides a strong counter selection effect in identifying such sources. We note that BSG star explosions are very rare in our Galaxy, and yet often enough to be able to contribute significantly to the observed CRs. When the SN shock does hit the wind shell, it slows down considerably; for radio supernovae (RSNe) the SN-shock speed is typically 0.1 c, while for the Supernova Remnants in M82 (Kronberg et al., 1981 [46], 1985, hereafter referred to as KBS85 [47]) the currently observed velocities of expansion (e.g., Fenech et al., 2008Fenech et al., , 2010) are all far below this value; this observation matches the expectation that when encountering the wind shell the SN plasma shock reverts to environment-dominated expansion (e.g., Cox 1972 [50]). That reversion causes the shocked region to be about 1/4 of the spherical volume encompassed by the shock, rather than 3/4, as in the shocked region in a wind (Biermann and Cassinelli 1993 [51], Biermann and Strom 1993 [52], referred to as CR-II and CR-III; this uses a strong shock, see Landau and Lifshitz 1959 [53], and Drury 1983 [14]). When a SN shock races through a stellar plasma wind, it produces radio emission at resolutions reachable with interferometry (i.e., sub-arcsecond radio imaging). This permits better specification of the magnetic field B(r), the stellar mass lossṀ , its shock radius r, and velocity V SN /c. We note that V SN is the velocity of the shock in the medium through which it propagates, here the wind of the star. As these shocks are observed to show about 0.1 c, we can neglect the prior velocity of the wind. From these three quantities, scale r, magnetic field B, and shock speed V SN /c we can derive two critical energies for the maximal energies particles can reach in shock acceleration, E ankle and E knee . We first show the maximal energy E ankle of a particle fitting into the space available. For a plasma simulation see Meli and Biermann 2006, [54], and for an analytical derivation see, e.g., Biermann 1993, referred to as CR-I, [55], or ASR18 [34].
where e is the elementary charge, Z is the charge of the particle under consideration (cgs units).
We note here that a Parker configuration of a magnetic field in a wind implies that any SN-shock racing through the wind encounters a perpendicular magnetic field (Parker 1958 [56]). We emphasize that this interpretation requires and predicts that all participating stars have about the same magnetic field strength B 0 in their SN-shock racing through the wind (see CR-I [55]). B 0 is the magnetic field strength at a nominal radius of r 0 , in observations of radio supernovae, usually taken to be 10 16 cm. This can be understood as the Larmor motion fitting into the space available. We have used here this plasma wind property which is (Parker (1958) [56]) This is consistent with the observations discussed here. Here B ϕ is the dominant component of the magnetic field (Parker 1958, Weber andDavis 1967 [56,57]). We interpret the associated energy as that turnup in the spectrum of Galactic CRs (GCRs) which is called the ankle. This corresponds to energies reached in perpendicular shocks (see Jokipii 1987, [58], CR-I, CR-IV, Meli and Biermann 2006, [39,54,55], and ASR18 [34]). The argument is simply that the Jokipii limit suggests κ ∼ r g V SN for the scattering coefficient, where r g is the Larmor radius, and V SN is again the shock speed. Using shock acceleration theory requires the time scale of acceleration ∼ {κ/V 2 SN } ∼ {r/V SN }, which derives from the fact, that in a perpendicular shock particles can only escape via strong turbulence, driven by the shock speed. This then cancels the factor containing the shock speed, and gives the spatial or Hillas limit (Hillas 1984 [59]) limit. In parallel shock acceleration the same argument gives then an additional factor of (V SN /c) 2 , since in the shock region the escape time scale is then r/c. This is because CRs can stream out along the field lines in the highly non-stationary and un-relaxed post-shock region, and the scattering coefficient does not carry the shock speed anymore (for another derivation of this expression see CR-I [55], using the Parker magnetic field configuration, see Parker 1958 [56]). It is obvious that this ankle energy is independent of radius r as long as the Parker configuration is maintained. That is one step in the argument which this paper proposes to check. The lower break energy, given by a parallel shock configuration, corresponds then to the knee, where the spectrum turns down: This energy is then also constant with radius r, as long as the SN-shock maintains its velocity V SN . For the observed numbers of RSNe (ASR18 [34]) the energies derived match the energies for ankle and knee directly to within the errors, without using any parameter adjustment. This expression follows for plane-parallel shock configurations (Drury 1983 [14], Völk and Biermann 1988 [60]). This energy is the ankle energy multiplied by (V SN /c) 2 . Some small fraction of the surface of the SN-shock can be expected to have a parallel configuration of the magnetic field in the form of temporary islands due to turbulence. The Parker solution also has a polar cap in its magnetic field configuration, where the radial magnetic field becomes locally dominant B r > B ϕ (see, e.g., Biermann et al., 2009 [61]). We note that to obtain these numbers we take the observed data and average the products r × B(r) and r × B(r) × (V SN /c) 2 . The numbers show that the ratio of these two averages is 10 −1.6 between them, instead of 10 −2 , as would be expected by multiplying separate averages, suggesting that these numbers might be weakly correlated.
However, as these observations show the situation typically around r 10 16 cm, due to opacity effects, we need to ascertain what happens at the larger scale, when the SN-shock actually hits the wind shell or other environment. We need to know whether the compact radio sources in M82 can be understood as the late stage development of the same kind of supernova explosions as those observed directly in other galaxies. This is possible to find out with the data from M82: We propose to test this concept with the observations of the SNRs in M82 (KBS85 [47]; Bartel [70]).

The Supernova Remnants in the Starburst Galaxy M82
Here we show the data on all supernova remnants (SNRs) in M82 that have published derived estimated magnetic field strengths, and for which the size was clearly established (Allen andKronberg 1998, Allen 1999 [65,66] with many further references therein): The flux density values listed in Table 1 for the non-thermal emission are referenced to a frequency of 1 GHz, and were obtained by a fit across very many frequencies, from 408 MHz all the way to 23. 46 GHz, using a screen model absorption and separating thermal from non-thermal emission. The spectral index given is that of the non-thermal straight power-law component. We note that this magnetic field strength determination posits that the combined energy density of CR electrons and heavier CR particles and the magnetic field strength is a minimum (see, e.g., Miley 1980 [71]). This corresponds roughly to equipartition between energetic particles and the magnetic field, and can be understood on the basis of stability arguments (such as in Parker 1966 [72] and Ames 1973 [73]). The magnetic field thus derived, in Allen and Kronberg (1998) [65], has to be scaled with the factor{(1 + k)/Φ} 2/7 , which is difficult to measure; here k is the ratio of the CR proton energy density to that of CR electrons, and Φ is the volume filling factor. At lower particle energy than a few GeV this ratio is expected to be smaller, since protons do not contribute any extra energy below their rest mass energy, but electrons continue to contribute due to their smaller rest mass. From above we adopt the uncertain filling factor Φ = 1/4, since clearly the SN-shock has encountered the outer shell of the wind, is observed to have slowed down, and so is not filling the entire spherical volume inside the radial scale given. We list here the magnetic field strength without any such correction. Table 1. Supernova remnants (SNRs) in starburst galaxy M82, based on Allen and Kronberg (their Table 5; 1998) [65], in turn based on Allen (1999) [66].

The CR Proton/Electron Ratio k
To determine the magnetic field strength based on radio data it is important to have an estimate of the proton/electron factor k. For the CR proton/electron ratio k, observed to be about 100 above a few GeV, we can make the following argument, based on the observations of young RSNe: The first process that happens in an electron-proton shock is the thermalization of these particles (see also Bell 1978b [10] and Drury 1983 [14] for analogous descriptions, developed before the new radio spectral data were known, and before the importance of drifts was recognized, see Jokipii 1987, [58]). Thermalization has to work at all radii, from below 10 16 cm out to parsec scale, rendering the density lower by a factor of 10 5 at least. The wind is magnetic with the basic configuration of the magnetic field a Parker type spiral (Parker 1958 [56]), so essentially yielding a perpendicular shock configuration. This process has been discussed in Spitzer (1962) [74], but his treatment is non-relativistic and uses isotropic phase space distributions; then simple electron-proton scattering is far too slow at high temperature. When the protons go through the shock they fill, in phase space, a torus with the radius in momentum m p V SN , and with the relative width of the inverse of the Mach-number, so corresponding to the upstream speed of sound. That entails that this anisotropic momentum distribution causes turbulence at all the scales from the maximum on down, over a range of the Mach-number. However, because the Mach-number is close to m p /m e the smallest scales in this range correspond to the waves matching the scales of the electrons with the maximal momentum m e V SN . This means that at all these scales the protons and electrons are scattered in momentum phase space. We have argued in ASR18 [34] that the spectrum of irregularities in the direct shock region is I(k) k ∼ 1 and so the power per log bin of wave-number is equal at all scales. It follows that thermalization occurs. Due to all this turbulence the time scale of thermalization is of order a small multiple of the thermal Larmor radius divided by the Alfvén speed, so highly adequate, by a factor of 2 E knee /(m p V 2 SN ), and a fortiori 2 E ankle /(m p V 2 SN ). The time scale for electrons might be m p /m e times slower due to the anisotropy of the scattering, but that is still quite adequate by many powers of ten. The time scale scales with the radius, and works similarly at all relevant radii, the condition noted above.
Thermalization yields a first characteristic electron energy γ e,sh c 2 m e . Beyond this energy, our interpretation of the radio observations then says that the electrons gain energy only by drifts (CR-I [55]); shock acceleration does not work yet, because the electrons do not see the shock. In ASR18 [34] we noted that the electrons only get Stochastic Shock Drift Acceleration (StSDA), giving them a spectrum of −3.06, from observations, say about −3 (this line was reasoning was also used for ions by Caprioli 2012 [75]). This drift-dominated spectrum extends to the momentum at which the energetic electrons reach a Larmor radius that equals that of the shock thermal protons, γ sh,inj c m e , (if another CR ion is dominant, corresponding to the Larmor radius of that ion, with mass number A dom and charge at injection Z dom ). From this second momentum electrons "see" the shock, and have the normal CR accelerated spectrum, say −7/3 for the most common CR component of CRs accelerated in a wind (see papers CR-I, CR-II, CR-III, CR-IV, and ASR18 [34,39,51,52,55], and the earlier references cited there). Protons attain the maximal energy content only at a few GeV (i.e., just above their rest mass energy). This applies to any particle number spectrum slightly steeper than −2. That implies that electrons have a higher energy content at that momentum where they still are just above thermal in the post-shocked region. We assume that at this stage electrons are in equipartition with the protons (or dominant ions) at the start, due to charge neutrality. We posit that in the initial post-shock thermal phase electric neutrality holds.
This means that and So between these two energies the energy density of CR electrons runs down by a factor of 2 β sh 1 Z dom (6) in the first step. Here we use the observations of a spectral index of about E −3 implying that the energy content above a given energy runs as E −1 , giving a factor of about 20, with the observed β sh = 0.1. We interpret this as due to Stochastic Shock Drift Acceleration (StSDA, first used in CR-I [55]).
There is an observational check on this, as the electron energies required to explain the radio observations of RSNe need to be above that first energy: We observe the emission at a frequency of 5 GHz at a magnetic field strength of order Gauss, requiring a Lorentz factor of γ e 30, while γ e,sh is about 10, so well below this. However, this argument implies that the emission should have a cutoff on the low frequency side, perhaps visible somewhat early in the emission phase and at a lower frequency, such as 400 MHz. It also follows that at that stage of the emission there are constraints on what ions could possibly dominate to explain the emission, implying that in this numerical example A dom 3, eliminating basically all common elements other than Hydrogen (We used A dom = 1 and Z dom = 1 in our numerical examples).
In the second step, using a spectrum of p −7/3 the electrons reduce their energy in that range where the proton contribution is not significant. They are below relativistic, so up to the energy of the proton rest mass (or dominant ion again) the factor is which is about 2, using the same data again, with protons. So the CR electron energy density is reduced by a combined factor of k 40 when their momentum reaches that of the relativistic protons, where this factor is determined by observations. Already Allen and Kronberg (1998) [65] had mentioned 40 as a possible and indeed plausible value for k.
The argument makes it obvious that other sources, for which β sh is smaller, contribute even less. On the other hand, if the shock speed is larger, i.e. the factor can be smaller and electrons can contribute more, consistent with many other well established arguments (e.g., Bell 1978b [10], Drury 1983 [14]). However, usually a parallel shock configuration is used (where the magnetic field is parallel to the shock normal). Here we emphasize perpendicular shocks and Stochastic Shock Drift Acceleration (StSDA), in a standard magnetic field MHD configuration in a SN-shock racing through a magnetic stellar wind as shown by Parker (1958) [56] (see also, e.g., CR-II, McClements et al., 1997 [76], and Dieckmann et al., 2000 [77]). That in turn implies that the combination of various sources contributing to electrons and protons can reduce the electron net contribution even more, down to the observed factor of 1 in 100. This line of reasoning starts with the energetic electrons having exactly the same energy density as the dominant ions in the shock, due to charge conservation at thermal energies in the thermal post-shock phase. To the degree that protons constitute a strong contribution from ISM-SN-CRs, which have a slightly steeper spectrum than the most common wind-SN-CR component as well as slower shocks, this suggests that at moderate energies the CR electrons have a spectrum slightly flatter than protons.

The Magnetic Field in the SN-Shocked Wind, and Implications
All this suggests that the proton/electron ratio k = 40 is appropriate as a fiducial parameter to derive the magnetic field strength from radio observation via the minimum energy condition. This leads to a correction factor of 4.29 over the magnetic field strength values listed in Allen and Kronberg (1998) [65], and repeated above. However the shock must have slowed down significantly from its free expansion phase through the wind. Using the ISM-Sedov volume fraction of the shocked shell of 1/4, we need to recognize that the magnetic field is likely to have been squashed by a new reverse shock pushed back through the material. This gives the magnetic field strength an extra factor of 4, effectively cancelling the factor of 4.29 mentioned above. The magnetic field strength which we wanted to test was the value before any additional enhancement, over the pure wind-shock mode. Therefore we should reduce the magnetic field strength obtained by a factor of 4, coming back to the original numbers by Allen and Kronberg (1998) [65]. This supports the argument by Allen and Kronberg (1998) [65] to give their magnetic field strength value without such an uncertain correction.
We emphasize, that the values of r × B(r) in this table show a rather narrow distribution (Allen 1999 [66]), strongly suggesting that this product is not only the same for various massive star explosions, but also that this product is maintained in their evolution. The average of the last column is 16.0: Using the calibration from Allen and Kronberg (1998) [65] and their flux density values we can also estimate the corresponding value of < log{B × r} > from the independent table of SNRs in the starburst galaxy M82 in Fenech et al. (2010) [49], who did not include a separation of thermal and non-thermal emission, nor any absorption. However, to within the errors, we obtain the same number with a larger scatter confirming the analysis done by Allen and Kronberg (1998) [65]. The uncertainty in < log{B × r} > due to the different distances to M82 used by Gendre et al. (2013) [70] and Allen and Kronberg (1998) [65] is 0.05. Just using the formula given for the magnetic field estimate in Miley (1980) [71] implies that {B × r} ∼ D 5/7 , where D is the distance used. Using the specific numbers here, a factor of 10 0.1 higher than in ASR18 [34], gives for the ankle energy now using eV for energy as is custom in the CR community. This specific number for the error was derived from the data given in paper ASR18 [34]. Thus we have two independent estimates of the product < log{B × r} >, one in exploded RSG stars as well as BSG stars, at a typical radius of 10 16 cm in various other galaxies, and one at a typical scale here of order pc, 3 × 10 18 cm in the starburst galaxy M82. The fact that these two numbers coincide allows the interpretation that the population of SNRs in the starburst galaxy M82 is just the later stage of the same explosions (see Allen 1999 [66]). We note that this behavior differs significantly both from a standard Sedov expansion into the ISM (see, e.g., CR-III [52]) as well as from an expansion model of a relativistic blob for compact radio source expansion (see, e.g., van der Laan 1966 [78] and Shklovskii 1960 [79]), for which B × r ∼ r +1 or ∼ r −1 , respectively. This is for the radial scale differences here a factor of either 10 +2.5 or 10 −2.5 . So the Parker solution is not only maintained across many powers of ten in radius (Allen 1999 [66]), but these explosions have the same magnetic field strength in the SN-shock region racing through the stellar wind of stellar explosions in our Galaxy, other galaxies, and in M82. It seems as we show below only BSG star winds can reach such a radius of parsec scale in the environment of M82. This view is by far the most straightforward, and simplest interpretation of the data. One might then ask why we do not observe such explosions in our Galaxy. These are very rare in our Galaxy. It is not all surprising that we have not seen one during the observations with modern tools. As noted size selection effects make the more common RSG star explosions very difficult to detect; BSG star explosions are quite rare compared to the time scale of modern observations. However, there is a possibility that the recent PeV source seen by H.E.S.S. (High Energy Stereoscopic System) is just such an explosion (Abramowski et al. (2016) [80]). We note that BSG explosions can contribute significantly during the CR storage time, which is much larger than the historical observation time scale, even at the highest energy.
Observations of the CR spectrum, mass composition and arrival direction anisotropy will be important from below the knee, such as with IceCube, HAWC and ISS-CREAM, between knee and ankle, and to higher energies above the ankle (see, e.g., Aab [87,88]): This will allow us to discern the origins of GCRs and the transition to extragalactic CRs with the Pierre Auger Observatory and the Telescope Array.
Two immediate questions remain: These are (i) whether the plasma wind maintains its magnetic profile B ϕ ∼ r −1 (Parker 1958, Weber and Davis 1967, [56,57]), and (ii) whether the SN-shock slows down (see Kronberg et al., 2000 [67] for variability). Simple estimates of the ram-pressure in RSG star winds and BSG star winds show that a RSG wind shell may be an order of magnitude smaller than a BSG star wind shell. That difference entails that in BSG star winds the SN-shock has the chance to race all the way to the wind shell before converting into a Sedov expansion (Cox 1972 [50]). In RSG star explosions this seems to be a more remote possibility, depending on the outside environment; a Sedov expansion implies conservation of energy (Sedov 1958 [89]). The sparse data on velocities (e.g., Fenech et al., 2008Fenech et al., , 2010) show that all measured velocities of the SNR expansion in M82 are below 0.1 c, usually far below. In the data listed above in Table 1 there is no hint of a correlation between size and spectral index, possibly suggestive of a changing shock speed and Mach number. Just from simple mass function arguments, the ratio of RSG star explosions to BSG star explosions is expected to be of order 1 to 2. However, in RSG star wind shells the density is far higher than in BSG star wind shells. This makes for more efficient cooling. It implies that we have four development phases: (i) First we have the racing phase, when the SN-shock goes through the plasma wind at constant velocity. (ii) For RSG star winds we may have an intermediate wind-Sedov phase, when the shock velocity slows down. (iii) Then we have the phase when the shock hits the outer shell and becomes a ISM-Sedov expansion. (iv) Finally the shock region transforms to a cooling shock. For RSG star explosions these stages are much shorter than for BSG star explosions. Consistent with such an expectation, the typical radius is about 1 pc, with only one source at about 0.1 pc, which may be a RSG star explosion, or a BSG star explosion in a very high pressure environment. During the ISM-Sedov phase the radio luminosity is constant, deviating as soon as cooling sets in. The observations appear consistent with this (Kronberg et al., 2000 [67]). The only exception with a well defined light-curve is 41.9 + 58 (Bartel et al., 1987 [62], Kronberg et al., 2000 [67]). It is the source suspected to be a recent mis-aligned GRB (Muxlow et al., 2005 [69]; see also below), and/or possibly a binary BH merger (ASR18 [34]).
We have also derived the required piston mass needed to keep driving the shock; the number we obtained from γ-ray line observations was about 0.1 M (ASR18 [34]). One question is whether this could be the 56  The computed mass accumulated in a BSG star wind can be written as whereṀ −5 is the mass loss of the star prior to the explosion in units of 10 −5 M yr −1 , and V W,8.3 is its wind velocity in units of 2000 km/s. The mass loss is directly determined from modelling the radio observations of the radio supernovae (the references in Tables 1 and 2 of ASR18 [34]). The wind velocity is typical for such stars. This suggests that a piston of 0.1 M could drive the SN-shock unimpeded to scales of order 10 parsec as long as the ram pressure of the wind is sufficient to overcome the environment. Another consistency test of our concept and its numbers is that the kinetic energy should significantly exceed the magnetic field energy. A piston of 0.1 M at 0.1 c corresponds to 10 51 erg, while the typical magnetic field strength of about 5 mGauss corresponds to 10 49.5 erg << 10 51 erg, using the volume fraction of 1/4 adopted above.
For the ram pressure we need to consider first the wind ram pressure against the outside medium, and then the SN-shock ram pressure. The wind ram pressure condition for matching the ISM pressure is where P ISM,−12 is the outside pressure in units of 10 −12 dyn cm −2 , a typical number for the outside pressure in our Galaxy. The higher pressure in M82 would allow smaller sizes for the wind shell, but of course the local pressure, where the SNRs are seen, may be lower than in the central region. A pressure 10 3 times higher would still allow the size scales observed for the SNRs in M82. The region containing the SNRs has been described as a ring of about 500 pc radius (Weliachew et al. 1984, Kronberg and Wilkinson 1975, and Kronberg et al. 1981, [46,95,96]). If the over-pressure inside that ring were two orders of magnitude higher, for instance as argued in KBS85 [47], then the maximal size would only be about 3 pc, near the upper end of the range of the SNR radii actually seen. The corresponding SN-shock ram pressure condition gives r < 10 2.7 pcṀ 1/2 where V SN,−1 is the shock speed in units of 0.1 c. Another check is when the mass accumulated outside within the extra radial range ∆ r slows down the shock's progression. That is given by ∆ r r = 10 −0.7 M pist,−1 r −3 pc n −1 (12) where M pist,−1 is the mass of the piston in units of 0.1 M , r pc is the radius in pc when the SN-shock hits the outer medium, and n is the outside medium density in units of cm −3 . That implies that once the SN-shock hits the wind shell and outside medium at full speed, the shock slows almost immediately to the Sedov stage, and only shortly later to the cooling phase. Adding in the mass accumulated in the wind does not modify this conclusion significantly. In a BSG star wind nothing earlier slows down the SN-shock. The knee energy predicted here is We conclude that, within the errors these numbers of the product of radial scale and magnetic field strength (Allen 1999 [66], here in Table 1 and in Tables 1 and 2 in ASR18 [34]) are the same. This result supports our argument that in these examples the SN-shock pushed out without losing significant strength, kept following the Parker regime (Parker 1958 [56], Allen 1999 [66]), and so maintained these energies. All or most of these examples in M82 are consistent with BSG star explosions. We note again that we identify massive exploding stars in many other galaxies, the starburst galaxy M82, and in our Galaxy as all being essentially the same, and assume that they produce the same CR spectrum-the data shown here and earlier certainly allow such an approach. All of the SNRs seen in M82 are consistent with being in the slow-down phase, which is after the SN-shock has hit the wind shell or other environment. We note that other physical concepts may also yield similar values for the quantity {B × r}, for instance for SN Ia explosions, where the transition from free expansion to Sedov expansion readily yields the same value, thereafter with a slow diminishing over time in the expansion.

The Origin of High Energy Galactic Cosmic Rays
As discussed above: RSG star explosions produce High Energy (HE) GCRs, i.e., energies near to, but with a kink down below the knee, and with a final cutoff below the ankle, with near normal abundances but not extending to the highest energies. BSG star explosions produce HE GCRs with all heavy elements enhanced up to the highest energies for GCRs (see, e.g., Todero Peixoto et al. 2015 [97], Thoudam et al. 2016 [98]). If our simple picture as deduced above, is correct, that RSG star explosions just terminate earlier at a radius about 1/10 of that of BSG star explosions. The anti-proton data (ASR18 [34]) in our Galaxy suggest that the RSG star explosions run their shock down in velocity, whereas the BSG star explosions hit their respective wind shell at full speed. We note that this particular distinction may be different in the galaxy M82 due to the higher environmental over-pressure (see Allen 1999 [66]). Therefore the characteristic energies, E ankle and E knee are maintained only for the BSG star explosions. It ensues that at the highest CR energies, RSG star explosions contribute far less than the BSG star explosions. Since in our model these energies scale with nuclear charge number Z, at a down-turn of the spectrum, the heavier elements go into dominance near the knee (CR-IV [39]); we note that this concept of 1993 required and predicted that the product {B × r} has a similar value in all relevant explosions of very massive stars exploding into their wind. At the final cutoff, at the ankle energy, protons will disappear very much earlier than Carbon, Oxygen or heavier elements such as Iron (see Figure 6 in CR-IV, as well as Figure 1 in ASR18 [34,39]). We emphasize again that the observed characteristic energies E knee and E ankle are directly reproduced by two independent data sets on RSNe and SNRs, as both are proportional to {B × r}: This seems confirmed by a third independent data set (Fenech et al., 2010 [49]). It is based on both RSNe in other galaxies, and SNRs in the starburst galaxies M82, without assigning any special parameter a particular value. We do assume that the CR population and its characteristics are the same in our Galaxy as in other galaxies, and in the starburst galaxy M82. The two data sets agree well within the errors. The spectrum both below and above the knee can be explained with a combination of Diffusive Shock Acceleration (DSA) and Stochastic Shock Drift Acceleration (StSDA), as shown first in CR-I, CR-II, CR-III, and CR-IV, rederived in ASR18 [34,39,51,52,55], all based on earlier work by Drury (1983), Jokipii (1987), Völk and Biermann (1988), [14,58,60] and others. One consequence is that most interaction happens in the wind-shell, when the SN-driven shock hits (e.g., Biermann et al., 2001 [99]). This concept can be tested with the spectra of secondary particles. Here we note that the concept requires that the γ-ray -spectra of the Galaxy should reflect this early interaction, showing evidence for a spectrum of mostly E −7/3 from the 4 π component in wind-SN-CRs, and a weaker polar cap component of E −2 . The Fermi data suggest two such spectral components (de Boer et al., 2017 [100]). HAWC and H.E.S.S. data are fully consistent with this prediction of a dominant 4 π component (e.g., Nayerhoda et al., 2018 [101]). At higher γ-ray energies the resulting spectrum is predicted to become flatter due to the emergence of the polar-cap component with E −2 in wind-SN-CRs (CR-IV [39]). Other new data can also be used for further tests, e.g., from AMS (Alpha Magnetic Spectrometer on the ISS) or HAWC, such as the spectra of various elements and secondary and primary contributions to Ni in cosmic rays (ASR18, Alfaro et al. 2017, Aguilar et al. 2017, 2018a).
Nonetheless several unanswered questions remain in this simple picture. We list some of them in the following. All the processes discussed up to here may happen as argued. However, they do not have to dominate, as there are many alternatives, especially considering SN Ia, as shown below.
(1) Very rarely, at a rate of order 10 −2 of massive stars, we do have Gamma Ray Bursts (GBRs) in any star forming galaxy, i.e. for order 100 very massive BSG star formed, there is one GRB. There is no question that they accelerate particles to high energy. Where is the contribution from GRBs in our own Galaxy? Could they also contribute below the ankle? Or is this already part of the extragalactic contribution, as argued in many papers, such as Usov (1995, 1996) [106,107]; Vietri (1995Vietri ( , 1996Vietri ( , 1998 [108][109][110]; Waxman (1995), Miralda-Escudé and Waxman (1996), Waxman and Bahcall (1997) [111][112][113]; Zhang and Mészáros (2004) [114]; Piran (2004) [115]? Based on the radio maps showing some evidence of a jet-and BH-spin-flip, we proposed in ASR18 [34] that the compact radio source 41.9 + 58 in the starburst galaxy M82 (KBS85 [47]) is actually the result of a recent merger of two stellar black holes (BHs). This could have happened concurrent with a GRB (see Muxlow et al., 2005 [69], who first argued for a misaligned GRB for this source). The argument then proceeds to propose that this BH merger could explain the Ultra High Energy CR (UHECR) particles detected by the Telescope Array (TA) detector (Abbasi et al., 2014(Abbasi et al., , 2018a [117][118][119]) that GRBs do not contribute significantly to the UHECR population. All this depends on the model assumptions, and so may not be in contradiction.
(2) There are many arguments in support of UHECRs coming from radio galaxies (e.g., Ginzburg and Syrovatskij 1963, Hillas 1984, Biermann and Strittmatter 1987, Rachen and Biermann 1993a, Rachen et al., 1993b, Meli et al., 2008 [125,126,128,129]). One of these candidates, TXS0506+056 is the most convincing, since it shows evidence for several neutrino events. To make HE neutrinos requires the prior acceleration of HE and UHECR particles first, so this tentative finding allows the speculation that mergers of super-massive black holes may lead to ubiquitous production of UHECR particles. This entails the notion that UHECR production is highly episodic. Do these sources also contribute below the ankle?
(3) BSG explosions can yield the energy, but can the CRs produced actually reach us? This is especially an issue for electrons and positrons: we have argued that we need electrons up to about 30 TeV in the source region; H.E.S.S. has now detected electrons at 20 TeV, published on its website [130], and in conference proceedings. The CR electron spectrum shown on the H.E.S.S. website is consistent with our 2009 prediction (Biermann et al. 2009 [61]) of at lower energies first E −10/3 , then E −3 , and last, at the highest energies, E −4 . These are all based on the idea that many sources contribute under several loss mechanisms (Kardashev 1962 [131]). At the highest energies the observed spectrum appears to be slightly flatter than E −4 , possibly due to local and/or recent nearby sources. We note that the synchrotron loss time of electrons at 20 TeV is of order 3 × 10 4 yrs (using a magnetic field strength of order 5 µGauss: Beck et al., 1996 [132]; a crude estimate of 5 µGauss was first made in the 1940s, then using equipartition with CRs as necessary for isotropy), implying a source rather close by (as has been argued by many, e.g., Yüksel et al., 2009Yüksel et al., , 2012), or a very anisotropic transport to allow larger distances. An alternate way of phrasing this speculation could be that the magnetic fields are extremely structured, allowing some transport over large distances.
(4) The relative contribution of CR particles from RSG star and BSG star explosions at low energy is about 4:1, following Binns et al. (e.g., Murphy et al., 2016 [135]). What is it at high energy? The anti-proton data suggest that the RSG star explosions shocks run out of steam, and so that RSG star explosions in our Galaxy contribute much less at the maximal energies than the BSG star explosions. This may imply that the piston mass in RSG star explosions is smaller than in BSG star explosions, consistent with a naive interpretation of the 56 [90][91][92][93][94]) be identified with the piston (see also, e.g., Biermann et al., 1992 [136] and Wesson et al., 2015 [137])? The numbers seem to match rather well, but have not been derived for BSG star explosions, and rely mostly on rather simple 1D models. It is aparent that we do not actually know the explosion mechanism for very massive stars, be it due to neutrinos (e.g., Bethe 1990 [138]), magnetic fields (e.g., Bisnovatyi-Kogan 1970, Bisnovatyi-Kogan et al., 2008 [139,140]), some combination of the two mechanisms (à la Seemann and Biermann 1997 [141]), or-another suggestion-quark deconfinement (Fischer et al., 2018 [142]).
(6) Do SN Ia explosions contribute anything to the observed CR particle population? The D 6 scenario seems currently to be favored (Shen et al., 2018a [143]; D 6 here stands for Dynamically Driven Double-Degenerate Double-Detonation): This argument has been strengthened by identifying some dwarf stars in the Gaia data as remaining from SN Ia explosions (Shen et al., 2018b, Raddi et al., 2018). The radio data (e.g., Dickel et al., 1991 [146]) and the interpretation of the radio and X-ray data (e.g., Chevalier 1982 [147], Archambault et al. 2017 [148]) demonstrate unequivocally that these explosions do accelerate energetic particles (see also Bykov et al., 2018 [149]).
Assuming a fraction of 0.1 ε B,−1 of the kinetic energy of the ejecta in SN Ia E SN,51 10 51 erg (e.g., Dwarkadas and Chevalier 1998;Mazzali et al., 2007;Townsley et al., 2009: [150-152]) goes into the amplification of magnetic fields following Lucek and Bell (2000), Bell and Lucek (2001), Bell (2004Bell ( , 2008 [40,[153][154][155], we obtain, in a Sedov expansion (Cox 1972) [50] Obviously, in this mechanism there is a corresponding energy density of the energetic particles. The mass M ISM reached by the SN shock can be reached by a typical ejected mass of about M ej 1 M (e.g., Dwarkadas and Chevalier 1998;Mazzali et al., 2007;Townsley et al., 2009;[150-152] (15) and can readily be reached in free expansion. At larger radii Sedov expansion implies a slow lowering of the quantity {B × r}. Therefore the quantity {B × r} at first rises linearly, as long as the shock is in free expansion, and so the shock velocity constant, to rise to its maximum: That is remarkably close to what we observe in M82, and is only weakly dependent on any single given parameter. Thereafter this quantity declines as r −1/2 , until the Sedov expansion shock reverts to a cooling shock, at a larger radius (see, e.g., Cox (1972) [50]).
Thus explosions of SN Ia would seem to be a viable alternative to explain the observations of M82. However, SN Ia explosions are not enhanced in their rate immediately in a starburst (Scannapieco and Bildsten 2005;Raskin et al. 2009; [156,157]). Only very massive stars explode nearly simultaneously, as only their evolutionary time scale is shorter than a typical starburst time of just a few million years (e.g., Biermann and Fricke 1977;Huchra 1977;Huchra et al., 1983; Alonso-Herrero et al., 2001; [158][159][160][161]). Very massive stars show the same magnetic field behavior, i.e., the derived values of {B × r} are consistent with the same number over several powers of ten in radius, from near 10 16 cm to parsec scale. How do we then exclude the mechanism proposed by Lucek and Bell (2000), Bell and Lucek (2001), Bell (2004Bell ( , 2008, [40,[153][154][155]? Since the SN shock is argued to run at undiminished velocity over the radial scale from 10 16 cm to parsec scale, the Bell-Lucek mechanism would in fact also produce a constant value of {B × r} across this radial range, matching all these observations; a condition for this would be a SN shock expansion into a Parker limit wind (Parker 1958, [56]), fulfilled for both RSG and BSG star explosions into their winds. However, that is contradicted by the interpretation of the observations, as SN shocks racing through the dense winds of Red Super Giant (RSG) stars showing the same magnetic field as those racing through the tenuous winds of Blue Super Giant (BSG) stars (see Bell 2004Bell , 2008; ASR18; [34,153,154]). A straightforward conclusion is that all these SNRs in M82 are produced by BSG star explosions, which come from the most massive normal stars. In ASR18 [34] some options are discussed for the Bell-Lucek mechanism to contribute to the observed magnetic fields. This would allow a possible understanding that both RSG and BSG stars show a similar magnetic field behavior. However, this line of reasoning assumes that all such SNRs and RSNe derive from the same class of stars, the most massive stars: The quantity {B × r} turns out to be similar (i) at its peak value for SN Ia explosions, using the Bell-Lucek mechanism, and (ii) for the most massive stars generally. Therefore the option is viable to assume that CRs derive in our Galaxy from SN Ia explosions and in M82 and other such galaxies from the most massive stars: However, SN Ia explosions only approximately give the same energy for knee and ankle due to the time dependence of {B × r}. Thus using many such explosions would imply considerable smearing at the knee and ankle energies, as well as a somewhat different number for the knee energy. This may be ruled out already (e.g., Thoudam et al., 2016, [98]), but should be carefully checked out to confirm or refute this option. What is the spread of individual knee energies among the population of contributing explosions? (7) Pulsars and pulsar wind nebulae are directly observed to accelerate particles to high energy (see, e.g., for a hint of neutrinos the TeVPA 2018 IceCube lectures; e.g., Aliu et al., 2014 [162]). What is their contribution? The high energy gamma data data clearly demonstrate that they produce energetic particles, but they cannot contribute easily to the observed positrons (Yüksel et al., 2009;Abeysekara et al., 2017 [133,163] [149], in ASR18 [34], and in earlier papers cited therein. They surely contribute, since they are a variant of GRBs; however, is their contribution mostly at the very highest energies? (9) Active single and binary stars: Any massive star with a powerful wind contributes (Seo et al. 2018 [165]). Micro-quasars are binary stars, as are all massive stars (Chini et al. 2012(Chini et al. , 2013). Micro-quasars surely produce energetic particles and cosmic rays (e.g., Heinz and Sunyaev 2002, Mirabel 2004). In sum the final evolutionary stages will contribute to energetic particles (e.g., Vulic et al., 2018, Abeysekara et al., 2018b [171,172]). Is their contribution discernible in the CR data? (10) Neutron star and stellar mass BH mergers: This may lead to a GRB. As observed they likely contribute to the observed CRs (Abbott et al., 2017a,b; Albert et al., 2017b [173][174][175]). However, if such CRs are produced in other galaxies, what path do they take through the combined magnetic fields of intergalactic space with its filaments, halo, wind and disk of our Galaxy (see, e.g., Pshirkov et al., 2011 [176], Mao et al., 2012 [177], Yüksel et al., 2012 [134], Kronberg 2013 [178])? (11) Transport out of the Galaxy? Radio data suggest that this CR transport is typically convective, i.e., in a galactic plasma wind (Heesen et al., 2018 [179], Krause et al. (2018) [180] and Miskolczi et al. (2018) [181]). They give a threshold condition in the form of star formation rate per area (see also Rossa and Dettmar 2003 [182]). This confirms some very old ideas, going back to Weber and Davis (1967) [57], Parker (1958) [56] and even earlier notions in the late 1940s.
(12) A CR contribution from the wind-shock of the Galactic plasma wind? This has been explored repeatedly (Jokipii and Morfill 1987 [45]), most recently by Merten et al. (2018) [183]. (13) Neutrino production in the Galaxy? The IceCube and Antares observations have begun to constrain the models (Albert et al., 2018 [184]). The γ-ray data already suggest plausible models, all consistent with the concept that interaction is dominant when the SN-shocks hit the wind-shell (Biermann et al., 2001). It suggests that interactions during CR transport through the ISM of the galaxy are less important (Nath et al., 2012 [185]), as had been suspected for a long time. Analogously to our comments on diffuse γ-ray emission there should be two neutrino spectra, one at E −7/3 from the 4 π component, and at higher energy a lower E −2 from the polar-cap component (CR-IV and de Boer et al. 2017, [39,100]). The corresponding neutrino emission spectra have a transition at an energy which is poorly determined right now. At TeV energies we cannot easily separate the RSG and BSG contributions where the latter may or may not dominate the neutrino emission. The stronger interaction in the RSG star explosions (leading to anti-protons in our proposed model) may push the transition energy to a higher value in the observations when summed over a line of sight; we note that this approach may clarify the anti-proton spectrum.
The overall diffuse Galactic neutrino emission is predicted to show a sharp turn-off of both components near 1/20 of the knee energy (for protons), so here IceCube should be able to test this prediction in the future. (14) What is the explosion mechanism, and do we obtain a BH every time, most of the time, or only occasionally when BSG stars explode? What is the nature of the piston? The innermost layers that get ejected may come from a zone just outside the budding black hole. Is this process the source of the strong magnetic fields observed, coming out with the piston? Are there some stars that produce a black hole without any significant explosion (see Mirabel 2017a,b,c [186][187][188])? Do we have a chance to learn about BH physics from these explosions?

Summary
The main task that we have set ourselves was to check whether the characteristic energies implied by the numbers for radius, magnetic field strength, and shock velocity deduced from observed Radio Supernova (RSN) data (listed in ASR18 [34]) can be confirmed with the observations of the compact sources, interpreted as Supernova Remnants (SNRs) in the starburst galaxy M82 (KBS85, Allen andKronberg 1998, Allen 1999, [47,65,66] as well as earlier and later references): The answer seems to be positive. In fact we can confirm not only that the Parker magnetic field solution (B × r = const) is maintained over many powers of ten in radial range of the SN-shock racing through the stellar wind for massive stars, but that they all have the same magnetic field strength, in our Galaxy, other galaxies and in starburst galaxies such as M82. The energies for the downturn, called the knee, and the upturn, called the ankle, seen in the Cosmic Ray (CR) spectrum can be derived directly from the data, again without any free parameter. That supports the interpretation that the observed CRs in the energy range from below the knee all the way to the ankle may derive directly from Blue Supergiant (BSG) star explosions. It also suggests that the general CR spectrum is the same in our Galaxy, other galaxies and in starburst galaxies such as M82. Given that BSG star explosions are suspected to produce black holes, and to drive the piston which gives rise to CRs, can CRs and their detailed properties teach us about BH formation?