\documentclass[12pt]{article}
\usepackage{geometry} % see geometry.pdf on how to lay out the page. There's lots.
\geometry{letterpaper}
%\geometry{a4paper}
\geometry{margin=1in}
% or letter or a5paper or ... etc
% \geometry{landscape} % rotated page geometry
\usepackage{times}
\usepackage{amsmath}
\usepackage{hyperref}
% See the ``Article customise'' template for come common customisations
\title{Diffuse Galactic gamma-ray emission; sources and challenges}
\author{L. O'C. Drury}
\date{Draft work in progress - for comments, not distribution!} % delete this line to display the current date
%%% BEGIN DOCUMENT
\begin{document}
\maketitle
\tableofcontents
\section{Introduction}
The diffuse gamma-ray emission from the Galactic disc is the most powerful tracer we have for the global distribution of cosmic rays in the Galaxy, the only other signal with any comparable diagnostic power being the radio synchrotron emission tracing the low-energy cosmic ray electrons.
In the distant future neutrino telescopes may challenge this, but for the moment this is the best we have.
Some constraints can be placed on the low energy ion component from interstellar chemistry and the heating of ISM clouds, and global energetics and pressure balance considerations provide other weak constraints, but as far as the Galactic distribution of cosmic rays is concerned the key data sets to analyse are the gamma-ray and radio maps. The analysis requires of course key inputs from other astronomical domains; in particular, the distribution of interstellar gas and the interstellar magnetic field. It also requires theoretical inputs in the form of source and propagation models.
\section{Unresolved sources?}
One question, not entirely semantic, is what exactly do we mean by diffuse emission? Is the diffuse emission simply the emission left over when we subtract all the identified sources? This is clear enough if the sources are unresolved point sources, but less clear if the sources are themselves spatially extended (as expected and seen for Galactic SNRs for example).
In other fields of astronomy we have seen many cases where diffuse backgrounds have ultimately been shown with better instruments to contain a significant contribution from unresolved sources and there is no reason why the Galactic gamma-ray emission should be any different. That being said there are also very strong theoretical arguments that there should be a genuinely diffuse emission from the dust and gas of the ISM resulting from its irradiation by the general interstellar cosmic ray flux. There are three principal mechanisms for this emission. The dominant process at high energies (above about 1 GeV) is gamma-ray production from the decay of neutral pions produced by the hadronic interaction of cosmic ray ions with interstellar atoms. Cosmic ray electrons also contribute through Bremstrahlung as a result of collisions with interstellar atoms and via inverse Compton up-scattering of photons and these processes are important in making detailed fits to the emission, especially at lower energies.
\section{Sources}
The standard model for the origin of the Galactic cosmic rays is that they originate in strong interstellar shocks through diffusive shock acceleration. A major, but not the only, source of such shocks is of course supernova explosions and SNRs are widely believed, for good reasons, to be the source for the bulk of the Galactic Cosmic Rays (GCRs). Perhaps the strongest indirect argument for a link between supernovae and the origin of cosmic rays is the fact that estimates of the power needed to maintain the observed Galactic cosmic ray flux are typically a few to ten percent of the mechanical energy input into the Galxy by supernova explosions. This, together with the fact that there is no other obvious energy source capable of driving the cosmic ray acceleration, strongly suggests a link between the two phenomena, a suggestion reinforced by the observation of radio synchrotron emission from supernova remnants which clearly points to the acceleration of electrons to GeV energies. All of this was already well known in the 1960s and is well discussed, for example, in the influential monograph by Ginzburg and Syrovatsky \cite{GS}. However there are other sources of strong shocks, such as young star clusters, OB association winds, in-falling high-velocity clouds, superbubbles etc and it would be strange if they did not also make a contribution. Indeed HESS appears to have detected at least one young stellar cluster as a potential source \cite{Westerlund2}.
In addition to acceleration by strong ISM shocks there may well be additional sources, in particular of electrons (including positrons). Pulsars clearly produce nebulae of synchrotron emitting electrons and some of these may escape into the general ISM. The much-debated PAMELA excess may be an indication of such an additional local source of high-energy electrons. Magnetic reconnection is also a real possibility; although there is no well developed theory, it is clear that in the sun and the Earth's magnetosphere reconnection does occur and is associated with particle (mainly electron) acceleration. It would be remarkable if it did not also occur in the Galaxy and indeed some reconnection is probably essential for a Galactic dynamo to operate.
\subsection{Current status of Diffusive Shock Acceleration theory}
The origins of Diffusive Shock Acceleration (DSA) theory is usually traced to four seminal papers which appeared independently in 1977-1978, \cite{K, Ba, Bb, BO, ALS} although there is some ``pre-history", for example the inspired suggestion by Fred Hoyle as far back as 1960 that interstellar shocks might put their energy not just into gas heating, but also into magnetic fields or cosmic rays \cite{Hoyle} (but of course Hoyle had no specific mechanism in mind and his paper was purely a what-if speculation).
It is important to understand that DSA is in many ways best thought of as a {\em mesoscopic} theory. We can distinguish between two extreme scales; an {\em outer} scale, characteristic of the macroscopic dynamics of the system producing the shock, and an {\em inner} microscopic scale characteristic of the collisionless plasma processes operating in the shock itself. DSA theory operates on intermediate scales located between these two extremes, much in the way that the Kolmogorov theory of turbulence operates between an outer driving scale and an inner dissipation scale.
Just as the simple power-law of Kolmogorov turbulence is only clearly seen when the inner and outer scales are sufficiently separated, so the simple power-law of test-particle DSA is only seen when the two scales are far enough apart. This can be a problem for heliospheric shocks, which are relatively small and weak, but for typical SNR parameters the scale separation can easily be ten decades or more.
DSA belongs to the general class of Fermi acceleration mechanisms, defined by the fact that the energy source driving the acceleration is differential motion of a magnetised plasma. It is distinguished from other Fermi processes by the fact that it utilises the strong compression in a shock front to drive the acceleration, and this leads to a number of unique properties. For a very physical account of the mechanism one can do little better than read Bell's original papers. The following argument gives a quick short-cut to the main results.
\subsubsection{Elementary test-particle theory}
The propagation of energetic charged particles embedded in a magnetised moving plasma can be described, under rather general conditions, by the advection diffusion equation
\begin{equation}
{\partial f\over\partial t} + \vec U \cdot\nabla f -
{1\over 3} \nabla\cdot \vec U p{\partial f\over\partial p} = \nabla\left(\kappa \nabla f\right)
+{1\over p^2}{\partial\over\partial p}\left(p^2 D {\partial f\over\partial p}\right).
\end{equation}
The various terms in this equation represent, first a standard convective derivate, secondly the effect of adiabatic compression, thirdly spatial diffusion with coefficient $\kappa$, and fourthly diffusion in momentum space with coefficient $D$. The basic assumption is that the energetic particles are scattered off irregularities in the magnetic field sufficiently rapidly that their distribution can be taken as close to isotropic in the local fluid frame and thus the full phase space distribution function reduces, in momentum space, to a scalar function of $p=\left | \vec p\right |$.
\begin{equation}
f(\vec x, \vec p, t) \to f(\vec x, p, t).
\end{equation}
This requires that the particle speed be very much greater than the typical fluid speed, which is certainly the case for relativistic particles in a non-relativistic flow. Under the same conditions the remaining transport processes are well represented, at the macroscopic level, as diffusive fluxes driven by simple gradient terms. The adiabatic compression term can be considered as a simple consequence of Liouville's theorem; if the system is compressed in physical space, it has to respond by expanding in momentum space, and the factor of one third comes about because under the isotropy assumption it expands equally in all three momentum dimensions. The momentum space diffusion term represents classical Fermi acceleration is normally small; we will neglect it from here on and retain only the spatial diffusion (this amounts to assuming that the flow speeds are fast relative to the local Alfven speed).
It is useful to look at the total acceleration flux of particles up through a given energy (or momentum) level,
\begin{equation}
\Phi(p) = \int {4\pi p^3\over 3} f(p) (-\nabla\cdot\vec U) \, d^3x.
\end{equation}
In the simple case of a plane shock all the compression is concentrated in a delta distribution located at the shock front and we deduce that there is a localised flux of particles upwards in momentum space at the shock front of magnitude
\begin{equation}
\Phi(p) = {4\pi p^3\over 3} f_0(p) \left(U_1-U_2\right)
\end{equation}
where $f_0(p)$ is the isotropic part of the phase space density of accelerated particles at the shock and $U_1$ and $U_2$ are the up-stream and down-stream velocities of the plasma with respect to the shock. We can now write down a simple particle conservation equation for the number of particles in the ``acceleration region'', which we take to be a box extending one diffusion length up-stream and down-stream of the shock,
\begin{equation}
L \approx {\kappa_1\over U_1} + {\kappa_2\over U_2}.
\end{equation}
This equation is
\begin{equation}
{\partial\over\partial t}
\left[4\pi p^2 f_0(p) L\right]
+{\partial\Phi\over\partial p}
= -4\pi p^2 f_0(p) U_2 +Q(p)
\end{equation}
and simply says that the number of particles in the acceleration region changes as a result of the divergence of the acceleration flux in momentum space, injection from a source $Q(p)$, and advection away downstream at the downstream flow velocity $U_2$. Substituting for $\Phi$ one finds,
\begin{equation}
L{\partial f\over\partial t}
+{U_1-U_2\over 3} p{\partial f\over\partial p} + U_1 f = {Q\over 4 \pi p^2}
\end{equation}
from which all the elementary results of DSA follow. We see immediately that in a steady state and away from injection
\begin{equation}
{U_1-U_2\over 3} p{\partial f\over\partial p} + U_1 f=0
\end{equation}
and thus the spectrum is a simple power law
\begin{equation}
f(p) \propto p^{-3U_1/(U_1-U_2)}
\end{equation}
with an exponent determined purely by the shock compression. We can also see that the acceleration time scale is simply
\begin{equation}
t_{\rm acc} = {3 L\over U_1 - U_2}
\end{equation}
and the equation can be written in characteristic form as
\begin{equation}
t_{\rm acc} {\partial f\over\partial t}
+ p {\partial f\over\partial p} =
- {3 U_1\over U_1-U_2} + {3Q\over 4 \pi p^2 (U_1 - U_2)}
\end{equation}
Note that while this ``box'' model equation is a very useful approximation and captures much of the physics it is not exact. Its major defect is that it artificially sharpens all spectral features.
\subsubsection{Nonlinear theory}
This elementary theory in which the accelerated particles are treated as test particles is clearly inadequate for cases where, as in GCR production by SNR shocks, the accelerated particles are assumed to carry a significant part of the energy dissipated in the shock. On the intermediate scales, what appears at the outer scale to be an unresolved shock, in the nonlinear version of DSA acquires an upstream precursor where the incoming flow is decelerated by the rising accelerated particle pressure gradient before a weakened subshock which is in turn resolved at the inner scale of the plasma kinetic level. This is shown schematically in Fig.~1.
Over the last twenty years or so a great deal of effort has been put into attempting to understand how this modification of the shock on intermediate scales affects, on the one hand the acceleration, and on the other the bulk dynamics of the flow. Early work was based mainly on two-fluid models, now largely of historical and didactic interest. The main thrust of more recent studies has been either to use semi-analytic approaches (Eichler, Malkov, Blasi, Gabici and co-workers), Monte-Carlo simulations (Ellison, Jones, Baring and co-workers) or full-scale numerical studies (Dorfi, Duffy, Jones, Kang and co-workers); for recent reviews see \cite{MD2, DB}. The good news is that there is now reasonable agreement between all these approaches and that in this sense the problem of nonlinear DSA can be said to be `understood'.
The starting point for the semi-analytic theories is the observation that on the intermediate scales we can treat the shock as a quasi-stationary structure with fixed mass and momentum fluxes,
\begin{eqnarray}
\rho U &=& A\\
A U + P_G + P_C &=& B
\end{eqnarray}
where $\rho$ is mass density, $U$ is flow speed, $P_G$ is the thermal pressure and
\begin{equation}
P_C = \int 4\pi p^2 {p v\over 3} f(p) dp
\end{equation}
is the pressure of the non-thermal particle population. We can still write down a particle conservation equation,
\begin{equation}
{\partial \Phi\over \partial p} = - 4\pi p^2 f_0(p) U_2
\end{equation}
however $\Phi$ is now distributed over the entire precursor region and is no longer localised at one point.
The key idea, as emphasised by Malkov in his seminal paper \cite{M97}, is that if we make some {\it Ansatz} relating the distribution throughout the precursor, $f(x,p)$, to the distribution at the subshock, $f_0(p)$, the momentum balance and number conservation equations become two coupled integro-differential equations for the two functions $f_0(p)$ and $U(x)$. The simplest such assumption, which leads to a remarkably simple set of ordinary differential equations, was originally proposed by Eichler and states that all particles penetrate a fixed distance $h$ upstream which is an increasing function of momentum and have a flat distribution downstream of this distance;
\begin{equation}
f(x,p) = \begin{cases}
f_0(p) & x> -h(p),\\
0 & x< -h(p).
\end{cases}
\end{equation}
Another more natural {\it Ansatz} would be to assume an exponential profile of the type familiar from test-particle theory,
\begin{equation}
f(x,p) = \begin{cases}
f_0(p) \exp\int_0^x {U(x')\over\kappa(x',p)}dx'
& x<0,\\
f_0(p) & x > 0).
\end{cases}
\end{equation}
This is quite similar to that originally proposed by Malkov,
\begin{equation}
f(x,p) = f_0(p) \exp\int \left(-{1\over3}{\partial \ln f_0\over \partial \ln p}\right){U(x)\,dx\over\kappa(x,p)}
\end{equation}
which can be motivated as follows.
The problem is that the compression is distributed over the precursor and no longer concentrated at the shock. Let us consider the extreme case where we distribute the compression uniformly everywhere so that we have a linear and convergent velocity field which we can take to have the form $U=-x$. Because particles then see the same compression no matter where they are the spatial diffusion and the acceleration decouple and it is easy to see that there must be a simple analytic solution. All particles are accelerated with
\begin{equation}
\dot p = {1\over 3} p
\end{equation}
and thus the total number of particles at any given momentum level must scale as $p^{-3}$ to give a constant acceleration flux $4\pi p^2 N(P) \dot p$. It is also clear, from the Gaussian nature of the Green's function for diffusion, that for a diffusion coefficient that is a function only of $p$ the spatial distribution of the particles should have a Gaussian profile and thus we can write,
\begin{equation}
f(x,p) = p^{-3} {1\over\sqrt{2\pi\sigma}}\exp\left(-{x^2\over 2\sigma}\right).
\end{equation}
It is easy to verify that this identically satisfies
\begin{equation}
{\partial f\over\partial t} + U {\partial f \over\partial x}
-{1\over 3}{\partial U\over\partial x} p{\partial f\over\partial p}
= {\partial\over\partial x} \left(\kappa {\partial f\over\partial x}\right)
\end{equation}
for $U=-x$ if
\begin{equation}
\sigma + {1\over 6} p {\partial\sigma\over\partial p} = \kappa.
\end{equation}
For the special case of $\kappa\propto p$
we get
\begin{eqnarray}
f(x,p) &=& p^{-7/2} \exp\left( -{7\over 6}{ x^2\over 2 \kappa}\right)\\
&=& f_0(p) \exp\int \left(-{1\over3}{\partial \ln f_0\over \partial \ln p}\right){U(x)\,dx\over\kappa(x,p)}
\end{eqnarray}
and thus Malkov's {\it Ansatz} is exact for a Bohm scaling of the diffusion coefficient and in the extreme limit where the compression is uniformly distributed.
A further improvement, due to Blasi, is to propose an {\it Ansatz} which interpolates naturally between Malkov's form, appropriate to strong modification, and the exponential form, which must be the test-particle limit,
\begin{equation}
f(x,p) = f_0(p) \exp\int \left(-{1\over3}{\partial \ln f_0\over \partial \ln p}\right)
\left(1 - {1\over r_{\rm tot}}\right)
{U(x)\,dx\over\kappa(x,p)}
\end{equation}
where $r_{\rm tot}$ is the total compression in the structure.
However there is a fundamental problem. All these approaches have to make the assumption that the mesoscopic structure of the shock can be treated as quasi-stationary. Unfortunately we have very good reasons for thinking that this is not the case. There exists a whole zoo of instabilities which can make the intermediate scale shock structure non-steady and even turbulent. Perhaps the best one can hope for is that even though the true structure is time-variable and fluctuating, in a time-averaged sense it may remain close to the quasi-static models assumed by the non-linear theories.
Actually this is not so implausible; one remarkable aspect of diffusive shock acceleration which impresses everyone working in the field is just how robust the fundamental physical process is.
The news is not all bad however. One consequence of the various instabilities is that one has the possibility to substantially amplify the local magnetic field in the shock precursor. This is highly desirable because the maximum energy to which one can accelerate is essentially given just by the magnetic field strength times the product of the outer velocity and length scales (this follows from very simple dimensional arguments; detailed studies of time-dependent acceleration introduce numerical factors which reduce this number by a factor of order 10, but leave the scaling intact). With conventional ISM and SNR parameters it has been known since the work of Lagage and Cesarsky \cite{LC} that it is difficult (though not impossible) to accelerate to the `knee' region of the GCR spectrum at about $3\times 10^{15}\,\rm eV$ and that the acceleration for all species should more naturally cut-off at a rigidity of about $10^{14} \,\rm V$ (unless one uses ideas such as significant cross-field diffusion in nearly perpendicular shocks \cite{J}).
\subsection{Magnetic Field Amplification}
The idea that the magnetic field could be substantially amplified by instabilities and the maximum energy thereby increased was first seriously advanced in quantitative form by Bell and Lucek \cite{LB, BL} although there are a number of anticipations of the idea in the literature \cite{McKV, VMcK}. Bell initially used a standard plasma dispersion analysis to show that in addition to the well-known resonant excitation of Alfven waves by cosmic ray streaming, under conditions appropriate to a shock precursor there should also be a very strong non-resonant instability, a result confirmed by subsequent analyses \cite {R, Z}. The, at first sight rather surprising, explanation of this instability is that it is essentially driven by the current carried by the cosmic rays. The high-energy protons (and heavy ions) accelerated at the shock penetrate into the upstream medium and, being positively charged, represent a positive cosmic ray current flowing into the upstream region. Of course the plasma reacts to preserve charge neutrality by pulling some low-energy thermal electrons with the cosmic rays so that there is an almost equal and opposite `return current'. The key point is that the carriers of the return current, being predominantly low energy electrons, are much more strongly magnetised than the high-energy particle that drive the current in the first place and form part of the MHD plasma. In discussing the dynamics of MHD systems we normally rewrite the $\vec j \times \vec B$ force in terms of magnetic pressure and field line tension by using
\begin{equation}
\nabla\times \vec B = {1\over\mu_0}\vec j
\end{equation}
where $\vec j$ is the total current. If we decompose this into a cosmic ray part and a `thermal' component,
\begin{equation}
\vec j = \vec j_{\rm CR} + \vec j_{\rm th}
\end{equation}
then it is clear that the Lorentz force acting on the MHD plasma is in fact
\begin{equation}
\vec j_{\rm th} \times \vec B = \left( \mu_0 \nabla\times\vec B - \vec j_{\rm CR}\right)\times \vec B
\end{equation}
so that if we treat the background plasma as a MHD system there is an additional force acting on it which is just the reaction to the Lorentz force exerted on the cosmic rays. On scales where the cosmic rays are scattered this force can be shown to reduce to the usual cosmic ray pressure gradient,
\begin{equation}
\left<-\vec j_{\rm CR} \times \vec B\right> \approx \nabla P_{\rm CR}
\end{equation}
but on scales where they are unmagnetised and stream freely it appears as a simple current driven term. Thus on scales smaller than the gyroradius of the driving particles one can in fact greatly simplify the analysis and just think of driving a constant current through the plasma. Clearly this will cause instabilities, and it is not hard to see that one expects the field lines to coil up into helical structures and expanding loops \cite{Z}. It is of course true that the total number of high-energy particles is small, and for a standard ``equal energy per decade'' $N(E)\propto E^{-2}$ the number per decade falls as $E^{-1}$. However the distance over which the current is driven grows with energy at least as fast as $E$ and thus the potential work done against the selfinductance of the circuit does not decrease.
Observationally there is now good evidence for such amplified fields. Actually this is not such a recent idea, Cowsik and Sarkar \cite{CS} inferred amplified fields of at least $80\,\mu\rm G$ in Cas A as far back as 1980, but the really convincing evidence has come from recent high-resolution X-ray observations. These have shown that in essentially all young historical supernova remnants, eg \cite{Bamba}, there are very narrow limb-brightened rims which exhibit a hard non-thermal X-ray spectrum best explained as synchrotron emission from ultra-relativistic electrons. The extreme narrowness of these rims (and the associated recent observation by Uchiyama et al \cite{U} of time-variability on scales of one year, so that the structures are small scale in both space and time) is most naturally explained by strong magnetic fields causing both fast acceleration of the electrons and rapid synchrotron cooling downstream of the shock. For a good review see Vink \cite{V}. It is of course true that the observations mainly show evidence for strong fields downstream of the shock, and as pointed out by Jokipii such fields can also be generated by the strong post-shock vorticity induced if the shock is propagating in a lumpy medium full of small scale density fluctuations without any need for upstream instabilities \cite{GJ}.
There is one obvious problem with claiming that this field amplification leads to acceleration to higher energies in that, if the above return current driven explanation is correct, the instability generates power in the field fluctuations only on scales smaller than the gyroradius of the driving particles. Malkov and Diamond \cite{MD} attempt to address this by invoking an inverse cascade driven by scattering off large scale inhomogeneities in the flow to transfer power to the larger scales. Another possibility is that the nonlinear interaction of the expanding current loops may produce regions of sufficiently concentrated field that the most energetic particles are brought back into the scattering regime (the gyroradius drops in the amplified field, diffusive acceleration becomes efficient and particles have to be accelerated to higher energy to escape in a boot-strap process).
\section{Propagation}
A major open problem is the precise nature of Galactic propagation. There are at least two very uncomfortable issues. One is the question of whether or not reacceleration is important. In reacceleration models for propagation the inferred production spectrum for the cosmic rays is quite soft with a differential energy exponent of -2.3 or so whereas without reacceleration it is harder with an exponent of around -2.0. The latter fits much better with theoretical expectations from acceleration theory, and the softer index is hard to reconcile with our current understanding of diffusive shock acceleration.
The other serious issue is the difficulty of reconciling the rather fast escape from the Galaxy inferred at low energies from secondary to primary rations (Boron to Carbon, sub-Iron to Iron etc) with the very low and essentially energy-independent anisotropy observed at around 10 to 100 TeV. This is substantially easier in reacceleration models, but even there it is not easy.
The current reference standard, the GALPROP model, is basically a computational version of the Ginzburg diffusion model from the 1960s with a few additions such as reaccleration. The only serious alternative, which seems to me to be worth much more attention that it has received to date, is the dynamical model of Zirakashvilli et al \cite{Z96, P97} which invokes an extended Galactic outflow with an embedded wound up spiral magnetic field.
\section{Conclusions}
There is a real need for an alternative model to GALPROP which in my view should be a dynamical model along the lines of Zirakashvili et al. If nothing else it is dangerous to use only one standard model, especially when there are a number of serious issues surrounding the basic assumptions of GALPROP (essentially static diffusion in a small stationary halo). With the much improved data becoming available from Fermi and possibly new radio data from Lofar (and in the future SKA) the time is ripe for such a development.
\section{Acknowledgments}
This note was written while a participant in the Kavli Institute for Theoretical Physics programme Astroplasmas09 supported in part by the National Science Foundation under Grant No. PHY05-51164. Stimulating discussions with the other participants are gratefully acknowledged.
\begin{thebibliography}{99}
\bibitem{Westerlund2}
Aharonian et al (HESS collaboration) (2007)
A\&A 467 1075-1080. Also
arXiv:astro-ph/0703427
\bibitem{ALS}
W. I. Axford, E. Leer and G. Skadron,
Proc. 15th ICRC (Plovdiv), {\bf 11} 132-137 (1977).
\bibitem{Bamba}
A. Bamba et al, {\it Ap.J.} {\bf 621} 793--802 (2005)
\bibitem{Ba}
A. R. Bell, {\it MNRAS} {\bf 182} 147-156 (1978)
\bibitem{Bb}
A. R. Bell, {\it MNRAS} {\bf 182} 443-455 (1978)
\bibitem{BL}
A. R. Bell and S. G. Lucek, {\it MNRAS} {\bf 321} 433--438 (2001).
\bibitem{BO}
R. D. Blandford and J. P. Ostriker, {\it Ap. J.} {\bf 221} L29--32 (1978)
\bibitem{CS}
R. Cowsik and S. Sarkar, {\it MNRAS} {\bf 191} 855--861 (1980)
\bibitem{DB}
P. Duffy and K. Blundell, {\it Plasma Physics and Controlled Fusion}, {\bf 47} B667-B678 (2005).
\bibitem{GJ}
J. Giacalone and J. R. Jokipii, {\it Ap.J.} {\bf 663} L41--44 (2007)
\bibitem{GS}
V. L. Ginzburg and S. I. Syrovatskii, {\it The origin of cosmic rays}, Gordon and Breach, New York, 1969.
\bibitem{Hoyle} Hoyle, F, (1960) MNRAS 120 338.
\bibitem{J}
J. R. Jokipii, {\it Ap.J.} {\bf 313} 842--846 (1987)
\bibitem{K}
G. F. Krymskii, {\it Akademiia Nauk SSSR, Doklady} \bf 234 \rm 1306--1308 (1977)
\bibitem{LC}
P. O. Lagage and C. J. Cesarsky, {\it A\&A} {\bf 125} 249--257 (1983).
\bibitem{LB}
S. G. Lucek and A. R. Bell, {\it MNRAS} {\bf 314} 65--74 (2000)
\bibitem{M97}
Malkov, M. (1997) {\it Ap.J.} {\bf 485} 638.
\bibitem{MD2}
M. A. Malkov and L. O'C. Drury, {\it Reports on Progress in Physics}, {\bf 64}, 429-481 (2001)
\bibitem{MD}
M. A. Malkov and P. H. Diamond, {\it Ap.J.} {\bf 654} 252--266 (2007)
\bibitem{McKV}
J. McKenzie and H. J. V\"olk
Proc 17th ICRC (Paris) {\bf 9} 242--245 (1981).
\bibitem{P97}
Ptuskin, V. S.; Voelk, H. J.; Zirakashvili, V. N.; Breitschwerdt, D.
(1997)
Astronomy and Astrophysics, v.321, p.434-443
\bibitem{R}
B. Reville, J. G. Kirk, P. Duffy and S. O'Sullivan, {\it A\&A} {\bf 475} 435--439 (2007).
\bibitem{U}
Y. Uchiyama et al, {\it Nature} {\bf 449} 576--578 (2007)
\bibitem{VMcK}
H. J. V\"olk and J. McKenzie
Proc 17th ICRC (Paris) {\bf 9} 246--249 (1981).
\bibitem{V}
J. Vink, arXiv:astro-ph/0601131 (2006)
\bibitem{Z}
V. N. Zirakashvili, V. S Ptuskin and H. J. V\"olk, arXiv0801.4486 (2008)
\bibitem{Z96}
Zirakashvili, V. N.; Breitschwerdt, D.; Ptuskin, V. S.; Voelk, H. J. (1996)
Astronomy and Astrophysics, v.311, p.113-126
\end{thebibliography}
\end{document}