Propagating Cosmic Rays with exact Solution of Fokker-Planck Equation

TL;DR

This paper presents an exact solution to the Fokker-Planck equation for cosmic ray propagation, enabling accurate modeling across ballistic, diffusive, and transdiffusive regimes, improving interpretation of cosmic ray spectra.

Contribution

It introduces a simplified, exact propagator for the Fokker-Planck equation, enhancing cosmic ray transport models beyond traditional Gaussian approximations.

Findings

Provides a complete description of cosmic ray propagation regimes.

Derives a simplified exact propagator for practical use.

Improves interpretation of cosmic ray spectral anomalies.

Abstract

Shortfalls in cosmic ray (CR) propagation models obscure the CR sources and acceleration mechanisms. This problem became particularly obvious after the Fermi, Pamela, and AMS-02 have discovered the electron/positron and $p/$He spectral anomalies. Most of the CR models use diffusive propagation that is inaccurate for weakly scattered energetic particles. So, some parts of the spectra affected by the heliospheric modulation, for example, cannot be interpreted. I discuss and adopt an exact solution of the Fokker-Planck equation arXiv1610.01584, which gives a complete description of a ballistic, diffusive and transdiffusive (intermediate between the first two) propagation regimes. I derive a simplified version of an exact Fokker-Planck propagator that can easily be employed in place of the Gaussian propagator, currently used in major Solar modulation and other CR transport models.

Equations

$$\frac{\partial f}{\partial t}+v\mu \frac{\partial f}{\partial x}=\frac{\partial }{\partial \mu }\left(1-{\mu }^2\right)D\left(\mu ,E\right)\frac{\partial f}{\partial \mu }.$$

$$\frac{\partial f_0}{\partial t}-{\kappa }_2\frac{{\partial }^2{f}_0}{\partial x^2}=-{\kappa }_4\frac{{\partial }^4{f}_0}{\partial x^4}+{\kappa }_6\frac{{\partial }^6{f}_0}{\partial x^6}-\dots ,$$

$$D\left(E\right)t\to t,\frac{D}{v}x\to x$$

$$\frac{\partial f}{\partial t}+\mu \frac{\partial f}{\partial x}=\frac{\partial }{\partial \mu }\left(1-{\mu }^2\right)\frac{\partial f}{\partial \mu }$$

$$M_{ij}\left(t\right)=⟨{\mu }^i{x}^j⟩={\int }_{-\infty }^{\infty }dx{\int }_{-1}^1{\mu }^i{x}^{j}fd\mu /2$$

$$M_{ij}\left(t\right)=M_{ij}\left(0\right)e^{-i\left(i+1\right)t}+{\int }_0^t{e}^{i\left(i+1\right)\left(t^{\prime }-t\right)}\times$$

$$\left[j{M}_{i+1,j-1}\left(t^{\prime }\right)+i\left(i-1\right)M_{i-2,j}\left(t^{\prime }\right)\right]d{t}^{\prime }$$

$$f_0\left(x,t\right)={\int }_{-1}^{1}f\left(\mu ,x,t\right)d\mu /2,$$

$$M_j\equiv M_{0,j}$$

$$f_{\lambda }\left(t\right)={\int }_{-\infty }^{\infty }{f}_0\left(x,t\right)e^{\lambda x}dx=\sum _{n=0}^{\infty }\frac{{\lambda }^{2n}}{\left(2n\right)!}{M}_{2n}\left(t\right)$$

$$f_0\left(x,t\right)\approx \frac{1}{4y}\left[\mathrm{erf}\left(\frac{x+y}{\Delta }\right)-\mathrm{erf}\left(\frac{x-y}{\Delta }\right)\right].$$

$$M_2=\int x^2{f}_0\left(x,t\right)dx,M_4=\int x^4{f}_0\left(x,t\right)dx$$

$$M_2=\frac{t}{3}-\frac{1}{6}\left(1-e^{-2t}\right)$$

$$M_4=\frac{1}{270}{e}^{-6t}-\frac{t+2}{5}{e}^{-2t}+\frac{1}{3}{t}^2-\frac{26}{45}t+\frac{107}{270}$$

$$y={\left[\frac{45}{2}\left(M_2^2-\frac{1}{3}{M}_4\right)\right]}^{1/4}$$

$$\Delta =\sqrt{2M_2-\sqrt{10}\sqrt{M_2^2-\frac{1}{3}{M}_4}}$$

$$f_0\left(x,t\right)=\sqrt{\frac{3}{2\pi t}}{e}^{-3x^2/2t}$$

Full text

cosmic rays\sepsupernova remnants \sepcosmic ray propagation \sepinterstellar medium\PACS98.38.Mz \sep98.70.Sa \sep96.50.S- \sep95.35.+d

I Lacuna in CR Transport Model

The cosmic ray (CR) propagation in turbulent environments, such as the interstellar medium (ISM) or Heliosphere, has been actively researched for more than half a century Jok1971RvGSP . Time asymptotically, CRs propagate diffusively; after several collisions, they “forget” their initial velocities and enter a random walk process. However, in astrophysical objects with infrequent particle collisions, there may not be enough time or room for even a few collisions. In such systems, the focus shifts to earlier propagation phases, which are better described as ballistic rather than diffusive propagation. The question is, what is in between these two regimes and how long it lasts?

The transition from ballistic to diffusive transport regime has always been a challenge for the theory. At the same time, it is often the key to understanding the CR sources. Since the particle mean free path (m.f.p) usually grows with energy, some part of their spectrum almost inescapably falls into a transient category where neither ballistic nor diffusive approximation applies. I will call this regime transdiffusive and argue that it lasts for long enough to compromise both the ballistic and diffusive model predictions. During this propagation phase, CR protons accelerated in supernova remnants (SNR), for example, may reach a nearby molecular cloud, making themselves visible by interacting with its dense gas MDS_11NatCo ; AgileW44_11 . The CR protons of lower energies would instead be diffusively confined to the SNR shell and evade detection. Due to a high CR intensity near the source, however, their confinement here must be due to self-generated Alfven waves. At a minimum, this problem should be treated at a quasilinear level MetalEsc13 , as opposed to the linear CR transport, considered throughout this paper. Another example is the propagation of solar energetic particles to 1 AU. Also, in this case, the m.f.p. of some particles is comparable to, or even exceed, 1 AU, so neither the diffusive nor ballistic approximation applies Klassen2016A&A ; Bian2017ApJ .

Galactic CRs ultimately propagating through the Heliosphere to the observer cannot always be propagated back to their source within simple diffusion or ballistic paradigms, so their spectra cannot be fully understood. This problem is particularly relevant to striking anomalies in the CR spectra and composition, which are becoming a general trend in the CR observations. Besides the $e^{+}/e^{-}$ anomaly, there is a $\sim 0.1$ difference in rigidity indices of proton and He. Although the explanations are available (see, e.g., MDSPamela12 ; MalkPositrons2016 , and a companion paper in this volume), the low-energy parts of these spectra are strongly affected by the solar modulation. Curvature and gradient drifts in the Heliospheric magnetic field are mostly treated by considering particle propagation along the field line as diffusive, e.g. Potgieter2013 , which we will show to be inaccurate for sufficiently energetic particles with long m.f.p.

II Governing Equation

The Fokker-Planck (FP) equation is a minimalist model suitable for the CR transport. An ambient magnetic field justifies a 1D treatment, while its fluctuating part supports the particle scattering in pitch angle. The simplest form of FP equation for the CR distribution function $f$ is the following:

[TABLE]

Here $x$ is directed along the local magnetic field, $\mu$ is the cosine of the particle pitch angle, $v,E$ are the particle velocity and energy, conserved in interactions with quasi-static magnetic turbulence. $D$ is the scattering rate (collision frequency).

One propagation scenario that Eq.(1) describes very well comes about through an instant release of a cloud of particles into a scattering medium. Again, Galactic SNRs, widely believed to generate CRs with energies up to $\sim 10^{15}$eV, must accelerate them in SNR shock waves with a subsequent release into a turbulent ISM. The question then is how exactly the particle density (the isotropic component of $f$) propagates along a magnetic flux tube that intersects the SNR shell. The goal is to achieve the simplicity of diffusive description (e.g., Jokipii66 and below) which is a well-known derivative of Eq.(1). As emphasized earlier, the diffusive treatment is inadequate in the preceding ballistic and transdiffusive propagation phases, while the latter is often the key for probing into the source.

II.1 Restricting Propagation Models by Limiting Cases

Because of the difficulties in reducing the FP equation to a manageable isotropic form, a framework for such reduction limited by the extreme cases of ballistic and diffusive propagation is helpful. We derive both regimes directly from Eq.(1), by eliminating angular dynamics.

In the ballistic case, which strictly applies to times shorter than the collision time $t\ll t_{c}\sim 1/D$, one can neglect the r.h.s. altogether. The solution then follows from integrating along the particle trajectories, $x-\mu vt=const$ (Liouville’s theorem), with a conserved pitch angle, $\mu=const$. The solution is simply $f\left(x,\mu,t\right)=f\left(x-v\mu t,\mu,0\right)$.

Consider an isotropic point source: $f\left(x,\mu,0\right)=\nicefrac{{1}}{{2}}\delta\left(x\right)\Theta\left(1-\mu^{2}\right),$ where $\delta$ and $\Theta$ denote the Dirac’s delta and Heaviside unit step functions, respectively. From the above solution for $f\left(x,\mu,t\right)$, one obtains the ballistic expansion in form of the second moment, $\left\langle x^{2}\right\rangle=v^{2}t^{2}/3$ by integrating $x^{2}f=\nicefrac{{1}}{{2}}x^{2}\delta\left(x-v\mu t\right)\Theta\left(1-\mu^{2}\right)$ over $x$ and $\mu$. The result describes a free escape with the mean square velocity $v/\sqrt{3}$, while the maximum particle velocity (along $x$) is $v$. The pitch angle averaged particle distribution, $f_{0}\left(x,t\right)=\left(2vt\right)^{-1}\Theta\left(1-x^{2}/v^{2}t^{2}\right)$, is best described as an expanding ’box’ of decreasing height. Among earlier attempts to reduce $f$ to its pitch angle-averaged part, $f\left(t,x,\mu\right)\to f_{0}\left(t,x\right)$, an approach leading to a “telegraph” equation, can be readily tested using the above box solution. We will briefly discuss this approach later and show that it is inconsistent with the ballistic limit of $f_{0}$ obtained directly from the FP equation. Needless to say that the exact solution of Eq.(1), presented further in this paper, converges to the above-described box distribution at $t\ll t_{c}$.

The second, well studied propagation regime is diffusive. It dominates at $t\gg t_{c}\sim 1/D$ and is treated in a way opposite to the above-described ballistic regime, Jok1971RvGSP . The r.h.s. of Eq.(1) is now the leading term, thus implying that the particle distribution is close to isotropy, $\partial f/\partial\mu\to 0$. Working to higher orders in anisotropic corrections $\sim 1/D$, and averaging the equation over $\mu$, one obtains the following equation for $f_{0}\left(x,t\right)$ MS2015

[TABLE]

with $\kappa_{2n}\sim 1/D^{n}$. The last equation results from an asymptotic (Chapman-Enskog) expansion of the problem in $1/D$ under the scattering symmetry: $D\left(-\mu\right)=D\left(\mu\right)$. It is valid only for $t\gg t_{c}\sim 1/D$, and all the short-time-scale, ballistic propagation effects are intentionally eliminated (cf. elimination of secular terms in perturbative treatments). A failure to do so results in a second order time derivative in Eq.(2) (already mentioned telegraph term) which is illegitimate unless $t\gg t_{c}$. Nevertheless, the telegraph equation has been putting forward over the last 50 years as a viable tool for describing the CR propagation from the ballistic to diffusive phases.

Meanwhile, the r.h.s. of eq.(2) provides small hyperdiffusive corrections that may be omitted at $t>t_{c}$, as the higher spatial derivatives quickly decay because of the smoothing effect from the diffusive term on its l.h.s. These corrections do not shed much light on the ballistic and transdiffusive propagation regimes, probably unless the series is summed up with no truncation. The latter requirement derives from a method whereby an exact solution of the parent FP equation 2016arXiv161001584M is evaluated. The evaluation consists in summing up an infinite series of moments $\left\langle x^{2n}f_{0}\left(x,t\right)\right\rangle$ that are evidently connected with the infinite series of coefficients $\left\{\kappa_{2n}\right\}$ in eq.(2). Conversely, by including just one (or several) hyperdiffusive correction outside of their validity range, $t\gg t_{c}$, one may even decrease the accuracy of the diffusive approximation. It can also be shown MS2015 that within its validity range, a truncated version of eq.(2), with $\kappa_{2n}=0$ for $n>2$, can be mapped onto the telegraph equation. It follows that neither a truncated hyperdiffusive approach nor the telegraph equation (a subset of the former) cannot adequately reproduce the FP solution at times shorter than $t\gg t_{c}$. This was recently demonstrated in Ref.Litvin2016PhPl , by a numerical integration of Eq.(1). The results of this work are illustrated for $t=t_{c}$ in Fig.1. We will quantify the constraint $t\gg t_{c}$, repeatedly stressed above, by comparing the full FP solution with its diffusive limit (see 2016arXiv161001584M for more details).

The primary failure of the diffusive approach is an unrealistically fast (acausal) propagation, which is especially pronounced during the ballistic and transdiffusive phases. Mathematically, the approximation violates an upper bound $\left|x\right|\leq vt$ that immediately follows from Eq.(1) for a point source solution, discussed above. There have been attempts to overcome this problem, but no adequate ab initio description of particle spreading that would cover ballistic and diffusive phases was elaborated. The most persistent such attempt is based on the telegraph equation discussed above. It has a misleading impact on the field of CR propagation for that simple reason that the solution of this equation is inconsistent with its parent FP equation. We obtained this simple result by considering the ballistic propagation phase directly from eq.(1) (see MS2015 ; M_PoP2015 ; 2016arXiv161001584M for more discussion).

It follows that there are no viable analytical tools to address the earlier phases of particle propagation, except to possibly sum up the series of hyperdiffusive terms or just to solve the FP equation directly. Below, we take the second option.

III Exact Solution of FP equation

The energy dependence of the particle scattering frequency enters eq.(1) only as a parameter, i.e., $D\left(E\right)$. The possible pitch-angle dependence of $D$ typically scales as $D\left(\mu\right)\propto\left|\mu\right|^{q-1}$ Jok1971RvGSP , thus being suppressed in an important case $q=1$, where $q$ is the power-law index of magnetic turbulence. Under these, quite realistic assumptions, the FP equation can be solved exactly 2016arXiv161001584M . To describe this solution, it is convenient to rewrite Eq.(1) using dimensionless time and length units according to the following transformations

[TABLE]

Instead of Eq.(1) we thus have

[TABLE]

This equation contains no parameters, thus precluding any direct asymptotic expansion in a small parameter, unless it enters the problem implicitly through the initial condition $f\left(x,\mu,0\right)$. In particular, if one is using Eq.(2) ($1/D$- type expansion), not only should the initial distribution be close to isotropy, but it should also be spatially broad. The latter condition will prevent a high anisotropy from arising via the second term on the l.h.s. of Eq.(4). Hence, the problem of a point source spreading (Green’s function, or fundamental solution) can not be treated using conventional $1/D$ expansion, until $f$ becomes quasi-isotropic, that is broadened to $x\gtrsim 1$.

The exact solution of Eq.(4) can be obtained using a fully resolvable infinite set of moments of $f\left(\mu,x\right)$

[TABLE]

for any integer $i,j\geq 0$. The lowest moment $M_{00}$ is automatically conserved by Eq.(4) (as being proportional to the number of particles) and we normalize it to unity, $M_{00}=1$. All the higher moments can be explicitly obtained from the following recurrence relation

[TABLE]

[TABLE]

Focusing on a point source (fundamental) solution, we assume the initial distribution $f\left(x,\mu,0\right)$ to be symmetric in $x$ and isotropic in $\mu$ which eliminates the odd moments. Furthermore, the initial spatial width must then also be set to zero, $M_{02}\left(0\right)=\left\langle x^{2}\right\rangle_{0}=0$.

From the mathematical point of view, only a full set (first two moments have been calculated by G.I. Taylor Taylor01011922 ) of moments in Eq.(6) provides a complete solution $f\left(x,\mu,t\right)$ of Eq.(4) given the initial value, $f\left(x,\mu,0\right)$ that determines the matrix $M_{ij}\left(0\right)$ in Eq.(6). Moreover, to adequately reproduce the ballistic and transdiffusive phases the series of moments cannot be truncated. Considering the fundamental solution, we will focus on the isotropic part of particle distribution

[TABLE]

as only this part contributes to the particle number density. To obtain the fundamental solution we impose the initial condition $f_{0}\left(x,0\right)=\delta\left(x\right)$. The matrix elements that represent $f_{0}$ are, therefore, $M_{0,j}$, which we denote $M_{j}$:

[TABLE]

Note, that $M_{ij}$ with $i>0$ are not small and remain essential for calculating the full set of the moments $M_{j}$. To link them to $f_{0}$, we use the moment-generating function

[TABLE]

where we omitted the odd moments irrelevant to the fundamental (symmetric in $x$) solution. The above expansion may be cast in a familiar Fourier transform of f_{0}$$\left(x,t\right) by setting $\lambda=-ik$.

Since expressions for the moments $M_{2n}$ are becoming cumbersome with growing $n$, an exact form the Green’s function $f_{0}\left(x,t\right)$, which can be recovered from eq.(8) by inverting the Fourier integral, is also not simple. Therefore, in the next section, we derive a new simplified version of the exact FP propagator that was recently obtained in Ref. 2016arXiv161001584M .

IV Two-moment Fokker-Plank Propagator

The infinite series entering the moment generating function $f_{\lambda}$ in eq.(8) has been summed up by considering the cases of small and large values of $t$ and $\lambda t$. Despite the multiplicity of limiting cases associated with these two independent quantities, all the expressions for the sum $f_{\lambda}\left(t\right)$ are surprisingly similar. They can be unified under a single *approximate (but valid *for all $x$ and $t$) expression for $f_{0}$:

[TABLE]

It has been obtained in 2016arXiv161001584M from an inverse Fourier transform, $f_{\lambda}\left(t\right)\mapsto f_{0}\left(x,t\right)$, after summing up the series for $f_{\lambda}$ in eq.(8). The two independent functions of time, $y\left(t\right)$ and $\Delta\left(t\right)$ can be expressed through the moment $M_{2}\left(t\right)$, which we calculate exactly from eq.(6). The solution $f_{0}$ with $y$ and $\Delta$ so obtained compares very well with the numerical FP solution. The disadvantage of this single-moment representation of $y$ and $\Delta$ is that it requires some (fairly minor, though) changes in $y\left(M_{2}\right)$ and $\Delta\left(M_{2}\right)$, between the cases $t\lesssim 1$ and $t\gtrsim 1$ 2016arXiv161001584M .

Here we suggest an alternative representation of the functions $y$ and $\Delta$. Although they lead to a slightly less accurate value of $f_{0}$ at $t\sim 1$ in eq.(9), but are the same for arbitrary $t:$ $0<t<\infty$. The idea behind this method of determination of $y$ and $\Delta$ is very simple. As the general form of the solution given in eq.(9) must arguably be the same for all $t$, we find the functions $y\left(t\right)$ and $\Delta\left(t\right)$ by requiring that $f_{0}\left(x,t\right)$ exactly satisfies the following two relations

[TABLE]

Recall, that we know exact values for all moments $M_{n}$ from eq.(6). Here, we will only use $M_{2}$ and $M_{4}$, which satisfy the initial conditions, $M_{2}\left(0\right)=0$ and $M_{4}\left(0\right)=0$:

[TABLE]

[TABLE]

Substituting $f_{0}$ from eq.( 9) into eqs.(10), we find

[TABLE]

[TABLE]

The FP solution, cast in a simplified form of eq.(9), is not more difficult than the familiar diffusive solution. If we ignore, for a moment, the time dependence of $y$ in the error functions, the FP solution appears as the solution of a conventional diffusion problem with an initial particle density evenly distributed between $-y<x<y$, and zero otherwise. The essential difference is only in the form of $y\left(t\right)$ and $\Delta\left(t\right)$. The first notable aspect of this solution is that at $t\ll 1$ it exactly corresponds to an ’expanding box’ ballistic solution described in Sec. II.1. Indeed, since $\Delta\propto t^{2}$ and $y\approx t$ for $t\ll 1$, the difference of the two error functions yields $2\Theta\left(1-x^{2}/t^{2}\right)$ and $f_{0}$ in eq.(9) is the same as the expanding box solution obtained in Sec.II.1. The telegraph solution, on the contrary, is inconsistent with this expansion regime as it contains two (nonexistent in the FP solution) singular components at the two propagating fronts, let alone the front positions and the overall profile, Fig1.

The width of the propagating fronts at $x=\pm y$, determined by $\Delta\left(t\right)$, behaves as follows, Fig.2. At small $t\ll 1$, when the box is expanding ballistically, i.e. $y\approx t$, the wall thickness $\Delta\approx 2t^{2}/3\sqrt{5}.$ After gradually proceeding through the transdiffusive phase, these quantities become $y\approx\left(11t/6\right)^{1/4}$ and $\Delta\approx\left(2t/3\right)^{1/2}$ for $t\gg 1$. Accordingly, the expression in eq.(9) converges (rather slowly, though) to:

[TABLE]

which is the diffusive asymptotic solution of the pitch angle averaged FP equation, given by eq.(2) with $\kappa_{2}=1/6$ and all the hyperdiffusive coefficients $\kappa_{2n}=0$ for $n>1$. Summarizing this section, the two-moment single formula representation of the FP solution in eq.(9) has correct asymptotic limits at $t\to 0,\infty$, both obtained independently. The remaining deviations from the numerical solution at $t\sim 1$ are minor and more than compensated by the simplicity of eq.(9) and its validity for all $0<t<\infty$.

V Conclusions

The exact solution of FP equation obtained in 2016arXiv161001584M is transformed into a simple form that accurately evolves the pitch angle averaged particle distribution $f_{0}\left(x,t\right)$, uniformly in $-\infty<x<\infty$ and $0\leq t<\infty$.

The overall CR propagation can be categorized into three phases: ballistic ($t<1)$, transdiffusive ($t\sim 1$) and diffusive ($t\gg 1$), (time in units of collision time $t_{c}$). In the ballistic phase, the source expands as a “box” of size $\Delta x\propto\sqrt{\left\langle x^{2}\right\rangle}\propto t$ with thickening “walls” at $x=\pm y\left(t\right)\approx\pm t$ of the width $\Delta\propto t^{2}$. The next, transdiffusive phase is marked by the box’s walls thickened to a sizable fraction of the box $\Delta\sim\Delta x\sim y$ and its slower expansion, Fig. 2. Finally, the evolution enters the conventional diffusion phase, in which $\Delta x\sim\Delta\propto\sqrt{t}$, while the walls are completely smeared out, as $y\propto t^{1/4}$, so $y\ll\Delta$.

In constraining earlier FP-based models for the CR propagation, the exact FP solution reveals the following:

• •

the conventional diffusion approximation can be safely applied but, only after 5-7 collision times, depending on the accuracy requirements

• •

a popular telegraph approach, originally intended to cover also the earlier propagation phases at $t\lesssim 1$, is inconsistent with the exact FP solution (see also 2016arXiv161001584M )

• •

no signatures of (sub) super-diffusive propagation regimes are present in the exact FP solution

The latter regimes are occasionally postulated, e.g., in studies of diffusive shock acceleration (DSA), in the form of a power-law dependence of particle dispersion $\sqrt{\left\langle x^{2}\right\rangle}\propto t^{\alpha}$, with $1/2<\alpha<1$ (superdiffusion) or $0<\alpha<1/2$ (subdiffusion). The exact FP propagation leads to $\sqrt{\left\langle x^{2}\right\rangle}$ that smoothly changes from the ballistic ($\alpha\to 1$) to diffusive ($\alpha\to 1/2$) propagation with no dwelling at any particular value of $\alpha$ between these limits. However, certain types of scattering fields in shock environments, e.g., MD06 , may result in both superdiffusive (Lévy flights) and subdiffusive (long rests) transport anomalies. Such fields are, however, less generic than those leading to an isotropic scattering considered in this paper. They should perhaps be justified on a case-by-case basis.

Acknowledgements

This work was supported by the NASA Astrophysics Theory Program under Grant No. NNX14AH36G.

?refname?

• (1)

M. A. Malkov, Exact solution of the Fokker-Planck equation for isotropic scattering arXiv:1610.01584, doi:doi:10.1103/PhysRevD.95.023007(2017).

• (2)

J. R. Jokipii, Propagation of cosmic rays in the solar wind., Reviews of Geophysics and Space Physics 9 (1971) 27–87.

doi:10.1029/RG009i001p00027.

• (3)

M. A. Malkov, P. H. Diamond, R. Z. Sagdeev, Mechanism for spectral break in cosmic ray proton spectrum of supernova remnant W44, Nature Communications 2/1195.

arXiv:1004.4714, doi:10.1038/ncomms1195.

• (4)

A. Giuliani, M. Cardillo, M. Tavani, et al., Neutral Pion Emission from Accelerated Protons in the Supernova Remnant W44, *Astrophys. J. * *Lett. *742 (2011) L30.

arXiv:1111.4868, doi:10.1088/2041-8205/742/2/L30.

• (5)

M. A. Malkov, P. H. Diamond, R. Z. Sagdeev, F. A. Aharonian, I. V. Moskalenko, Analytic Solution for Self-regulated Collective Escape of Cosmic Rays from Their Acceleration Sites, *Astrophys. J. *768 (2013) 73.

arXiv:1207.4728, doi:10.1088/0004-637X/768/1/73.

• (6)

A. Klassen, N. Dresing, R. Gómez-Herrero, B. Heber, R. Müller-Mellin, Unexpected spatial intensity distributions and onset timing of solar electron events observed by closely spaced STEREO spacecraft, Astronomy and Astrophys.593 (2016) A31.

doi:10.1051/0004-6361/201628734.

• (7)

N. H. Bian, A. G. Emslie, E. P. Kontar, The Role of Diffusion in the Transport of Energetic Electrons during Solar Flares, *Astrophys. J. *835 (2017) 262.

arXiv:1612.09456, doi:10.3847/1538-4357/835/2/262.

• (8)

M. A. Malkov, P. H. Diamond, R. Z. Sagdeev, Proton-Helium Spectral Anomaly as a Signature of Cosmic Ray Accelerator, Physical Review Letters 108 (8) (2012) 081104.

arXiv:1110.5335, doi:10.1103/PhysRevLett.108.081104.

• (9)

M. Malkov, P. Diamond, R. Sagdeev, Positive charge prevalence in cosmic rays: Room for dark matter in the positron spectrum, Physical Review D 94 (6) (2016) 063006.

• (10)

M. S. Potgieter, Solar modulation of cosmic rays, Living Reviews in Solar Physics 10 (1) (2013) 3.

doi:10.12942/lrsp-2013-3.

URL http://dx.doi.org/10.12942/lrsp-2013-3

• (11)

J. R. Jokipii, Cosmic-Ray Propagation. I. Charged Particles in a Random Magnetic Field, *Astrophys. J. *146 (1966) 480–+.

doi:10.1086/148912.

• (12)

M. A. Malkov, R. Z. Sagdeev, Cosmic Ray Transport with Magnetic Focusing and the "Telegraph" Model, *Astrophys. J. *808 (2015) 157.

arXiv:1502.01799, doi:10.1088/0004-637X/808/2/157.

• (13)

Y. E. Litvinenko, P. L. Noble, Comparison of the telegraph and hyperdiffusion approximations in cosmic-ray transport, Physics of Plasmas 23 (6) (2016) 062901.

doi:10.1063/1.4953564.

• (14)

M. A. Malkov, Cosmic ray confinement and transport models for probing their putative sources, Physics of Plasmas 22 (9) (2015) 091505.

arXiv:1503.01466, doi:10.1063/1.4928941.

• (15)

G. I. Taylor, Diffusion by continuous movements, Proceedings of the London Mathematical Society s2-20 (1) (1922) 196–212.

[

](http://arxiv.org/abs/http://plms.oxfordjournals.org/content/s2-20/1/196.full.pdf+html)arXiv:http://plms.oxfordjournals.org/content/s2-20/1/196.full.pdf+html,

doi:10.1112/plms/s2-20.1.196.

• (16)

M. A. Malkov, P. H. Diamond, Nonlinear Shock Acceleration Beyond the Bohm Limit, *Astrophys. J. *642 (2006) 244–259.

doi:10.1086/430344.