COSMIC RAY TRANSPORT WITH MAGNETIC FOCUSING AND THE "TELEGRAPH" MODEL

The problem addressed here is fundamental and not new to the cosmic ray (CR) transport studies. It can be formulated very plainly: how does one describe CR transport by only their isotropic component after the anisotropic one has decayed by scattering on magnetic irregularities? Suppose the angular distribution of CRs is given by the function $f(\mu ,t,z)$, obeying an equation from which the rapid gyro-phase rotation is already removed (drift approximation, e.g., Vedenov et al. 1962; Kulsrud 2005)

$$\begin{eqnarray}&&\displaystyle \frac{\partial f}{\partial t}+v\mu \displaystyle \frac{\partial f}{\partial z}=\displaystyle \frac{\partial }{\partial \mu }\left(1-{\mu }^{2}\right){\mathcal{D}}(\mu )\displaystyle \frac{\partial f}{\partial \mu }.\end{eqnarray} \tag{ 1 }$$

Here, z is the local coordinate along the ambient magnetic field, μ is the cosine of the particle pitch angle, and ${\mathcal{D}}$ is the pitch angle diffusion coefficient. Now, we take the next step in simplifying the transport description and seek an equation for the pitch-angle averaged distribution

$$\begin{eqnarray*}&&{f}_{0}(t,z)\equiv \displaystyle \frac{1}{2}{\int }_{-1}^{1}f(\mu ,t,z)d\mu \equiv \langle f\rangle .\end{eqnarray*}$$

The basic solution to this problem has been known for at least half a century (e.g., Jokipii 1966 and references therein). To the leading order in $1/{\mathcal{D}}$ (assuming the characteristic scale and time of the problem being longer than the particle mean free path and collision time) it can be obtained straightforwardly by averaging Equation (1)

$$\begin{eqnarray*}&&\displaystyle \frac{\partial {f}_{0}}{\partial t}=-\displaystyle \frac{v}{2}\displaystyle \frac{\partial }{\partial z}\left\langle\left(1-{\mu }^{2}\right)\displaystyle \frac{\partial f}{\partial \mu }\right\rangle,\end{eqnarray*}$$

and substituting $\partial f/\partial \mu \ll {f}_{0}$ from Equation (1), as

$$\begin{eqnarray}&&\displaystyle \frac{\partial f}{\partial \mu }\approx -\displaystyle \frac{v}{2{\mathcal{D}}}\displaystyle \frac{\partial {f}_{0}}{\partial z}.\end{eqnarray} \tag{ 2 }$$

Thus, the diffusion equation for ${f}_{0}$ takes the following form:

$$\begin{eqnarray}&&\displaystyle \frac{\partial {f}_{0}}{\partial t}=\displaystyle \frac{{v}^{2}}{4}\displaystyle \frac{\partial }{\partial z}\left\langle\displaystyle \frac{1-{\mu }^{2}}{{\mathcal{D}}}\right\rangle\displaystyle \frac{\partial {f}_{0}}{\partial z}.\end{eqnarray} \tag{ 3 }$$

A questionable point of course is neglecting $\partial f/\partial t$ in favor for $v\partial f/\partial z$ in Equation (2). It is somewhat justified by the small parameter ${{\mathcal{D}}}^{-1}\ll 1$ in the final result, given by Equation (3), hopefully making $\partial f/\partial t$ small. On the other hand, this is true for $\partial {f}_{0}/\partial t$, but not necessarily for $\partial f/\partial t$, since the latter may also contain the rapidly decaying anisotropic part, $\tilde{f}=f-{f}_{0}$, of the initial CR distribution. For ${\mathcal{D}}t\gtrsim 1$, however, $\tilde{f}$ must die out, and neglecting $\partial f/\partial t$ appears plausible for the long-term CR transport. At the same time, $\partial {f}_{0}/\partial t$ is large when ${f}_{0}$ is initially very narrow in z, such as in the fundamental solution. In the sequel, these aspects of the CR propagation will be key in choosing an appropriate asymptotic reduction method.

However convincing the justification, the CR diffusion model encounters the problem of a superluminal, or simply "too-fast" particle propagation. Although rather common for diffusive models, the problem is largely ignorable as long as the number of such particles remains small. There are cases, however, such as the propagation of ultra high-energy CR, where this problem must be addressed (Aloisio et al. 2009). Various attempts, starting as early as the 60s, e.g., (Axford 1965), have been made to devise a better transport equation for CRs. Unfortunately, in our view, they lack mathematical rigor and clarity and sometimes lead to inconsistent results.

In the most recent telegraph model, due to Litvinenko & Schlickeiser (2013), a higher order in $1/{\mathcal{D}}\ll 1$ term was included by retaining $\partial f/\partial t$, dropped in the simplest derivation above. This strategy gave rise to an additional ${\partial }^{2}{f}_{0}/\partial {t}^{2}$ term in the "master" equation. This additional term transforms Equation (3) into a "telegraph" equation:

$$\begin{eqnarray}&&\displaystyle \frac{\partial {f}_{0}}{\partial t}+T\displaystyle \frac{{\partial }^{2}{f}_{0}}{\partial {t}^{2}}=\displaystyle \frac{\partial }{\partial z}k\displaystyle \frac{\partial {f}_{0}}{\partial z}+\displaystyle \frac{k}{L}\displaystyle \frac{\partial {f}_{0}}{\partial z}\end{eqnarray} \tag{ 4 }$$

with

$$\begin{eqnarray}k & = & \displaystyle \frac{{v}^{2}}{4}\left\langle\displaystyle \frac{1-{\mu }^{2}}{{\mathcal{D}}}\right\rangle,\\ T & = & \left\langle{\left[{\displaystyle \int }_{0}^{\mu }d\mu /{\mathcal{D}}\right]}^{2}\right\rangle/\left\langle\displaystyle \frac{1-{\mu }^{2}}{{\mathcal{D}}}\right\rangle,\\ {L}^{-1} & = & -{B}^{-1}\partial B/\partial z.\end{eqnarray} \tag{ 5 }$$

For the sake of comparison with earlier telegraph equation results that will be made in Section 4, we have added here the magnetic focusing effect (the last term on the right-hand side (rhs) with B $(z)$ being the magnetic field), not initially included in Equation (1). Equation (4) is just a linear equation that can be solved immediately. The fundamental solution to Equation (4) that starts off from a $\delta (z)$ distribution, instantaneously released at t = 0, is as follows (e.g., Goldstein 1951; Axford 1965; Schwadron & Gombosi 1994; Litvinenko et al. 2015, ${L}^{-1}=0$, for simplicity)

$$\begin{eqnarray}{f}_{0} & = & \displaystyle \frac{1}{2}{e}^{-t/2T}\left\{\delta \left(z-\sqrt{\displaystyle \frac{k}{T}}t\right)+\delta \left(z+\sqrt{\displaystyle \frac{k}{T}}t\right)\right.\\ & & +\displaystyle \frac{H\left(\sqrt{k/T}t-| z| \right)}{2\sqrt{{kT}}}\left[{I}_{0}\left(\displaystyle \frac{1}{2}\sqrt{\displaystyle \frac{{t}^{2}}{{T}^{2}}-\displaystyle \frac{{z}^{2}}{{kT}}}\right)\right.\\ & & +\left.\left.\displaystyle \frac{t}{\sqrt{T\left({{kt}}^{2}-{{Tz}}^{2}\right)}}{I}_{1}\left(\displaystyle \frac{1}{2}\sqrt{\displaystyle \frac{{t}^{2}}{{T}^{2}}-\displaystyle \frac{{z}^{2}}{{kT}}}\right)\right]\right\}.\end{eqnarray} \tag{ 6 }$$

Here ${I}_{\mathrm{0,1}}$ denote the modified Bessel functions and the H denotes the Heaviside unit function.

One promising aspect of the telegraph equation is that it allows for a ballistic mode of CR propagation when the initial conditions empower the higher-order derivative terms to dominate (at least in the early phase of evolution). If, in addition, T has a proper value, the bulk speed of CRs may also be realistic. For example, this speed was derived in Earl (1973) to be $v/\sqrt{3}$, which has also been used earlier by Axford (1965). This is just the rms velocity projection of an isotropic, one-sided CR distribution on the z-axis, which appears to resolve the issue with the superluminal propagation. What is worrisome here, is that this bulk speed essentially requires a one-sided "isotropic" CR distribution, which, of course, is highly anisotropic overall, contrary to the basic assumption of most treatments. Thus, we need to start with isotropic initial distribution, but, taking it as narrow and symmetric in z (say Gaussian), the higher-order derivative terms will, again, dominate in Equation (4) and the single CR pulse will split into two, propagating in opposite directions at the speeds $\pm \sqrt{k/T}$ (as $\partial {f}_{0}/\partial t=0$ at t = 0, according to Equation (1)). Neglecting $\partial {f}_{0}/\partial t$ (which is justified for ${f}_{0}$  being sufficiently narrow in z), the solution is simply ${f}_{0}(z,t)=F\left(z-\sqrt{k/T}t\right)+F\left(z+\sqrt{k/T}t\right)$, where $2F(z)=f$ $(z,0)$. This result casts doubts on whether the telegraph term can be dominant under the assumption of frequent CR scattering (asymptotic expansion in small $1/{\mathcal{D}}$) because the initially sharp profile does not spread. Bringing $\partial {f}_{0}/\partial t$ back into the equation will only damp, but not spread, the profile, as is clearly seen from the solution given in Equation (6)), where $F\left(z\pm \sqrt{k/T}t\right)=\delta \left(z\pm \sqrt{k/T}t\right)$. Besides, the bulk CR speed $\sqrt{k/T}=v/\sqrt{3}$ for an isotropic scattering ${\mathcal{D}}=1$, though implicitly confirmed in the recent derivation of the telegraph equation for an arbitrary ${\mathcal{D}}(\mu )$ by Litvinenko & Schlickeiser (2013), is not universally accepted. Gombosi et al. (1993), Pauls et al. (1993), and Schwadron & Gombosi (1994), using simplified forms of $D(\mu )$, advocate the value $\sqrt{5/11}v$ for the propagation speed.

The last result is consistent with our calculations below, but with strong reservations regarding the telegraph equation set out later in the paper. Here we merely note that the solution of the telegraph equation specifically considered by Litvinenko & Schlickeiser (2013), which does not have the property of splitting the initially narrow pulse into two, does not conserve the total CR number $N=\int {f}_{0}{dz}$. It starts off from N = 0 which is unphysical. Because the equation has no source on its rhs to conserve N, two δ-pulses in Equation (6) are necessary. Those have been added to the treatment by Litvinenko & Noble (2013), but the δ-pulses have not been shown on their plots, for obvious reasons. Thus, the comparison of this solution with the solution of the original scattering problem is rather misleading. The disagreement of the propagation speed $\sqrt{k/T}$ is also critical because the solution in Equation (6) is cut off at a point z moving with this speed. For $z\lt \sqrt{k/T}t$, the profile is close to a Gaussian (for $t\gg T$, where T is the scattering time); thus, small variations in the speed can produce significant variations in the solution. We will also return to this later.

Another disadvantage of the telegraph Equation (4) is that it is no longer an evolution equation and requires the time derivative ${\partial }_{t}{f}_{0}$ as an initial condition. Although this can be inferred from the angular distribution at t = 0 using Equation (1), the "telegraph" description of CR transport is not self-contained. We show below that the T-term in Equation (4) is subdominant in an asymptotic series for ${\mathcal{D}}t\gtrsim 1$, thus representing transients in the CR transport. Strictly speaking, it should be omitted in the asymptotic transport description along with the small hyper-diffusion term $\sim {\partial }^{4}{f}_{0}/\partial {z}^{4}$, particularly if the term $\sim {\partial }^{3}{f}_{0}/\partial {z}^{3}$ does not vanish. The latter was not included in Equation (4) because it was obtained by applying an insufficient direct iteration to Equation (1) in (Litvinenko & Schlickeiser 2013), or assuming symmetric scattering, ${\mathcal{D}}(-\mu )={\mathcal{D}}(\mu )$, when it indeed vanishes. This symmetry restriction was relaxed in Pauls et al. (1993).

There are several reasons why we undertake the derivation of the master equation to a higher (fourth) order of approximation for an arbitrary ${\mathcal{D}}(\mu )$. First, it is necessary to clarify the role of the telegraph term entertained in the literature as an allegedly viable alternative to the standard diffusion model. Second, it is important to obtain the transport coefficients that are valid for arbitrary ${\mathcal{D}}(\mu )$, that is for an arbitrary spectrum of magnetic fluctuations. As we will show, the previous such derivation due to Litvinenko & Schlickeiser (2013), does not include the third order term, while including only one fourth-order term, while there are more such terms (see Equation (27)). Furthermore, the diffusion Equation (3) supplemented by a convective term $u(z)\partial {f}_{0}/\partial z$ for the case of the bulk fluid (scattering center) motion with velocity u, has long been and remains the main tool for diffusive shock acceleration (DSA) models. An accurate assessment of the next nonvanishing term, not included in Equation (3), is thus utterly important for the DSA, particularly because claims are being made about the necessity to include the telegraph term in the CR transport. In most DSA applications, it is crucial to allow not only for an arbitrary fluctuation spectrum ${\mathcal{D}}$ but for its dependence upon ${f}_{0}$ as well. This dependence directly affects the particle spectrum and acceleration time. We will discuss these aspects briefly in Section 6.

In the next section, the basic transport equation with magnetic focusing is introduced and the shortcomings of a reduction scheme based on direct iterations are demonstrated. The appropriate asymptotic method is elaborated in Section 3. Apart from what we already discussed regarding the telegraph equation, the objective of Section 3 is to create a framework that is also suitable for nonlinear (e.g., Ptuskin et al. 2008; Malkov et al. 2010b) and quasi-linear (Fujita et al. 2011; Malkov et al. 2013) versions of CR transport, which are important for both the DSA and for the subsequent escape of the accelerated CR. In these settings, the CR pressure is high enough to strongly modify, at least, the pitch-angle diffusion coefficient ${\mathcal{D}}$ and possibly the shock structure itself (Malkov et al. 2010b). In Sections 4 and 5, the implications of our results for the telegraph model and for the long-time CR propagation are discussed, while Section 6 concludes the paper.

Energetic particles (e.g., CRs) in a magnetic field, slowly varying on the particle gyro-scale, are transported according to the following gyro-phase averaged equation, (e.g., Vedenov et al. 1962; Jokipii 1966; Kulsrud 2005)

$$\begin{eqnarray}&&\displaystyle \frac{\partial f}{\partial t}+v\mu \displaystyle \frac{\partial f}{\partial z}+v\displaystyle \frac{\sigma }{2}\left(1-{\mu }^{2}\right)\displaystyle \frac{\partial f}{\partial \mu }=\displaystyle \frac{\partial }{\partial \mu }\nu D(\mu )\left(1-{\mu }^{2}\right)\displaystyle \frac{\partial f}{\partial \mu }.\end{eqnarray} \tag{ 7 }$$

Here v and μ are the particle velocity and pitch angle, z points in the local field direction, $\sigma =-{B}^{-1}\partial B/\partial z$ is the magnetic field inverse scale, and ν is the pitch angle scattering rate, while $D(\mu )\sim 1$ depends on the spectrum of magnetic fluctuations. As the fastest transport is assumed to be in μ, we introduce the following small parameter,

$$\begin{eqnarray}&&\varepsilon \equiv \displaystyle \frac{v}{l\nu }\equiv \displaystyle \frac{\lambda }{l}\ll 1,\end{eqnarray} \tag{ 8 }$$

where λ is the particle mean free path and l is a characteristic scale that should be chosen depending on the problem considered. One option is the scale of $B(z)$, in which case $l\sim {\sigma }^{-1}$. If the CR source is present on the rhs of Equation (7), its scale can be taken as l. Finally, l can be the scale of an initial CR distribution. Strictly speaking, the shortest of these scales should be taken as l. The problem with the initial CR distribution is that, in the most interesting case of the fundamental solution, this scale is zero. Therefore, over the initial period of CR spreading, before the actual CR scale $l(t)\sim f/(\partial f/\partial z)$ exceeds the m.f.p. λ, direct asymptotic expansions in small remain inaccurate. The goal here is to choose the least inaccurate out of all possible expansion schemes. At a minimum, it should be the one that does not introduce additional singularities, apart from the initial delta function $\delta (z)$, that must spread out under the particle recession and collisions. Therefore, while taking $l=\mathrm{const}$ in Equation (8), and assuming $l\gg \lambda$, caution will be exercised during the initial phase of the CR relaxation when the terms with higher spatial and time derivatives are large, even if they contain small factors ${\varepsilon }^{n}\ll 1$. By measuring time in ${\nu }^{-1}$, z in l, and simply replacing $\sigma l\to \sigma \sim 1$, the above equation transforms as follows.

$$\begin{eqnarray}&&\displaystyle \frac{\partial f}{\partial t}-\displaystyle \frac{\partial }{\partial \mu }D(\mu )\left(1-{\mu }^{2}\right)\displaystyle \frac{\partial f}{\partial \mu }\\ &&\quad =-\varepsilon \left(\mu \displaystyle \frac{\partial f}{\partial z}+\displaystyle \frac{\sigma }{2}\left(1-{\mu }^{2}\right)\displaystyle \frac{\partial f}{\partial \mu }\right).\end{eqnarray} \tag{ 9 }$$

A suitable scheme for asymptotic reduction of the above equation using $\varepsilon \ll 1$ is due to Chapman and Enskog, suggested in the development of earlier ideas by Hilbert (a good discussion of the history of this method with mathematical details is given by Cercignani 1988). Originally, it was applied to the Boltzmann equation in a strongly collisional regime. Similar approaches have been used in plasma physics, e.g., in regards to the hydrodynamic description of collisional magnetized plasmas (Braginskii 1965) and the problem of run-away electrons (Gurevich 1961; Kruskal & Bernstein 1964).

Regardless of the asymptotic scheme, Equation (9) suggests to seek f as a series in

$$\begin{eqnarray}&&f={f}_{0}+\varepsilon {f}_{1}+{\varepsilon }^{2}{f}_{2}+...\equiv {f}_{0}+\tilde{f}\end{eqnarray} \tag{ 10 }$$

where

$$\begin{eqnarray}&&\langle f\rangle ={f}_{0},\ \ \mathrm{with}\ \ \ \langle \cdot \rangle =\displaystyle \frac{1}{2}{\int }_{-1}^{1}(\cdot )d\mu ,\end{eqnarray} \tag{ 11 }$$

so that $\tilde{\langle f\rangle }=\left\langle{f}_{n\gt 0}\right\rangle=0$. The equation for ${f}_{0}$, which is the main ("master") equation of the method, takes the following form.

$$\begin{eqnarray}\displaystyle \frac{\partial {f}_{0}}{\partial t} & = & -\varepsilon \left(\displaystyle \frac{\partial }{\partial z}+\sigma \right)\langle \mu f\rangle \\ & = & -\displaystyle \frac{{\varepsilon }^{2}}{2}\left(\displaystyle \frac{\partial }{\partial z}+\sigma \right)\displaystyle \sum _{n=1}^{\infty }{\varepsilon }^{n-1}\left\langle\left(1-{\mu }^{2}\right)\displaystyle \frac{\partial {f}_{n}}{\partial \mu }\right\rangle.\end{eqnarray} \tag{ 12 }$$

We see from this equation that, similarly to the case of Lorentz's gas in an electric field (Gurevich 1961; Kruskal & Bernstein 1964), ${f}_{0}$ depends on the "slow time" ${t}_{2}={\varepsilon }^{2}t$ rather than on t. Indeed, the two problems are similar in that they describe diffusive expansion of particles in phase space. The expansion occurs in the z-direction for the CR diffusion problem and in energy for runaway electrons. The expansion is driven by a rapid isotropization in pitch angle plus the convection in the z-direction, or acceleration in the electric field direction, for the CR transport and electron runaway, respectively.

The slow dependence of ${f}_{0}$ on time in Equation (12) may suggest that one should attribute the time derivative term in Equation (9) to a higher order approximation (thus moving it to the rhs). Such ordering has been employed by Litvinenko & Schlickeiser (2013) and the term $\propto {\partial }^{2}{f}_{0}/\partial {t}^{2}$ has been produced in Equation (12). Obviously, a continuation of this process would result in progressively higher time derivatives of ${f}_{0}$, corresponding to shorter and shorter times in the initial relaxation. These transient phenomena will be removed using the Chapman–Enskog asymptotic reduction scheme in the next section.

Unlike ${f}_{0}$, $\tilde{f}$ in Equation (10) does depend on a "fast" time t. Therefore, it is illegitimate to attribute the first term on the left-hand side (lhs) of Equation (9) to any order of approximation different from that of the second term, notwithstanding its fast decay for $t\gtrsim 1$. Thus, using Equations (9)–(10), we must apply the following ordering

$$\begin{eqnarray}&&\displaystyle \frac{\partial {f}_{n}}{\partial t}-\displaystyle \frac{\partial }{\partial \mu }D(\mu )\left(1-{\mu }^{2}\right)\displaystyle \frac{\partial {f}_{n}}{\partial \mu }\\ &&\qquad =-\mu \displaystyle \frac{\partial {f}_{n-1}}{\partial z}-\displaystyle \frac{\sigma }{2}\left(1-{\mu }^{2}\right)\displaystyle \frac{\partial {f}_{n-1}}{\partial \mu }.\end{eqnarray} \tag{ 13 }$$

The above expansion scheme is sufficient to recover the leading order of ${f}_{0}$ evolution from Equation (12) by substituting $\partial {f}_{1}/\partial \mu \approx -{(2D)}^{-1}\partial {f}_{0}/\partial z$, obtained from the last equation for $t\gtrsim 1$. However, this scheme is not suitable for determining f_n for $n\geqslant 2$ to submit to Equation (12). Indeed, as it may be seen from Equation (13), the solubility condition for ${f}_{2}$ at $t\gg 1$ is $(\partial /\partial z+\sigma )\left\langle\left(1-{\mu }^{2}\right)\partial {f}_{1}/\partial \mu \right\rangle$ $\approx -(\partial /\partial z+\sigma )\left\langle\left(1-{\mu }^{2}\right)/2D\right\rangle\partial {f}_{0}/\partial z=0$. This is clearly too strong of a restriction. The reason for this inconsistency of the direct asymptotic expansion is that ${f}_{0}$ depends on time much slower than ${f}_{n\gt 0}$; thus, a slow time ${t}_{2}={\varepsilon }^{2}t$ needs to be taken into consideration. The Chapman–Enskog method has been developed for such cases, and we will make use of it in the next section.

As we have seen, the asymptotic reduction of the original CR propagation problem, given by Equation (9), to its isotropic part cannot proceed to higher orders of approximation using a simple asymptotic series in Equation (10) and requires a multi-time asymptotic expansion. In the Chapman–Enskog method, the operator $\partial /\partial t$ is expanded instead. Its purpose is to avoid unwanted higher time derivatives to appear in higher orders of approximation. This is very similar to, e.g., a secular growth in perturbed oscillations of dynamical systems. To eliminate the secular terms, one seeks to alter (also expand in a small parameter) the frequency of the zero-order motion, which is similar to the $\partial /\partial t$ expansion. One example of such an approach may be found in a derivation of hydrodynamic equations for strongly collisional, but magnetized plasmas, starting from the Boltzmann equation (Mikhailovsky 1967). The classical monograph by Chapman & Cowling (1991; Ch.VIII) gives another example of a subdivision of the $\partial /\partial t$ operator for solving the transport problem in a non-uniform gas-mixture. Expanding $\partial /\partial t$ operators eliminates secular terms, such as the telegraph term. Perhaps more customary today, and equivalently, is to introduce a hierarchy of formally independent time variables (e.g., Nayfeh 1981) $t\to {t}_{0},\ {t}_{1},...$, so that

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}=\displaystyle \frac{\partial }{\partial {t}_{0}}+\varepsilon \displaystyle \frac{\partial }{\partial {t}_{1}}+{\varepsilon }^{2}\displaystyle \frac{\partial }{\partial {t}_{2}}...\end{eqnarray} \tag{ 14 }$$

Instead of Equation (13), from Equation (9), we have

$$\begin{eqnarray}&&\displaystyle \frac{\partial {f}_{n}}{\partial {t}_{0}}-\displaystyle \frac{\partial }{\partial \mu }D(\mu )\left(1-{\mu }^{2}\right)\displaystyle \frac{\partial {f}_{n}}{\partial \mu }\\ &&\qquad =-\mu \displaystyle \frac{\partial {f}_{n-1}}{\partial z}-\displaystyle \frac{\sigma }{2}\left(1-{\mu }^{2}\right)\displaystyle \frac{\partial {f}_{n-1}}{\partial \mu }-\displaystyle \sum _{k=1}^{n}\displaystyle \frac{\partial {f}_{n-k}}{\partial {t}_{k}}\\ &&\qquad \equiv {{\mathcal{L}}}_{n-1}[f]\left({t}_{0},...,{t}_{n};\mu ,z),\right.\end{eqnarray} \tag{ 15 }$$

where the conditions ${f}_{n\lt 0}=0$ are implied. The solution of this equation should be sought in the following form

$$\begin{eqnarray}&&{f}_{n}={\bar{f}}_{n}\left({t}_{2},{t}_{3},...;\mu \right)+{\tilde{f}}_{n}\left({t}_{0},{t}_{1},...;\mu ),\right.\end{eqnarray} \tag{ 16 }$$

where ${\tilde{f}}_{n}$ and ${\bar{f}}_{n}$ are chosen to satisfy, respectively, the following two equations.

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\tilde{f}}_{n}}{\partial {t}_{0}}-\displaystyle \frac{\partial }{\partial \mu }D(\mu )\left(1-{\mu }^{2}\right)\displaystyle \frac{\partial {\tilde{f}}_{n}}{\partial \mu }={{\mathcal{L}}}_{n-1}\left[\tilde{f}\right]\left({t}_{0},...,{t}_{n};\mu ,z\right)\end{eqnarray} \tag{ 17 }$$

and

$$\begin{eqnarray}&&-\displaystyle \frac{\partial }{\partial \mu }D(\mu )\left(1-{\mu }^{2}\right)\displaystyle \frac{\partial {\bar{f}}_{n}}{\partial \mu }={{\mathcal{L}}}_{n-1}\left[\bar{f}\right]\left({t}_{2},...,{t}_{n};\mu ,z)\mathrm{.}\right.\end{eqnarray} \tag{ 18 }$$

The solution for ${\tilde{f}}_{n}$ is as follows

$$\begin{eqnarray}&&{\tilde{f}}_{n}=\displaystyle \sum _{k=1}^{\infty }{C}_{k}^{(n)}\left({t}_{0}\right){e}^{-{\lambda }_{k}{t}_{0}}{\psi }_{k}(\mu )\end{eqnarray} \tag{ 19 }$$

and it can be evaluated for arbitrary n by expanding both sides of Equation (17) in a series of eigenfunctions of the diffusion operator on its lhs:

$$\begin{eqnarray*}&&-\displaystyle \frac{\partial }{\partial \mu }D(\mu )\left(1-{\mu }^{2}\right)\displaystyle \frac{\partial {\psi }_{k}}{\partial \mu }={\lambda }_{k}{\psi }_{k}.\end{eqnarray*}$$

For D = 1, for example, ${\psi }_{k}$ are the Legendre polynomials with ${\lambda }_{k}=k(k+1)$, $k=0,1,...$. The time dependent coefficients ${C}_{k}^{(n)}$ are determined by the initial values of ${\tilde{f}}_{n}$ (the anisotropic part of the initial CR distribution) and the rhs of Equation (17), that depends on ${\tilde{f}}_{n-1}$, obtained at the preceding step. It is seen, however, that ${\tilde{f}}_{n}$ exponentially decay in time for $t\gtrsim 1$ and we may ignore them³ because we are primarily interested in evolving the system over time $t\gtrsim {\varepsilon }^{-2}\gg 1$ and even longer. Starting from n = 0 and using Equation (15), for the slowly varying part of f, we have

$$\begin{eqnarray}&&\displaystyle \frac{\partial {f}_{0}}{\partial {t}_{0}}=0.\end{eqnarray} \tag{ 20 }$$

The solubility condition for ${f}_{1}$ (obtained by integrating both sides of Equation (15) in μ) also gives a trivial result

$$\begin{eqnarray}&&\displaystyle \frac{\partial {f}_{0}}{\partial {t}_{1}}=0,\end{eqnarray} \tag{ 21 }$$

so the last two conditions are consistent with the suggested decomposition in Equation (16), since from Equation (18) with n = 1, we have

$$\begin{eqnarray}&&{\bar{f}}_{1}=-\displaystyle \frac{1}{2}W\displaystyle \frac{\partial {f}_{0}}{\partial z}\end{eqnarray} \tag{ 22 }$$

and, thus both ${\bar{f}}_{0}$ and ${\bar{f}}_{1}$ are, indeed, independent of ${t}_{0}$ and ${t}_{1}$. We have introduced the function $W(\mu )$ here by the following two relations

$$\begin{eqnarray}&&\displaystyle \frac{\partial W}{\partial \mu }=\displaystyle \frac{1}{D},\ \ \ \langle W\rangle =0.\end{eqnarray} \tag{ 23 }$$

The solubility condition for ${f}_{2}$ yields the nontrivial and well-known (e.g., Jokipii 1966) result, which is actually the leading term of the $\partial {f}_{0}/\partial t$ expansion in $\varepsilon \ll 1$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {f}_{0}}{\partial {t}_{2}}=\displaystyle \frac{1}{4}\left(\displaystyle \frac{\partial }{\partial z}+\sigma \right)\kappa \displaystyle \frac{\partial {f}_{0}}{\partial z},\end{eqnarray} \tag{ 24 }$$

where

$$\begin{eqnarray*}&&\kappa =\left\langle\displaystyle \frac{\left(1-{\mu }^{2}\right)}{D}\right\rangle.\end{eqnarray*}$$

The solubility conditions for ${f}_{3}$, ${f}_{4},\ldots$ will generate the higher order terms of our expansion, which, after some algebra, can be manipulated into the following expressions for the third and fourth orders of approximation

$$\begin{eqnarray}&&\displaystyle \frac{\partial {f}_{0}}{\partial {t}_{3}}=-\displaystyle \frac{1}{4}\left(\displaystyle \frac{\partial }{\partial z}+\sigma \right)\left(\displaystyle \frac{\partial }{\partial z}+\displaystyle \frac{\sigma }{2}\right)\left\langle\mu {W}^{2}\right\rangle\displaystyle \frac{\partial {f}_{0}}{\partial z}\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}\displaystyle \frac{\partial {f}_{0}}{\partial {t}_{4}} & = & \displaystyle \frac{1}{8}\left(\displaystyle \frac{\partial }{\partial z}+\sigma \right)\left\{{\left(\displaystyle \frac{\partial }{\partial z}+\displaystyle \frac{\sigma }{2}\right)}^{2}\left\langle{W}^{2}\left({U}^{\prime }-\kappa \right)\right\rangle\right.\\ & & +\left.\displaystyle \frac{1}{2}\left(\displaystyle \frac{\partial }{\partial z}+\sigma \right)\displaystyle \frac{\partial }{\partial z}\left\langle\displaystyle \frac{{[\kappa (1-\mu )+U]}^{2}}{D\left(1-{\mu }^{2}\right)}\right\rangle\right\}\displaystyle \frac{\partial {f}_{0}}{\partial z}.\end{eqnarray} \tag{ 26 }$$

We have denoted

$$\begin{eqnarray*}&&U\equiv {\int }_{1}^{\mu }\displaystyle \frac{1-{\mu }^{2}}{D}d\mu ,\end{eqnarray*}$$

and ${U}^{\prime }=\partial U/\partial \mu$. The pitch-angle diffusion coefficient $D(\mu )$ and magnetic focusing σ are considered z-independent for simplicity, a limitation that can be easily relaxed by rearranging the operators containing $\partial /\partial z$ in Equation (26). We can proceed to higher orders of approximation ad infinitum since terms containing $\left\langle\left(1-{\mu }^{2}\right)\partial {f}_{n}/\partial \mu \right\rangle$ can be expressed through ${f}_{n-1},\ {f}_{n-2},\ldots$. According to Equations (20)–(21), of interest is the evolution of ${f}_{0}$ on the timescales ${t}_{2}\gtrsim 1$ or $t\gtrsim {\varepsilon }^{-2}$; thus, as we already mentioned, the contributions of ${\tilde{f}}_{n}(\mu )$ to all of the solubility conditions, similar to those given by Equations (24)–(26), have to be dropped (because they become exponentially small) and only ${\bar{f}}_{n}(\mu )$-contributions should be retained. Using Equations (20)–(21) and (24)–(26) to form the combinations ${\varepsilon }^{n}{\partial }^{n}{f}_{0}/\partial {t}_{n}$ and summing up both sides, on the lhs of the resulting equation, we simply obtain $\partial {f}_{0}/\partial t$ (see Equation (14)). Therefore, the evolution of ${f}_{0}$ up to the fourth order in takes the following form.

$$\begin{eqnarray}\displaystyle \frac{\partial {f}_{0}}{\partial t} & = & \displaystyle \frac{{\epsilon }^{2}}{4}{\partial }_{z}^{\prime }\left\{\kappa -\epsilon {\partial }_{z}^{\prime\prime }\langle \mu {W}^{2}\rangle -\displaystyle \frac{{\epsilon }^{2}}{2}\left[{K}_{1}{\left({\partial }_{z}^{\prime\prime }\right)}^{2}\right.\right.\\ & & -\left.\left.{K}_{2}{\partial }_{z}^{\prime }{\partial }_{z}\right]\right\}\displaystyle \frac{\partial {f}_{0}}{\partial z},\end{eqnarray} \tag{ 27 }$$

where ${\partial }_{z}^{\prime }={\partial }_{z}+\sigma$, ${\partial }_{z}^{\prime\prime }={\partial }_{z}+\sigma /2$, and

$$\begin{eqnarray}&&{K}_{1}=\left\langle{W}^{2}\left(\kappa -{U}^{\prime }\right)\right\rangle,\ \ {K}_{2}=\displaystyle \frac{1}{2}\left\langle\displaystyle \frac{{[\kappa (1-\mu )+U]}^{2}}{D\left(1-{\mu }^{2}\right)}\right\rangle.\end{eqnarray} \tag{ 28 }$$

The above algorithm allows one to obtain the master equation to arbitrary order in $\varepsilon$. By construction, in no order of approximation will higher time derivatives emerge, as has been devised by Chapman and Enskog. We have truncated this process at the fourth order, ${\varepsilon }^{4}$. As we show in the next section, this is the lowest order required to relate the above result to the telegraph equation. It also gives the first non-vanishing correction to the standard CR diffusion model in an important case $\left\langle\mu {W}^{2}\right\rangle=0$, which is fulfilled, in particular, for $D(-\mu )=D(\mu )$. Higher order terms can be calculated at the expense of a more involved algebra, but we argue below that such calculations would not change the results significantly.

In contrast to the telegraph equation given by Equations (4) and (5), that has been derived by Litvinenko & Schlickeiser (2013), using a direct iteration of Equation (7) with no explicit ordering of the emerging terms, Equation (27) is derived to the ${\varepsilon }^{4}$ order of approximation with an ${\varepsilon }^{n}$ factor labeling each term. Yet, it has no second order time derivative, which is inconsistent with Equation (4). Below, we demonstrate that Equation (27) can still be converted to the telegraph form, however, with additional terms absent from Equations (4)–(5). Although Equation (27) is obtained by a broadly applicable Chapman–Enskog method, its reduction to the telegraph form below is more restrictive and should be taken with a grain of salt for the reasons we discuss later.

Several versions of the telegraph equation have been obtained using different methods but, unfortunately, many of them do not offer clear ordering, as, e.g., Equation (4) derived by Litvinenko & Schlickeiser (2013). In an earlier treatment by Earl (1973), an eigenfunction expansion was truncated with no transparent assessment of discarded terms. As we mentioned already, many treatments do not systematically eliminate short timescales, which are irrelevant to the long-time evolution of the isotropic part of the CR distribution. In principle, this is acceptable if the reduction scheme is based on an exact solution of the original equation, to include all required orders of approximation into the master equation. Such an approach, along with a nearly exhaustive analysis of the previous work has been presented in Schwadron & Gombosi (1994). Their treatment, however, is limited, by necessity, to a relatively simple $D(\mu )$ (i.e., power-law in μ).

To clarify the role of the higher order terms in Equation (27), we note that the rhs of this equation only represents the first three non-vanishing contributions from an infinite asymptotic series (in $\varepsilon \ll 1$) which we would obtain by continuing the reduction process described in the preceding section. This series may or may not converge to some linear (integral) operator in z. From a practical standpoint, the maximum order term that needs to be retained is either the first non-vanishing term or else it introduces a new property to the solution, such as symmetry breaking. Precisely the last aspect has been highlighted by Litvinenko & Schlickeiser (2013) who used the telegraph term to calculate the skewness of a CR pulse. From this angle, we examine the third and the fourth order below separately.

Third order equation. To this order, Equation (27) rewrites

$$\begin{eqnarray}&&\displaystyle \frac{\partial {f}_{0}}{\partial \tau }=V\displaystyle \frac{\partial {f}_{0}}{\partial z}+{\kappa }_{1}\displaystyle \frac{{\partial }^{2}{f}_{0}}{\partial {z}^{2}}-\varepsilon \left\langle\mu {W}^{2}\right\rangle\displaystyle \frac{{\partial }^{3}{f}_{0}}{\partial {z}^{3}}.\end{eqnarray} \tag{ 29 }$$

We have introduced the slow time $\tau ={\varepsilon }^{2}t/4$, as a natural timescale for the reduced system, and the following notation

$$\begin{eqnarray*}&&V=\sigma \left(\kappa -\displaystyle \frac{1}{2}\varepsilon \sigma \left\langle\mu {W}^{2}\right\rangle\right),\ \ \ {\kappa }_{1}=\kappa -\displaystyle \frac{3}{2}\varepsilon \sigma \left\langle\mu {W}^{2}\right\rangle.\end{eqnarray*}$$

Note, that using τ instead of t makes the terms lower by two orders in , but describing them in the sequel we will use their original order as it stands in Equation (27), rather than Equations (30) or (32).

Equation (29) can be solved using a Fourier transform and integral representations of Airy functions. We consider some basic properties of this solution using the moments of ${f}_{0}$. In particular, it is seen from this equation that the skewness of a CR pulse, propagating at the bulk speed V, arises in this order of approximation. Indeed, upon a Galilean transform to the reference frame moving with the speed $-V$, $z\to {z}^{\prime }=z+V\tau$, the above equation rewrites

$$\begin{eqnarray}&&\displaystyle \frac{\partial {f}_{0}}{\partial \tau }={\kappa }_{1}\displaystyle \frac{{\partial }^{2}{f}_{0}}{\partial {z}^{\prime 2}}-\varepsilon \left\langle\mu {W}^{2}\right\rangle\displaystyle \frac{{\partial }^{3}{f}_{0}}{\partial {z}^{\prime 3}}\end{eqnarray} \tag{ 30 }$$

so that the last term generates an antisymmetric component of ${f}_{0}\left({z}^{\prime }\right)$, even if ${f}_{0}\left({z}^{\prime }\right)$ is an even function of z' initially. By normalizing ${f}_{0}$ to unity,

$$\begin{eqnarray*}&&{\bar{f}}_{0}={\int }_{-\infty }^{\infty }{f}_{0}{{dz}}^{\prime }=1,\end{eqnarray*}$$

and assuming the coefficients in Equation (30) to be constant, for the moments of ${f}_{0}\left({z}^{\prime },t\right)$

$$\begin{eqnarray*}&&\bar{{z}^{\prime n}}={\int }_{-\infty }^{\infty }{f}_{0}{z}^{\prime n}{{dz}}^{\prime },\end{eqnarray*}$$

we obtain

$$\begin{eqnarray*}&&\displaystyle \frac{d}{d\tau }\bar{{z}^{\prime }}=0,\ \ \ \displaystyle \frac{d}{d\tau }\bar{{z}^{\prime 2}}=2{\kappa }_{1},\ \ \ \displaystyle \frac{d}{d\tau }\bar{{z}^{\prime 3}}=6{\kappa }_{1}\bar{{z}^{\prime }}+6\varepsilon \left\langle\mu {W}^{2}\right\rangle.\end{eqnarray*}$$

With no loss of generality, we may set $\bar{{z}^{\prime }}=0$ and, in addition, $\bar{{z}^{\prime 3}}=0$ at $\tau =0$, so that the skewness of the CR distribution changes in time as follows:

$$\begin{eqnarray*}&&S\equiv \displaystyle \frac{\bar{{z}^{\prime 3}}}{{\left(\bar{{z}^{\prime 2}}\right)}^{3/2}}=\displaystyle \frac{6\varepsilon \left\langle\mu {W}^{2}\right\rangle\tau }{{\left(\bar{{z}_{0}^{\prime 2}}+2{\kappa }_{1}\tau \right)}^{3/2}},\end{eqnarray*}$$

where "0" at $\bar{{z}^{\prime 2}}$ refers to its value at $\tau =0$. The skewness remains small and its maximum

$$\begin{eqnarray*}&&{S}_{\mathrm{max}}=\displaystyle \frac{2\varepsilon \left\langle\mu {W}^{2}\right\rangle}{{\kappa }_{1}\sqrt{3\bar{{z}_{0}^{\prime 2}}}}\end{eqnarray*}$$

is achieved at $\tau =\bar{{z}^{\prime 2}}/{\kappa }_{1}$. Not surprisingly, the skewness increases with decreasing initial spatial dispersion of CRs. Indeed, according to Equation (9), a narrow $f(z)$ generates strong pitch-angle anisotropy which, in combination with asymmetric pitch-angle scattering ($\left\langle\mu {W}^{2}\right\rangle\ne 0$), generates the spatial skewness of the CR pulse. For a τ that is not too small, the explicit form of the solution of Equation (30) may be easily written down by using, e.g., a Fourier transform in z' and the steepest descent estimate of its inversion

$$\begin{eqnarray}&&f\left({z}^{\prime },\tau \right)\simeq \displaystyle \frac{C}{\sqrt{\tau }}\mathrm{exp}\left[-\displaystyle \frac{{z}^{\prime 2}}{4{\kappa }_{1}\tau }\left(1-\displaystyle \frac{\varepsilon {z}^{\prime }}{2{\kappa }_{1}^{2}\tau }\left\langle\mu {W}^{2}\right\rangle\right)\right].\end{eqnarray} \tag{ 31 }$$

It is quite possible, however, that even this small effect does not occur because of the pitch angle scattering symmetry, that is $\left\langle\mu {W}^{2}\right\rangle=0$. In this case the solution remains diffusive and, to obtain corrections to it and to see where the telegraph term might come from, the next approximation needs to be considered.

Fourth order equation: the telegraph term. In the absence of ${\varepsilon }^{3}$ terms, that is, when $\left\langle\mu {W}^{2}\right\rangle=0$, Equation (27) takes the following form

$$\begin{eqnarray}&&\displaystyle \frac{\partial {f}_{0}}{\partial \tau }+{\varepsilon }^{2}\left({K}_{1}-{K}_{2}\right)\left(\sigma +\displaystyle \frac{1}{2}\displaystyle \frac{\partial }{\partial z}\right)\displaystyle \frac{{\partial }^{3}{f}_{0}}{\partial {z}^{3}}\\ &&\qquad ={\kappa }_{2}\displaystyle \frac{{\partial }^{2}{f}_{0}}{\partial {z}^{2}}+{V}_{2}\displaystyle \frac{\partial {f}_{0}}{\partial z}+{\mathcal{O}}\left({\varepsilon }^{4}\right)\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray*}&&{\kappa }_{2}=\kappa -\displaystyle \frac{{\sigma }^{2}{\varepsilon }^{2}}{2}\left(\displaystyle \frac{5}{4}{K}_{1}-{K}_{2}\right),\ \ \ {V}_{2}=\sigma \left[\kappa -\displaystyle \frac{{\sigma }^{2}{\varepsilon }^{2}}{8}{K}_{1}\right].\end{eqnarray*}$$

This equation, obtained with the Chapman–Enskog method, and the telegraph Equation (4), obtained by a direct iteration method, differ from each other in the second term on the lhs. The remaining terms of the two equations are equivalent, though not identical due to the insignificant $\sim {\varepsilon }^{4}$ corrections included in the coefficients ${\kappa }_{2}$ and V on the rhs of Equation (32).

To understand how the conflicting terms on the lhs of both equations are related, we note that within the regular ordering scheme leading to Equation (32) the term in question must remain small, being nominally an ${\varepsilon }^{4}$-term. However, as the conflicting terms in both equations are the higher-order derivatives, they may stick out from their order of approximation if the solution strongly varies in space and time. In order to preserve the integrity of the overall solution  in such events a multi-(time)scale or matched asymptotic expansion method is normally applied. We will argue that the telegraph equation approach to the CR transport does not handle this situation properly, as opposed to the Chapman–Enskog approach. However, this is not to say that the two terms cannot be mapped to each other, when they are well in the validity range in the above sense. Note that in a number of other treatments the telegraph term was tacitly handled as one of the dominant terms. We pointed out that indeed, the term in question on the lhs of Equation (32) is small only insofar as the z-derivative does not change its order of approximation due to strong inhomogeneity, whose scale should not be less than the m.f.p., λ. Thus, assuming strong inequality ${\partial }_{z}\ll {\lambda }^{-1}$ (or simply ${\partial }_{z}\sim 1$ in our dimensionless variables), we can express the high order spatial derivatives using a zero order ($\varepsilon \to 0$) version of this equation with sufficient accuracy. The result reads

$$\begin{eqnarray}\displaystyle \frac{\partial {f}_{0}}{\partial \tau }+\displaystyle \frac{{\varepsilon }^{2}}{2{\kappa }^{2}}\left({K}_{1}-{K}_{2}\right)\displaystyle \frac{{\partial }^{2}{f}_{0}}{\partial {\tau }^{2}} & = & \left(\kappa -\displaystyle \frac{{\sigma }^{2}{\varepsilon }^{2}}{8}{K}_{1}\right)\displaystyle \frac{{\partial }^{2}{f}_{0}}{\partial {z}^{2}}\\ & & +{V}_{2}\displaystyle \frac{\partial {f}_{0}}{\partial z}+{\mathcal{O}}\left({\varepsilon }^{4}\right).\end{eqnarray} \tag{ 33 }$$

This equation is indeed equivalent to the telegraph equation by its form, but the equivalence requires not only $\varepsilon \ll 1$, but also smooth variation in z and τ, as not to raise the actual value of these terms significantly. Under these equivalence conditions, both the telegraph and the hyperdiffusive transport terms are just the corrections and may be safely ignored (especially if the ${\varepsilon }^{3}$ contribution is not empty). On the contrary, when the higher derivatives strongly enhance these terms, the equations cannot be mapped to each other and their solutions are disparate. One of them (or even both) may become less accurate than the underlying leading order (diffusive) approximation. This is a quite common situation in asymptotic expansions when the form of the next order term should be selected on the ground of the least possible singularity it introduces into the expansion (compare to small denominator, secular growth, etc., in mechanical problems, where the higher order approximations, if handled blithely, only aggravate disagreement with true solutions). The telegraph term correction to the diffusive approximation appears to come from that variety because it generates singular components (δ- and Heaviside functions) that are not only inconsistent with the strong pitch angle scattering and resulting spatial diffusion, but with the scatter-free limit of the parent differential equation itself. Therefore, the singular part of the telegraph solution is inherited from the derivation of the telegraph equation. We compare the telegraph and hyperdiffusion type corrections to the basic diffusive propagation further in the next section.

We start with a relatively minor aspect of the differences between the two models. As we stated in Section 1, the telegraph coefficient in Equation (33) is inconsistent with some of the earlier derivations. In the simplest case $D=1,$ for example, after proper rescaling of τ and z, it turns out to be smaller than the term T in Equation (4) by a factor of $11/15$. On the other hand, this is consistent with the respective result obtained by Gombosi et al. (1993), Pauls et al. (1993), and Schwadron & Gombosi (1994). While the above difference may be considered rather quantitative, in the general case of $D(-\mu )\ne D(\mu )$, the appropriate equation for describing CR transport is that given by the lower, ${\varepsilon }^{3}$-order, not included in Equation (4).

More importantly, the telegraph version of Equation (32) given by Equation (33) is valid only if the telegraph term ($\propto {\varepsilon }^{2}$, a fourth order term) remains small compared to the other terms and the original ordering in Equation (32) is not violated by strong variations of the solution in space and time, as we pointed out earlier. We signify this by the "slow" time $\tau \sim {\varepsilon }^{2}t$. In most other treatments, t is used instead, which formally makes the telegraph term in Equation (33) appear as a zero order term. It is not important, of course, whether the term is labeled by ${\varepsilon }^{2}$ or not; what is important is that it is treated as a subordinate term. Attempts to make it dominant a posteriori violates assumptions that are essential for its derivation. This point is demonstrated below by repeating a simple calculation of the CR pulse skewness that we already made earlier working to the ${\varepsilon }^{3}$ order.

Litvinenko & Schlickeiser (2013) suggested to study an asymmetry (skewness) of a CR pulse propagating along the field under the action of magnetic focusing using the telegraph equation. This and some other characteristics of the CR pulse, such as the kurtosis, can be easily analyzed using the primary Equation (32). The calculation of the pulse skewness essentially repeats the one already done at the ${\varepsilon }^{3}$ level, where it is generated by asymmetric scattering, $D(-\mu )\ne D(\mu )$. Thus, transforming Equation (32) to the reference frame moving with the speed $-{V}_{2}$, that is $z\to {z}^{\prime }=z+{V}_{2}\tau$, we obtain

$$\begin{eqnarray}&&\displaystyle \frac{d}{d\tau }\bar{{z}^{\prime }}=0,\ \ \ \displaystyle \frac{d}{d\tau }\bar{{z}^{\prime 2}}=2{\kappa }_{2},\ \ \ \displaystyle \frac{d}{d\tau }\bar{{z}^{\prime 3}}=6{\varepsilon }^{2}\sigma \left({K}_{1}-{K}_{2}\right),\end{eqnarray} \tag{ 34 }$$

where we have, again, assumed $\bar{{z}^{\prime }}=0$. The skewness thus evolves in time as follows

$$\begin{eqnarray*}&&S=\displaystyle \frac{6{\varepsilon }^{2}\sigma \left({K}_{1}-{K}_{2}\right)\tau }{{\left(\bar{{z}_{0}^{\prime 2}}+2{\kappa }_{2}\tau \right)}^{3/2}}.\end{eqnarray*}$$

Unless $S(\tau )$ reaches its maximum very early it is fairly small. Because ${\tau }_{\mathrm{max}}=\bar{{z}_{0}^{\prime 2}}/{\kappa }_{2}$, the maximum value

$$\begin{eqnarray}&&{S}_{\mathrm{max}}=S\left({\tau }_{\mathrm{max}}\right)=\displaystyle \frac{2{\varepsilon }^{2}\sigma \left({K}_{1}-{K}_{2}\right)}{{\kappa }_{2}\sqrt{3\bar{{z}_{0}^{\prime 2}}}},\end{eqnarray} \tag{ 35 }$$

so that an initially symmetric CR pulse develops significant asymmetry only if $\bar{{z}_{0}^{\prime 2}}\lesssim {\varepsilon }^{4}$. This, however, would require ${\tau }_{\mathrm{max}}\sim {\varepsilon }^{4}$ (or, equivalently, ${t}_{\mathrm{max}}\sim {\varepsilon }^{2}$), in strong violation of the requirement $t\gtrsim 1$, established in Section 3. Note that significant pulse asymmetry obtained by Litvinenko & Schlickeiser (2013) using the telegraph equation was based on the fundamental solution to this equation, that is ${z}_{0}^{\prime }\to 0$. We argued in Section 1 that the early propagation phase is not adequately described by the telegraph equation so that the pulse asymmetry might have been overestimated in the above paper. We specify the validity range of the telegraph equation below.

Starting from a long time regime $\tau \gt \varepsilon {z}^{\prime }$, similarly to the ${\varepsilon }^{3}$ result given in Equation (31), from Equation (32), we find

$$\begin{eqnarray}f\left({z}^{\prime },\tau \right) & \simeq & \displaystyle \frac{C}{\sqrt{\tau }}\mathrm{exp}\left\{-\displaystyle \frac{{z}^{\prime 2}}{4{\kappa }_{2}\tau }\left[1-{\varepsilon }^{2}\displaystyle \frac{\left({K}_{1}-{K}_{2}\right){z}^{\prime }}{2{\kappa }_{2}\tau }\right.\right.\\ & & \times \left.\left.\left(\sigma -\displaystyle \frac{{z}^{\prime }}{4{\kappa }_{2}\tau }\right)\right]\right\}.\end{eqnarray} \tag{ 36 }$$

This propagation regime is not much different from the regular diffusion (the first term in the square bracket), so both the hyperdiffusion and telegraph models produce similar results because they are largely equivalent in this regime. It is worthwhile to write the requirement for the agreement between the two models in physical units, which is simply

$$\begin{eqnarray*}&&{vt}\gt {z}^{\prime }.\end{eqnarray*}$$

The point $z={vt}$ is close to the cut-off in the telegraph solution, $z/t=\sqrt{k/T}$ (for $\sigma =0$), Equation (6). It follows then that unless ${z}^{2}\gg {kt}$ ($t\gg T$), the cut-off strongly changes the overall solution.

The opposite case $\tau \lt \varepsilon {z}^{\prime }$, which corresponds to the initial phase of pulse relaxation, is the key to understanding the difference between the Chapman–Enskog and the telegraph methods. In this regime, they are not equivalent because the ${\varepsilon }^{4}$ term in Equation (32) cannot be neglected to make the transition to the telegraph Equation (33). Indeed, when propagation starts with an infinitely narrow pulse, in the early phase of its relaxation the higher z-derivatives are still too large for such transformation. A spatially narrow pulse automatically generates strong pitch angle anisotropy which, in turn, results in rapid time variation, making the telegraph term also large. It is this regime where both methods become questionable and, in addition, their predictions deviate from each other both quantitatively and qualitatively. We need to check first whether they are relevant to this regime.

The phase $\tau \lt \varepsilon {z}^{\prime }$ is well described by the Chapman–Enskog approach down to $\tau \sim 1/{\varepsilon }^{2}$ (${t}_{0}\sim 1$, Section 3). The situation with the telegraph equation is more complex, as the conversion from the Chapman–Enskog expansion is invalid, while independent derivations rarely provide clear ordering. A rigorous derivation in Schwadron & Gombosi (1994) requires the same assumptions that we made when transforming the hyperdiffusive equation into the telegraph equation, that is ${\partial }_{\tau }\sim {\partial }_{z}^{2}$ (${\epsilon }_{\tau }\sim {\epsilon }_{\lambda }^{2}$ under their nomenclature). So, the telegraph equation appears to be a subset of the hyperdiffusion equation valid only under the above ordering. It is likely to break down in the $\tau \lt \varepsilon {z}^{\prime }$ regime, in other words, near its cut-off. We support this premise by the following considerations.

Recently, Effenberger & Litvinenko (2014) and Litvinenko et al. (2015) have carried out simulations of the full scattering problem, corresponding to Equation (9). The results deviate from the telegraph solution precisely at the early phase of the pulse propagation, when the hyperdiffusion and telegraph models disagree. The two δ-function pulses with sharp fronts in its solution given by Equation (6) are not seen in the simulations. This is understandable, as such features are inconsistent with the underlying scattering problem. They should have been smeared out by scattering earlier, since the spatial profile is shown at five collision times (Figure 2 in Effenberger & Litvinenko 2014). Moreover, the δ-function pulses $\delta \left(z\pm \sqrt{k/T}t\right)$ that are an integral part of the telegraph solution, as they maintain its normalization, are irrelevant to the primary Equation (1), even without collisions. Indeed, if ${\mathcal{D}}=0$, and the initial condition is $\delta (z)$ being constant in μ for $-1\lt \mu \lt 1$, the scatter-free solution is $\delta (z-\mu {vt})$. Hence, ${f}_{0}(z,t)\equiv \int {fd}\mu /2={(2{vt})}^{-1}H({vt}-| z| )$. Therefore, a certain property of the telegraph equation allows the $\delta \left(z\pm \sqrt{k/T}t\right)$ and sharp front components (although with modified speed, Equation (6)) to survive multiple collisions. As these components are inconsistent with the underlying scattering problem, this property of the equation must have been acquired during its derivation. Obviously, it is rooted in the hyperbolic (telegraph) operator ${\partial }_{t}^{2}-(k/T){\partial }_{z}^{2}$ that allows the singular profiles to propagate without spreading.

By contrast, the hyperdiffusion equation does not require singular components, but, on the contrary, smears them out. We present an approximate solution of Equation (32) after neglecting magnetic focusing in ${\varepsilon }^{4}$-order terms in Equation (32). We also neglect the regular diffusion compared to the hyperdiffusion, which is acceptable during an early phase of pulse relaxation, $\tau \lt \varepsilon {z}^{\prime }$. Finally, we assume ${z}^{\prime }\gt 0$, as the solution is an even function of z'. The asymptotic result is as follows:

$$\begin{eqnarray}f & = & \displaystyle \frac{2}{3}\sqrt{\displaystyle \frac{2}{\pi }}{(4h\tau )}^{-1/6}{z}^{\prime -1/3}\mathrm{exp}\left(-\displaystyle \frac{3}{8}{4}^{-1/3}\displaystyle \frac{{z}^{\prime 4/3}}{{h}^{1/3}{\tau }^{1/3}}\right)\\ & & \times \mathrm{cos}\left(\displaystyle \frac{{3}^{3/2}}{8}{4}^{-1/3}\displaystyle \frac{{z}^{\prime 4/3}}{{h}^{1/3}{\tau }^{1/3}}\right).\end{eqnarray} \tag{ 37 }$$

We have denoted the hyperdiffusion constant $h={\varepsilon }^{2}\left({K}_{1}-{K}_{2}\right)/2$. More about this result and further discussion of the two conflicting approaches can be found in the appendix. We see that there is a considerable slow down of the CR spreading compared to the conventional diffusion, ${z}^{\prime }\propto {t}^{1/2}$, that embodies a sub-diffusive propagation, ${z}^{\prime }\propto {t}^{1/4}$. This ameliorates the problem of acasual propagation in the diffusion regime, yet no sharp fronts or spikes develop. By contrast, the telegraph solution to the causality problem is to cut off the solution beyond certain distance ($| z| \gt \sqrt{k/T}t$), thus introducing an unphysical singularity. The immediately arising normalization problem is then "solved" by adding an even stronger singularity in form of two δ functions at the cut-off points.

Using the Chapman–Enskog method, we have extended the CR transport equation with magnetic focusing to the fourth order in a small parameter $\varepsilon =\lambda /l$ (CR mean free path to the characteristic scale of the problem). This analysis clarifies the nature of the telegraph transport equation, widely publicized in the literature as a promising alternative to diffusive propagation models. We have shown that the telegraph extension ($\propto {\partial }^{2}{f}_{0}/\partial {t}^{2}$) of the diffusion equation can be mapped from the (small) hyper-diffusive term ($\propto {\partial }^{4}{f}_{0}/\partial {z}^{4}$) of the regular Chapman–Enskog expansion, but the telegraph term, originating from an ${\varepsilon }^{4}$ term of the expansion, must remain subordinate to the main, diffusive transport and magnetic focusing (if present) contributions. This condition is met after $\nu t\gt z/\lambda$ collision times for an initially narrow (${\rm{\Delta }}z\lt \lambda$) CR distribution. Another important limitation of the telegraph equation is that, by contrast to the Chapman–Enskog equation and the original pitch-angle scattering equation, it is not self-contained requiring an initial condition for $\partial {f}_{0}/\partial t$ also, which needs information about the anisotropic part of the initial CR distribution. Furthermore, an attempt to proceed to higher orders in introduces progressively shorter timescales associated with "ghost" terms reflecting a quick relaxation of initial anisotropy or strong spatial inhomogeneity. By contrast, the classic Chapman–Enskog method is devised to eliminate short timescales, irrelevant to the evolution of the isotropic part of CR distribution ${f}_{0}$, which accurately describes this evolution after a few collision times, $\nu t\gt 1$.

We have derived the CR transport equation for an arbitrary pitch-angle scattering coefficient $D(\mu )$. This form of transport Equation (27), is suitable for describing CR acceleration and escape problems where the phenomenon of self-confinement (D is a functional of f, $D=D[f;\mu ,t]$) is critical (e.g., Ptuskin et al. 2008; Malkov et al. 2010b, 2013; Fujita et al. 2011). Accounting for magnetic focusing effects is required, e.g., for describing particle acceleration in CR-modified shocks with an oblique magnetic field. In this case, the field increases toward the shocks due to the pressure exerted by the accelerated CRs on the flow, thus producing a mirror effect. The particle drift velocity along the field associated with the mirror effect is (e.g., Equation (29)) $V\sim \kappa /{l}_{B}$ which, for the magnetic field variation scale being of the order of the shock precursor scale $\kappa /{U}_{\mathrm{sh}}$ and strong shock modification, almost automatically becomes comparable with the shock velocity ${U}_{\mathrm{sh}}$. This additional bulk motion of the accelerated CR (directed toward the shock) will affect their spectrum and acceleration time.

Furthermore, as the CR scattering in such environments (i.e., supernova remnant; SNR shocks) must be self-sustained by virtue of instabilities of the CR distribution (see, e.g., Bykov et al. 2013; Bell 2014 for the recent reviews), the above magnetic drift needs to be included in the CR stability analysis. It should be noted, however, that the results obtained in the present paper formally require a magnetic field ${B}_{0}$ that does not strongly change over the gyroradius of energetic particles. This is not to be expected in SNR shocks, especially if strong, CR current- and pressure-driven instabilities generate fields with $\delta B\gt {B}_{0}$. However, one may use the shock normal direction as the polar axis to calculate the pitch angle diffusion coefficient $D(\mu )$, needed for the description of the CR spatial transport. Also, such treatment will require a description of the gyro-motion and averaging by computing particle orbits beyond the standard quasi-linear description (Malkov & Diamond 2006), implied throughout this paper.

In conclusion, by comparison with the telegraph equation, the classic Chapman–Enskog hyper-diffusion equation consistently describes the long-term CR propagation in a self-contained, order-controlled fashion. Further improvement of the CR diffusion models should probably address the anisotropic component of the CR distribution. There are situations, such as ultra-high energy CR propagation, where the mean free path grows too long with the energy as to make the diffusive approach irrelevant and a rectilinear transport to dominate (Aloisio et al. 2009; cf. Levi flight regime described in aforementioned study by Malkov & Diamond 2006). Another interesting example is a sharp angular anisotropy $\sim 10^\circ$ in CR arrival directions discovered by the MILAGRO observatory (Abdo et al. 2008) and a number of other instruments, e.g., (Abbasi et al. 2011; Bartoli et al. 2013; Abeysekara et al. 2014; Desiati 2014). Unless this anisotropy is of a very local origin (such as the heliosphere; Lazarian & Desiati 2010; O'C. Drury 2013), it poses a real challenge to CR propagation models and clearly cannot be addressed within the diffusive approaches discussed in this paper (Drury & Aharonian 2008; Malkov et al. 2010a; Giacinti & Sigl 2012; Ahlers 2014; Malkov 2015). On the other hand, when the diffusive transport model is well within its validity range (weakly anisotropic spatially smooth CR distributions) neither the telegraph nor the hyper-diffusive term (both $\sim {\varepsilon }^{4}$) is essential to the CR transport and can be neglected.

Useful and stimulating comments of the Referee are gratefully acknowledged. The work of M.M. was supported by the NASA ATP-program under the grant NNX14AH36G. He is also indebted to the University of Maryland for hospitality and partial support during this work. Partial support from US DoE is also appreciated.

To derive the result given by Equation (37), we rewrite Equation (32) for a simple case of weak focusing ($\sigma f\lt {\partial }_{z}f$) and using the Galilean transform to the frame moving at the magnetic drift velocity ${V}_{2},$ ${z}^{\prime }=z+{V}_{2}\tau$:

$$\begin{eqnarray}&&\displaystyle \frac{\partial f}{\partial \tau }={\kappa }_{2}\displaystyle \frac{{\partial }^{2}f}{\partial {z}^{\prime 2}}-h\displaystyle \frac{{\partial }^{4}f}{\partial {z}^{\prime 4}},\end{eqnarray} \tag{ 38 }$$

where $h={\varepsilon }^{2}\left({K}_{1}-{K}_{2}\right)/2$ and ${\kappa }_{2}$ is the same as in Equation (32). The fundamental solution of Equation (38) can be written as an inversion of the Fourier image ${f}_{k}(\tau )$, assuming that ${f}_{k}(0)=1/2\pi$. Because we are going to find the Fourier inversion using asymptotic methods, we write an arbitrary constant C instead of this value:

$$\begin{eqnarray}&&f\left(\tau ,{z}^{\prime }\right)=C{\int }_{-\infty }^{\infty }{e}^{{{ikz}}^{\prime }-{k}^{2}\left({\kappa }_{2}+{{hk}}^{2}\right)\tau }{dk}.\end{eqnarray} \tag{ 39 }$$

We will determine the normalization constant C shortly from the requirement $f\to \delta \left({z}^{\prime }\right),\ \tau \to 0$. It is convenient to introduce the following notation:

$$\begin{eqnarray*}&&{\zeta }^{4}=4h\tau {k}^{4},\ \ \xi ={z}^{\prime }{(4h\tau )}^{-1/4}\gt 0.\end{eqnarray*}$$

We may limit our consideration to the $\xi \geqslant 0$ half-space, because the solution is an even function of z'. Focusing on a short time asymptotic regime $\tau \lt \varepsilon {z}^{\prime },$ which is opposite to the case considered earlier (Equation (36)), we neglect the diffusive term $\sim {k}^{2}$ in the exponent (as hyperdiffusion dominates) and rewrite Equation (39) as follows:

$$\begin{eqnarray}&&f\left(\tau ,{z}^{\prime }\right)=\displaystyle \frac{C}{{(4h\tau )}^{1/4}}{\int }_{-\infty }^{\infty }{e}^{i\xi \zeta -{\zeta }^{4}/4}d\zeta .\end{eqnarray} \tag{ 40 }$$

Note that the general evaluation of the integral in Equation (39), with the diffusive term included, is not difficult, but more cumbersome. The phase of the integral here has three saddle points, and the following two should be on the integration path

$$\begin{eqnarray*}&&{\zeta }_{\pm }=i{\xi }^{1/3}{e}^{\pm i\pi /3},\end{eqnarray*}$$

since the integrand reaches its maxima at these points. Thus, the integration runs from $-\infty$ through ${\zeta }_{+}$ to $+i\infty$ then to ${\zeta }_{-}$ and, finally to $+\infty$. The contributions from the two saddle points then yield

$$\begin{eqnarray}&&f=\displaystyle \frac{2C\sqrt{\pi /3}}{{(h\tau )}^{1/4}{\xi }^{1/3}}\mathrm{exp}\left(-\displaystyle \frac{3}{8}{\xi }^{4/3}\right)\mathrm{cos}\left(\displaystyle \frac{{3}^{3/2}}{8}{\xi }^{4/3}\right).\end{eqnarray} \tag{ 41 }$$

The normalization $\int {{fdz}}^{\prime }=1$ requires the constant $C=1/\sqrt{3}\pi$, which only insignificantly deviates from the exact value $C=1/2\pi$ that, in turn, follows from the integral representation of f in Equation (39) for $\tau \to 0$. Note that we could have replaced an oscillating exponential tail of this solution by zero beyond the point $\xi \gt {(4\pi )}^{3/4}{3}^{-9/8}$. Such modification of the hyperdiffusive solution would be in the spirit of the telegraph cut-off at $z/t=\sqrt{k/T}$, however, with an essential difference of being only a discontinuity in the solution derivative. The unphysical (oscillatory) behavior beyond the first zero point of the solution in Equation (41) results from neglecting the diffusion term in Equation (38), asymptotic methods used to calculate the integral in Equation (40), and from lacking higher order terms, $\sim {\varepsilon }^{n},\ n\gt 4$. Therefore, the solution can be improved systematically. In addition, it starts from a point source that is clearly inconsistent with the main approximation $\varepsilon \ll 1$. A somewhat broader initial profile will not develop an oscillatory tail, if convolved with the Green's function in Equation (39). We will not attempt to improve on this minor aspect of the solution here, as it becomes only weakly irregular if cut off at its first zero.

From the perspective of a general improvement of the asymptotic expansion considered in this paper, the derivation of Equation (34), for example, is robust in the following sense. The residual higher order terms in $\varepsilon \ll 1$, if included in Equation (32), will not change Equation (34) in any other way than small corrections to the coefficients ${\kappa }_{2}$ and ${V}_{2}$. Indeed, the higher z-derivatives coming from higher order terms, will vanish from the (first four) moment equations after integrating by parts. By contrast, continuing the telegraph approach to higher orders will generate terms with small parameters at higher time derivatives in all moment equations. These terms will become crucial during the initial relaxation of the CR distribution. The relaxation is associated with the CR anisotropy or strong initial inhomogeneity, that is with large $\tilde{{f}_{n}}$, Section 3. However, these decay over a short time $t\lesssim 1$. This is the time period when the telegraph or hyperdiffusive correction is large but its effect on the subsequent evolution ought to be limited because this time is short. The hyperdiffusive correction meets this requirement, as we argued using moment equations. To see whether the same is true for the telegraph correction, let us rewrite Equation (33) using the "fast" time, $t=\tau /{\varepsilon }^{2}$:

$$\begin{eqnarray}&&\left(1+{\tau }_{{\rm{T}}}\displaystyle \frac{\partial }{\partial t}\right)\displaystyle \frac{\partial {f}_{0}}{\partial t}=\displaystyle \frac{{\varepsilon }^{2}}{4}\kappa \displaystyle \frac{{\partial }^{2}{f}_{0}}{\partial {z}^{2}},\end{eqnarray} \tag{ 42 }$$

where we denoted ${\tau }_{{\rm{T}}}=2\left({K}_{1}-{K}_{2}\right)/\kappa \sim 1$ and assumed $\sigma =0$, to make the following simple argument. Namely, in the limit $\varepsilon \to 0$ there are two modes, of which the first is ${f}_{0}={f}_{0}(z)$. This is the main diffusion mode that slowly evolves in time when $0\lt \varepsilon \ll 1$, and, as we are interested in the evolution over the timescales $t\gtrsim {\varepsilon }^{-2}$, the telegraph term becomes $\sim {\varepsilon }^{4}$ and can be discarded. The second mode corresponds to a rapid decay of the initial distribution $\sim \mathrm{exp}\left(-t/{\tau }_{{\rm{T}}}\right)$ which is associated with the decay of initial anisotropy or strong inhomogeneity. If this mode is active ($\partial {f}_{0}/\partial t\ne 0$ in Equation (42)), then even the total number of particles N is not conserved automatically. Thus, turning to the moments of Equation (42), we need to impose the initial condition, $\partial N/\partial t=0$, to ensure the particle conservation. This probably means that $\varepsilon \to 0$ is a difficult limit for the telegraph reduction scheme. The singular components in the telegraph solution (6) appear to be primarily associated with the particle conservation problem. The initial relaxation phase ($t\lesssim 1$) perhaps, cannot be adequately described by the telegraph reduction scheme using an equation for ${f}_{0}$ alone because it does not properly "average out" an anisotropic component $\tilde{f}$, which is large during this period of time. The telegraph term is therefore to be understood as a "ghost" term reflecting rapid decay of such components. It follows that the rapidly changing part $\tilde{f}$ in the decomposition in Equation (16) needs to be retained in the short-time analysis along with ${f}_{0}$. Otherwise, the telegraph operator generates unphysical δ-pulses and sharp fronts, just to conserve the number of particles, as discussed earlier. These considerations are, however, not nearly complete. Further useful analysis of propagation modes in the context of the telegraph equation can be found in Schwadron & Gombosi (1994).

To conclude this appendix, we make yet another argument in disfavor of the telegraph equation that is partially related to the above considerations. A consistent asymptotic reduction method must be continuable to infinity in powers of small . The Chapman–Enskog scheme clearly is. The outcome will be a series of terms $\sim {\partial }_{z}^{n}{f}_{0}$ on the rhs of Equation (27) with just ${\partial }_{t}{f}_{0}$ on its lhs. To solve the resulting equation, only the initial distribution ${f}_{0}(0,z)$ is needed because the equation remains evolutionary and (generalized) parabolic, as its pitch-angle diffusion superset is. The telegraph equation, on the contrary, turns hyperbolic and non-evolutionary after the reduction from the superset equation. By continuing the telegraph reduction scheme to higher orders of approximation, progressively higher time derivatives will emerge (along with higher space derivatives). The resulting equations will thus be non-evolutionary and a growing set of initial time derivatives ${\partial }_{t}^{n}$ ${f}_{0}$ will then be needed to solve the initial value problem. These data can be extracted only from the full anisotropic distribution with recourse to the full (anisotropic) equation. Therefore, the telegraph equation is not self-contained and cannot be improved systematically. Any attempts to improve it will introduce shorter and shorter timescales that would require a return to the full anisotropic description.