An Exact Solution of the Fokker-Planck Equation for Isotropic Scattering

# An Exact Solution of the Fokker-Planck Equation for Isotropic Scattering

M.A. Malkov University of California San Diego, La Jolla, CA 92093
###### Abstract

The Fokker-Planck (FP) equation is solved analytically. Foremost among its applications, this equation describes the propagation of energetic particles through a scattering medium (in - direction, with being the - projection of particle velocity). The solution is found in terms of an infinite series of mixed moments of particle distribution, . The second moment () was obtained by G.I. Taylor (1920) in his classical study of random walk: (where is given in units of an average time between collisions). It characterizes a spatial dispersion of a particle cloud released at , with being its initial width. This formula distills a transition from ballistic (rectilinear) propagation phase, to a time-asymptotic, diffusive phase, . The present paper provides all the higher moments by a recurrence formula. The full set of moments is equivalent to the full solution of FP equation, expressed in form of an infinite series in moments . An explicit, easy-to-use approximation for a point source spreading of a pitch-angle averaged distribution (starting from i.e., Green’s function), is also presented and verified by a numerical integration of FP equation.

## I Introduction

Propagation of energetic particles, which we will call cosmic rays (or CR for short), particularly through magnetized turbulent media, has been actively researched in the astrophysical community for more than half a century. The time asymptotic solution of this problem is known to be diffusive. After several collisions particles “forget” their initial velocities and enter a random walk process. However, in astrophysical objects with infrequent particle collisions, there is not enough time or space for even a few collisions. In such systems, an early-time propagation is not random as particles “remember” their starting velocities and positions.

At times much shorter than the collision time, , most particles propagate with their initial velocities or their projections on the magnetic field direction, if present. This regime is called the ballistic, or rectilinear propagation. The question then is what happens next, namely at but before the onset of diffusion at ? What exactly is the value of , when the simple diffusive description becomes applicable? In other words, what is the extent of an intermediate phase when neither ballistic nor diffusive model applies? These are the questions we address below using a Fokker-Planck (FP) transport equation and its exact analytic solution.

A review of early results on CR propagation is contained, e.g., in Refs.Jok1971RvGSP (); Voelk73 (), while Ref.M_PoP2015 () covers some recent results relevant to the present paper. The FP model is extensively applied by the CR and heliophysics communities to a wide range of transport processes driven by small-step stochastic variations of particle velocity. This formulation is relevant, e.g., to the solar wind, through which solar energetic particles and CRs propagate to the Earth, and interstellar medium, through which CRs propagate from more distant sources. Under “collisions” in such media, one usually understands particle scattering off magnetic disturbances. Another example is the propagation of ultra-high-energy CRs from extragalactic sources. The transition from ballistic to diffusive transport regime, while being challenging for the theory, is key to understanding the nature of such sources. Since the particle mean free path usually grows with energy, part of their spectrum almost inescapably falls into a transient category where neither ballistic nor diffusive approximation applies. For the lack of better terms, we will call these particles transdiffusive.

Transdiffusive particles are likely to carry most of the information about their source. Indeed, the low-energy, diffusively propagating particles either do not reach us in time or merge into a featureless isotropic background. The highest energy particles, on the contrary, propagate ballistically. They may point back to their sources and are therefore invaluable. Unfortunately, they are exceptionally rare. On average, just one CR particle with energy eV is expected to arrive per century per square kilometer. Only a handful of such events has been registered over decades of observation. Their number grows towards lower energies, but their arrival direction becomes random too, thus making the source indiscernible on the sky. It follows that it is the transdiffusive propagation regime that may steer particles between Scylla of poor statistics and Charybdis of orbit scrambling and unveil the source. Understanding the transdiffusive propagation is then the key to an emerging field of CR astronomy.

## Ii Formulation of the problem and its significance

As there is no lack of motivation for studying the transdiffusive propagation regime, it is worthwhile to discuss the Fokker-Planck (FP) equation as a simple yet adequate mathematical model for this purpose. The FP equation is general in that it applies to both magnetized and unmagnetized media, as long as particles (or other entities, such as wave quanta) scatter randomly at small angles. The magnetic field, however, conveniently justifies the one-dimensional reduction of more realistic three-dimensional problems, as particles are usually bound to the field lines. After averaging out their gyromotion (typically unimportant) the particle phase space becomes two-dimensional. The resulting equation for their distribution function describes a one-dimensional spatial transport constrained by angular scattering:

 ∂f∂t+vμ∂f∂x=∂∂μ(1−μ2)D∂f∂μ. (1)

Here is the only spatial variable along which the particle concentration varies (local field direction or another symmetry axis), is the cosine of particle pitch angle to the axis, is the magnitude of particle velocity, conserved in interactions with quasi-static magnetic turbulence. is the scattering rate (collision frequency).

A major propagation scenario that Eq.(1) handles well comes about through an instant release of a small particle cloud into a scattering medium. For instance, galactic supernova remnants (SNR), widely believed to generate CRs with energies up to eV, must accelerate CRs in their shock waves and subsequently release them into the turbulent interstellar medium GrenierRev2015 (). Mathematically, the question then is how exactly the pitch-angle averaged particle distribution propagates along a magnetic flux tube that intersects the SNR shell. It is highly desirable 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 during the ballistic and transdiffusive propagation, while the latter is often the key for probing into the source. During these phases energetic CR protons (the main species) may reach a nearby molecular cloud, making themselves visible by interacting with its dense gas AharDV94 (). The CR protons of lower energies would instead be diffusively confined to the SNR shell and remain elusive for observers. Another example is the propagation of solar energetic particles to 1 AU. Also, in this case, the particle mean free path (m.f.p.) may be comparable to or exceed 1 AU, so the diffusion approximation fails again 111Note, however, that if the pressure of released particles is of the order or larger than the magnetic pressure of the ambient medium, the particles are self-confined by scattering off self-generated MHD waves. This kind of problems should be treated differently, MDS10PPCF (); MetalEsc13 (); NavaGabici2016 (); DAngeloBlasi2016 () from what is described further in this paper, as the scattering frequency would strongly depend on , whereas here we consider the case ..

### ii.1 Simple limiting cases

Before proceeding with the solution of the FP equation (1) in Sec.III, it is useful to characterize its limiting cases of ballistic and diffusive propagation. We deduce them directly from Eq.(1), by eliminating angular dynamics.

#### ii.1.1 Ballistic propagation regime

In the ballistic regime which strictly applies to times shorter than the collision time, , one can neglect the r.h.s. of the equation, and obtain the solution by integrating along the particle trajectories (Liouville’s theorem), with a conserved pitch angle, . The solution is simply .

It is sufficient to consider here a point source with initially isotropic distribution: where and denote the Dirac’s delta and Heaviside unit step functions, respectively. From the above solution for , one obtains the ballistic expansion in form of the second moment, by integrating over and . The result describes a free escape with the mean square velocity , while the maximum particle velocity (along ) is . The form of the pitch angle averaged particle distribution, , is characterized by an expanding ’box’ of decreasing height. Note that this result is inconsistent with the solution of the so called “telegraph” equation that has been put forward for CR propagation over the last 50 years and is discussed below at some length. By contrast, an exact solution of Eq.(1) obtained below converges to the above-described box distribution at .

#### ii.1.2 Diffusive (hyperdiffusive) propagation regime

The second simple propagation regime is diffusive which dominates at , and is treated in a way opposite to the above-described ballistic regime, Jokipii66 (). The r.h.s. of Eq.(1) is now the leading term, thus implying that the particle distribution is close to isotropy, . Working to higher orders in anisotropic corrections , and averaging the equation over , one obtains the following equation for MS2015 ()

 ∂f0∂t−κ2∂2f0∂x2=−κ4∂4f0∂x4+κ6∂6f0∂x6−…, (2)

with . This particular form of expansion is relevant under the scattering symmetry: . Otherwise, also the odd derivatives appear on its r.h.s. The latter situation is not considered here for simplicity. The last equation results from an asymptotic (Chapman-Enskog) expansion of the problem in . It is valid only for , 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) (“telegraph” term) which is illegitimate unless (see MS2015 (); M_PoP2015 () and below). Meanwhile, the r.h.s. of the above equation provides a small hyperdiffusive correction which may be omitted, as the higher spatial derivatives quickly decay because of the smoothing effect from the diffusive term on its l.h.s.

The most serious problem with the diffusive approximation is an unrealistically high propagation speed of particles which reach a given point faster than the maximum speed would allow (acausal solution, e.g. SchwadronTelegraph94 (); Berezinsky2007 (); AloisBerezSuperLum09 ()). Mathematically, the approximation violates an upper bound that immediately follows from Eq.(1) for a point source solution. There have been many attempts to overcome this problem, but no adequate ab initio description of particle spreading that would cover ballistic and diffusive phases was elaborated. Much popularity received the telegraph equation, mentioned above. Ian Axford AxfordTelegr65 () derived it semi-empirically and supported by studies of discontinuous random walk GOLDSTEIN01011951 () back in 1965. Many other derivations of telegraph equation from Eq.(1) and attempts to justify its application to various CR transport problems have appeared ever since, e.g., Earl1974 (); SchwadronTelegraph94 (); LitvinEffenSchlick15 (); Litvin2016PhPl (). However, the solution of the telegraph equation is arguably inconsistent with the FP equation. It is sufficient to mention a presence of singular components incompatible with the parent equation, eq.[1]. These are introduced to make up for nonconservation of the total number of particles which we also discuss below.

Both diffusive and “telegraph” approaches are aimed to extract the spatial particle distribution from Eq.(1). The elimination of pitch angle, however, has always been a challenge. While the time-asymptotic character of the diffusive approach is well understood, the idea behind the telegraph equation was to span also the transdiffusive evolution, preferably down to the ballistic phase. As we mentioned, the derivation of telegraph equation is inconsistent with a regular Chapman-Enskog asymptotic expansion (systematic elimination of pitch angle) and the ballistic propagation. The higher order Chapman-Enskog expansion, in turn, is valid only for time-asymptotic regimes, , and should be considered as a correction to the diffusive treatment. As in many other asymptotic expansions, when applied outside of their validity range, higher order terms often make the approximation less accurate which was recently demonstrated in Ref.Litvin2016PhPl (), using the numerical integration of Eq.(1). By contrast to the hyperdiffusive Chapman-Enskog expansion, the telegraph equation, as mentioned above, was intended to cover also the crossover phase, . Not surprisingly, its solution failed to provide an adequate fit to a full numerical solution, thus confirming its failure at .

Being most easily obtained from the hyperdiffusive corrections in Eq.(2), the telegraph equation inherits its validity range, . In this appearance, it constitutes just another form of small correction to diffusion. Indeed, considering the two terms on the l.h.s. of Eq.(2) as leading (which is required by its derivation!) and converting then the fourth spatial derivative (the leading term on the r.h.s) into a second time-derivative, we write: . By dropping higher - derivatives one recovers the telegraph equation

 ∂2f0∂t2−V2∂2f0∂x2+τ−1∂f0∂t=0. (3)

At first glimpse, it indeed captures a ballistic (wave-like) propagation of particle bunches at a speed . However, the number density in the bunches decays with time at a rate which is nonphysical. Unlike eqs.(1) and (2), this equation does not conserve the number of particles, , automatically. To conserve the total number of particles in a fundamental solution (Green’s function) for Eq.(3), two different solutions need to be added together. As the equation is linear, such addition is indeed legitimate. The first component of the solution is smooth within the characteristics of eq(3), and zero otherwise, thus it develops discontinuities at . The number of particles contained in this solution component grows in time from zero as they start spreading ballistically along the characteristics. The second solution, whose sole purpose is to compensate for the nonphysical multiplication of particles in the first component, needs to be taken in even more singular form, . The number of particles in this component decays as . However, this auxiliary component is inconsistent with the original equation (1). Indeed, the particle distribution of the form implies that all particles have the same pitch angle. But the operator of angular scattering on the r.h.s. of Eq.(1) would momentarily smear out any sharp pitch-angle distribution. Therefore, the components cannot persist in the parent FP equation and have been illegitimately acquired during its reduction 222Evidently, the above arguments do not apply to a telegraph equation derived from a Boltzmann equation with a simplified scattering operator of the form (so called BGK, or approximation). In this context, the telegraph equation has been studied in, e.g., Gombosi93 (); Zank2000 (); Webb2000 (); WebbTelegr06 (); Fedorov2016 ()..

We conclude that the telegraph equation can be accepted only at when the nonphysical - peaks die out. The telegraph equation can then be deduced from the hyperdiffusive expansion in Eq.(2), obtained, in turn, by using the canonical Chapman-Enskog approach. At earlier times, the telegraph solution does not match the actual solution of Eq.(1), e.g., Litvin2016PhPl (). Interestingly enough, only at merges the telegraph solution with diffusive, hyperdiffusive and direct numerical solution of the parent FP equation. Even then, nonphysical peaks stemming from the singular part of the solution remain well pronounced (Fig.5 in the above paper).

Notwithstanding the irrelevance of the telegraph solution to early phases of particle transport, singular - components do arise in a different context of the telegraph equation (apart from transmission lines, of course). It is associated with a discontinuous random walk. In this process, studied in Ref. GOLDSTEIN01011951 () and earlier by G.I. Taylor Taylor01011922 (), particles are allowed to move only at fixed velocities, positive or negative, say , which naturally results in a particle distribution. Statistically, this is similar to a coin tossing with only two possible outcomes. Under a continuous pitch-angle dependence relevant to the CR propagation, the discontinuous random walk restriction corresponds to the particle distribution, concentrated at , thus producing the singular components of the telegraph equation. As we noted, however, the underlying angular distribution is inconsistent with the solutions of FP Eq.(1). That is why the telegraph equation cannot be consistently derived from the FP equation except for . Being alternatively derived from the discontinuous random walk process, it is unsuitable for describing the transport of energetic particles since their spectrum in velocity projection on the travel direction (pitch angle) is a fundamentally continuous variable.

It follows that apart from the well established diffusive description of Eq.(1), though valid only for , there are no viable analytical tools to address the earlier phases of particle propagation. Therefore, an exact solution of the FP equation we tackle below is more than motivated.

## Iii Fokker-Planck Equation and its Solution

Turning to the solution of Fokker-Planck Eq.(1), we first discuss its scope. The energy dependence of the particle scattering frequency enters it only as a parameter, i.e., , and does not prevent us from solving the equation. By contrast, a possible pitch-angle dependence of (nonisotropic scattering) does. However, in media with magnetic irregularities, derives from a power index of the scattering turbulence, , if the interaction between particles and turbulence is resonant, e.g. VVSQL62 (); Kennel66 (); Jokipii66 (); Skill75a (); BlandEich87 (). In particular, for a power spectrum , where is the wave number, one obtains and the scattering is isotropic for an important case , as discussed below. Even more complex, anisotropic turbulence spectra, such as those derived by Goldreich and Shridhar goldr97 (), result in a flat ChandranGS00PhRvL (), except for relatively narrow regions near . These areas require special considerations in any event, as they are strongly affected by particle mirroring and resonance broadening Acht81a () (), as well as by the field aligned propagation Milagro10 () (. We may ignore them as they occupy only a small fraction of particle phase space. Besides, fluctuation spectra with the index have been obtained in Monte-Carlo studies of shock-accelerated particles Bykov2014ApJ (). The choice of -independent scattering coefficient has been advocated in Shalchi2009 () even for strong magnetic fluctuations, , where is the unperturbed field. Motivated by the above, we consider the case as it embraces many physically interesting situations and, at the same time, allows for an exact analytic solution of Eq.(1).

We now rewrite Eq.(1) using dimensionless time and length units according to the following transformations

 D(E)t→t,Dvx→x (4)

Instead of Eq.(1) we thus have

 ∂f∂t+μ∂f∂x=∂∂μ(1−μ2)∂f∂μ (5)

This equation contains no explicit parameters thus precluding any uniformly valid asymptotic expansion unless a small parameter enters the problem implicitly through the initial condition . In particular, if one is interested in an isotropic approximation that can be treated using Eq.(2), (- type expansion), not only should the initial distribution be close to isotropy, , but it should also be spatially broad, . The latter condition prevents a high anisotropy from arising in the course of time via the second term on the l.h.s. of Eq.(5). Hence, the fundamental problem of a point source spreading (Green’s function, or fundamental solution) can not be treated using a conventional expansion, until becomes quasi-isotropic, that is broadened to . This is another reason why the telegraph equation falls short in describing CR propagation from a point source.

Based on the above considerations, we tackle an exact solution of Eq.(5). The only restriction that we impose on the spatial distribution at , which holds up during its subsequent evolution, is a sufficiently rapid decay of at . Namely, we require that for and . This standard restriction guarantees the existence of all moments. Another standart restriction is the regularity of at : for .

Turning to the solution of Eq.(5) we introduce moments of in the form of the following matrix

 Mij(t)=⟨μixj⟩=∫∞−∞dx∫1−1μixjfdμ/2 (6)

for any integer . We will discuss conditions for the equivalence of the moments and the distribution when the solution for the matrix is obtained.

The lowest moment is automatically conserved by Eq.(5) (as being proportional to the number of particles) and we normalize it to unity

 M00=∫∞−∞dx∫1−1fdμ/2=1.

Multiplying Eq.(5) by and integrating by parts (where appropriate), we obtain the following matrix equation for the moments with :

 ddtMij+i(i+1)Mij−jMi+1,j−1=i(i−1)Mi−2,j (7)

This equation couples triads of matrix elements. The two elements are on the same antidiagonal (l.h.s), and the third one is on the next but one antidiagonal above it (r.h.s). One may get an impression that this infinite system of moment equations will require truncation to become useful. Actually, not only is this unnecessary but should be avoided at all cost, as we argue below. The set of moments in Eq.(7) can, in fact, be recursively obtained to any order with no truncation. Indeed, as and we can set for any , all the elements of the infinite matrix

on each antidiagonal can be found, working from left to right as shown by the arrows. Alternatively, the matrix of moments can be viewed as a triangle (akin to Pascal’s or Bernoulli’s triangles, for example) with its rows being the matrix antidiagonals, starting from the unity at the top of the triangle. The two moments on the next antidiagonal (triangle row), are easily found from Eq.(7) to be and Higher moments can be obtained inductively. So, in general, from Eq.(7) we find

 (8)

As may be seen from the above expressions for and all the higher moments in Eq.(8) are series in , where and are integral numbers. In particular, the next set of moments is on the third anti-diagonal, :

 M20=13,M11=16(1−e−2t)M02=M02(0)+t3−16(1−e−2t) (9)

Interesting in a point source (fundamental) solution, we assume the initial distribution to be symmetric in and isotropic which eliminates the odd moments. Furthermore, the initial spatial width must then also be set to zero, . Note that in Eq.(9) coincides with the respective random walk result obtained by G.I. Taylor Taylor01011922 (), which we discussed earlier. However useful for understanding the transition between ballistic and diffusive phases of particle propagation, alone does not, of course, resolve the FP equation. From the mathematical point of view, only a full set of moments in Eq.(8) provides a complete solution of Eq.(5) given the initial value, that determines the matrix in Eq.(8). Moreover, to adequately reproduce the ballistic and transdiffusive phases the series of moments cannot be truncated.

In general, the equivalence between an arbitrary distribution and its full set of moments is not guaranteed automatically, but can be established for Eq.(5) with its solution in the form of Eq.(8) using Hamburger’s theorem, e.g., reed1975methods (). The theorem assumes upper bounds on the moments in the form with the constants and being independent of . According to Eq.(8), for any fixed the moments grow with not faster than . Although for small higher powers of are present, they also have an upper bound . Therefore, the condition for Hamburger theorem is satisfied.

Being interested in the fundamental solution of the FP equation, we will focus on the moments that correspond to the isotropic part of particle distribution

 f0(x,t)=∫1−1f(μ,x,t)dμ/2, (10)

as only this part contributes to the particle number density. Note that under the fundamental solution we understand here the solution for that it isotropic at and . It constitutes the pitch-angle averaged distribution and has been the target of most reduction schemes applied to the FP Eq.(5). The matrix elements that represent are, therefore, Note, that with are nevertheless not small and remain essential for calculating the full set of the moments . To link them to , we use a standard moment-generating function

 fλ(t)=∫∞−∞f0(x,t)eλxdx=∞∑n=0λ2n(2n)!M0,2n(t) (11)

where we omitted the odd moments irrelevant to the fundamental (symmetric in ) solution. The above expansion will be converted into a Fourier transform of by setting . To find its inverse, that is the function , we will use the inverse Fourier transform

 f0(x,t)=12π∫∞−∞dkeikx∞∑n=0(−1)nk2n(2n)!M0,2n(t) (12)

To lighten the algebra, we will continue to use instead of for a while. For practical use, we need to simplify the series entering each moment , starting from those shown in Eq.(9). The higher moments (a few of them can be found in Appendix A) contain more terms and quickly become unmanageable without computer algebra. In the next section, we will sum up the series in Eq.(12) by extracting the dominant terms from each moment in the sum, depending on . It is crucial to sum up all the moments with no truncation, as we mentioned and will further discuss in the next section 333Some of the earlier studies of FP equation also pursued moment approaches, e.g., Shalchi2006 (). However, a set of moments similar (to ) has been truncated typically at or no larger than two. Besides, no spatial profile for similar to Eq.(12) has been given in the above reference, so we make no comparison with our results. .

To emphasise the role of the anisotropic part of particle distribution, we provide an explicit expression for the moments in Eq.(12) through the lower order, dipolar, moments

 M0,2n(t)=2n∫t0M1,2n−1(t′)dt′ (13)

The isotropic part of the solution of Eq.(5) is thus given by eqs.(12), (13) and (8).

## Iv Simplified forms of the solution

Equations (12) and (13) provide an exact closed form solution of the Fokker-Planck Eq.(5). The calculation of the moments is relatively straightforward. Using Eq.(8) and integrating by parts one obtains to any order in form of polynomials in and : with the constant matrix elements that can also be recursively obtained from Eq.(8). However, the expressions for grow rapidly in length with , and some computer algebra is virtually required to calculate the series in the Fourier integral in Eq.(12). Therefore, for practical use, a simplified approximation of the series of Fourier integrals is desirable. Before turning to its derivation, we note that the CR distribution develops two moving sharp fronts in the profile , as discussed in Sec.II.1.1. Sharp fronts are dominated by the contributions from in the Fourier integral given by Eq.(12). This reaffirms our earlier statement that the series in should not be truncated at any finite .

### iv.1 Ballistic and transdiffusive phases

We begin summing up the infinite series entering the moment generating function given by Eq.(11) for time relevant for a ballistic and (poorly understood) transdiffusive propagation. Each term of the series is calculated using a three-term expansion in powers of at each power of :

 fλ=1+λ22!t23(1−23t+t23+…)+λ44!t45(1−43t+5845t2+…)+λ66!t67(1−63t+4315t2+…)+… (14)

The first two terms in all parenthetical expressions suggest to introduce the following variable instead of :

 t′=t(1−t3). (15)

The “retarded time” is related to the second moment as it slows down similarly during the transition from the ballistic to transdiffusive phase, , for . We will use this relation between and later. Meanwhile, the two leading terms in at each power of in Eq.(14) can be written as powers of

 λntn(1−n3t+…)≈λnt′n

and summed up straightforwardly for . The remaining terms ( in the parentheses in eq.[14]) can also be summed up (see Appendix A). To the same order in , also valid for large but still restricted by , we rewrite the series in Eq.(14) as follows

 fλ≈1λt′sinh(λt′)eλ2Δ2/4 (16)

where

 Δ={2t′1/2t3/2/3√3,λt≪12t′1/2t3/2/3√5,λt≫1 (17)

The full coverage of for is vital to a correct description of sharp fronts, implying ( in Fourier integral). The additional variable in the above solution changes between its limits at and only insignificantly which offers an opportunity to use the moment generating function in Eq.(16) for a full description of the solution for , including the stage when the sharp fronts smear out (see below). However, the transition between the two limiting cases is important for such description and discussed further in Sec.IV.3.

Using the Fourier spectral parameter and performing the inverse Fourier transform according to Eq.(12), from Eq.(16) we obtain

 f0(x,t)≈14t′[erf(x+t′Δ)−erf(x−t′Δ)] (18)

Here erf stands for the error function and, again, the essential steps of the derivation can be found in Appendix A. By construction, this simple result constitutes an approximate, pitch angle averaged Green’s function for Eq.(5), which, as it should, satisfies the condition for . As we mentioned and will argue further In Sec.IV.3, Eq.(18) reaches far beyond its formal validity range, , but some conclusions can be drawn immediately from its present form.

The spreading of an infinitely narrow initial peak, , proceeds in the following phases. During a ballistic phase of propagation, that is strictly at , the particle distribution is best represented by an expanding ’box’ of the height in the region . Its edges propagate in opposite directions along the ’trajectories’, The box walls are initially much thinner than the box itself, , since . However, as the expansion progresses the box walls thicken to its total size and the process transitions into a diffusive phase that we consider in the next section. As the expansion progresses in a retarded time, , it slows down. Evidently, this process cannot be followed far enough in time using Eq.(18), already for the reason that becomes negative for . Clearly, this problem occurs from the limitation of the approximate summation of the series in Eq.(14). As we will show in Sec.IV.3, Eq.(18) has a potential for much better approximation of the exact solution at , if supplemented by more accurate expressions for and . To better understand the transdiffusive regime at we first consider the opposite limit of diffusive particle propagaion at .

### iv.2 Transition to diffusive propagation

For all terms containing powers of in the expansion given by Eq.(11) can be discarded and only the highest powers of need to be retained. Upon extracting such terms from each by using the solution for the moments in eqs.(8) and (13) (see also the expressions of the first few moments in Appendix A), we may write the moment generation function in Eq.(11) as

 fλ(t)=∞∑n=0λ2n(2n)!M0,2n(t)≈∞∑n=0(2n−1)!!(2n)!λ2n(t3)n (19)

By writing , the latter series can be summed up straightforwardly to yield

 fλ=eλ2t/6 (20)

After replacing and performing an inverse Fourier transform we obtain the conventional diffusive solution

 f0(x,t)=12π∫∞−∞dkeikx−k2t/6=√32πte−3x2/2t (21)

Now that we have obtained for both small and large , we turn to the most interesting transdiffusive regime between these two limiting cases. Before doing so, we note that the last expression for , valid for , can be improved by adding the next order terms to the series in Eq.(19). At the same time, the and results given by eqs.(18) and (20), respectively, already allow us to synthesize them into an approximate, uniformly valid (for all ) propagator. We, therefore, consider it in the next subsection and defer a rigorous matching of the three propagation regimes to a future study.

### iv.3 Unified CR Propagator

Although our derivation of the simplified CR propagator in Eq.(18) relies on the condition , this condition can be relaxed, provided that the two time-dependent variables of the propagator, and , are properly redefined. As we will see, the propagator then accurately describes the solution for all , including . To examine this premise we rewrite Eq.(18) accordingly

 f0(x,t)≈14y[erf(x+yΔ)−erf(x−yΔ)]. (22)

So far, we only know and at (see Sec.IV.1), so they are left to be determined for . As may be seen from Eq.(22), are the coordinates of two fronts propagating in opposite directions from the initial pulse at , while is the growing front width. Our hope that Eq.(22) approximates the true solution better than its prototype in Eq.(18), is based on the following observations:

1. It recovers a correct expansion for small , both for sharp and smooth fronts ( and cases), considered in Sec. IV.1

2. It recovers the correct solution for in Eq.(21).

3. It conserves the total number of particles, for arbitrary and

The last statement can be verified by a direct substitution of from Eq.(22)

 ∫∞−∞f0(x,t)dx=12√πy∫∞−∞dx[∫(x+y)/Δ0−∫(x−y)/Δ0]e−x′2dx′=1,

with the last relation following, e.g., from doing the outer integral by parts in and then changing the integration variable in the both remaining integrals, . The result (1) was derived in Sec.IV.1. The statement (2) can be verified by expanding the r.h.s. of Eq.(22) in small which recovers the diffusive solution in Eq.(21) under the condition .

Based on the above considerations and the results of Secs.IV.1 and IV.2, we require and to behave at small and large , respectively, as follows:

 Δ=⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩23√5t2,t≪1,ballistic23√3t2,t≪1,transdiffusive√2t/3,t≫1,diffusive (23)
 y={t′,t≪1o(√t),t≫1

The - values for the ballistic and transdiffusive regimes (both at differ from each other only by a factor of . Furthermore, a trend of changing the and time dependencies from, respectively, and to a slower growth at larger has already been rigorously established in Sec.IV.1 in form of a retarded time, . For larger this relation should be improved by summing terms with higher powers of in Eq.(14). However, an obvious similarity of the time dependence of and suggests that the moment (and possibly some higher moments) may better describe the quantities and entering the propagator in Eq.(22). Indeed, as opposed to , for example, that has been consistently derived only for , the moments are calculated exactly for arbitrary . The key observation here is that the propagator in Eq.(12) depends on time not explicitly but only through the moments, . As they all monotonically grow with time, a single moment, such as , suffices to describe the time dependence of the propagator. On a practical note, all these moments are built upon the powers of and , so that they all can be explicitly expressed through and from the moment generating function, given by Eq.(11). Another promising avenue to explore along these lines in a future study is to convert the moment expansion in Eq.(11) into a cumulant expansion.

In this paper, we will use only the second moment to extend our expressions for and , obtained for and , over the remainder of the time axis. The particular choice of is supported by its ability to correctly (exactly) describe the particle spatial dispersion over all the three stages of their propagation. As we are considering the Green function solution, we must set , so the moment takes the following simple form

 M02=t3−16(1−e−2t) (24)

Obviously, we cannot recover all the three cases of in Eq.(17) by making in the above formulae simply proportional to , which is correct only for . For , that is in Eq.(21), should be taken proportional to instead. The transition from ballistic to transdiffusive phase (both within ) can, in turn, be calculated by applying the saddle point method to the Fourier integral in Eq.(11). Again, this exercise is planned for a separate publication. Here we take a practical approach based on a simple interpolation between different and , and verify the results numerically. So, based on our representation of in Eq.(17) we use the following fit, that should work for not too large

 Δ=2√5M02(t){1+18[1+tanh(t−ttrΔt)]} (25)

Here and are the transition time and its duration between ballistic and trans-diffusive phases. The factor approximately corresponds to the depth of the transition (between two top lines in Eq.(23)). As for , a simple choice consistent with our result would be . However, since should decay faster than for , we introduce a small correction at

 y=√3M02[1−A(t−t′tr)Θ(t−t′tr)] (26)

where denotese the Heaviside unit function.

The exact solution of the FP equation is still given in the form of an infinite series in Eq.(12), containing terms obtained recursively. It can hardly be reduced to an expression simple enough and capable of reproducing the solution with very high accuracy at the same time. To attain an arbitrary accuracy, the series needs to be summed up numerically. The approximate propagator in Eq.(22) appears as a plausible alternative in this regard, as it is almost as simple as the standard Gaussian propagator. To demonstrate that it is also sufficiently accurate, we plot from Eq.(22) with its two input parameters, and , from eqs.(25) and (26), respectively. This choice leads to a somewhat less accurate result for than the rigorously obtained and in eqs.(17) and (18), but we adhere to Eq.(22) because it is accurate for all .

Shown in Fig.1 is the FP solution at different times, starting from an infinitely narrow initial pulse, . The solution profiles are plotted using the simplified analytic propagator from Eq.(22), the standard diffusive solution from Eq.(21), and a direct numerical integration of the FP equation, as described in Appendix B. By its derivation, the analytic propagator is expected to be accurate for , which is successfully confirmed by its comparison with the numerical results. By contrast, the diffusive solution is particularly inacurate and acausal at the early times of evolution which was also not unexpected. Somewhat surprising is that its deviation from the exact solution continues after quite a few collisions, that is at .

The analytic propagator, on the contrary, demonstrates a very good agreement with the numerical solution over the total integration time 444Moreover, it successfully reproduces the true FP solution for all , when an obvious additional constraint is imposed in Eq.(26). On the other hand, all the three solutions merge into one Gaussian after , Fig2, so there is no need to discuss this range. . It does show a rather expected but surprisingly minor deviation at , given that no efforts have yet been made to address the overlap region between the and rigorous approximations, using which the series in the exact solution has been summed up. By considering the overlap region, the agreement can be systematically improved which we plan for a future study. The diffusive regime begins surprisingly late. Even at , the diffusive solution remains noticeably acausal.

Surface plots of afford additional insight into the particle propagation. Fig.3 shows them for the three functions from Fig.1. The surface plots have an advantage of conveying the information about the time history of particle intensity observed at a fixed distance from the source. To obtain this time history, it is sufficient to make a slice of the surface plot along the plane It should be remembered, however, that both and variables entering the solution for also depend on particle energy (velocity), Eq.(4). With this in mind, the time histories obtained for different may be compared, for example, with the measurements of solar energetic particles at 1AU for different energies Roelof+2002ApJ (); Klassen2016A&A (). Propagation of nearly relativistic electrons is suggested, by these and several other observations, to be in a scatter-free (ballistic) regime, which is supported by a highly anisotropic, beam-like, electron distributions. From the standpoint of the FP solution obtained in this paper, such regime takes place over times in Fig.3. The slices along the planes, that show at relatively small are characterized by a very steep rise and long decay after a maximum which is consistent with the observations. However, an inspection of the anisotropic components of distribution shows that they still exceed the value of near the sharp fronts on the distribution profile even at times somewhat beyond , as may also be inferred from Fig.1. In this situation, it might be more appropriate to speak about the transdiffusive rather than ballistic propagation. Also, anisotropic components of the distribution need to be extracted from the exact solution for a meaningful comparison with the observations.

## V Conclusions

In this paper, an exact solution for a Fokker-Planck equation, given in the form of Eq.(1), is obtained. The primary focus has been on the isotropic part of the particle distribution describing an evolution of particle number density during their propagation from a point source (Green’s function). The propagation is fundamentally simple. It can be categorized into three phases: ballistic (, transdiffusive () and diffusive (), with being the collision time. In the ballistic phase, the source expands as a “box” of size with gradually thickening “walls.” The next, transdiffusive phase is marked by the box’s walls thickened to a sizable fraction of the box and its slower expansion, Figs. 1 and 3. Finally, the evolution enters the conventional diffusion phase, Fig.2.

The exact closed-form solution of Eq.(5) is given by Eq.(12) in form of an infinite series in moments of particle distribution that are, in turn, obtained recursively from eqs.(8) and (13). Eq.(22) and Appendix A provide simplified summation formulas for the series. They adequately describe a point source spreading through ballistic, transdiffusive and diffusive phases of particle transport.

No signatures of a well-known solution to the telegraph equation with the same initial condition (e.g., Litvin2016PhPl ()) are present in the exact solution. The signatures are expected in the form of two sharp peaks attached to the oppositely propagating fronts. Their absence confirms the earlier conclusion MS2015 () that the telegraph equation is inconsistent with its parent Fokker-Planck equation except for a late diffusive phase ( During this phase of particle propagation, however, a common diffusive reduction of Fokker-Planck equation, e.g. Jokipii66 (), suffices.

Signatures of superdiffusive propagation regime are also absent in the exact solution. Such regime is often postulated, e.g., in studies of diffusive shock acceleration (DSA) Perrone2013SSRv (), in the form of a power-law dependence of particle dispersion , with . According to an exact result shown in Eq.(24), the time dependence smoothly changes from the ballistic, , propagation to diffusive one, , without dwelling at any particular value of inbetween. This only means, of course, that a simple small-angle scattering model behind the FP equation does not lead to the superdiffusive transport. More complicated scattering fields in shock environments, e.g., MD06 (), may result in both superdiffusive (Lévy flights) and subdiffusive (long rests) transport anomalies MetzlerKla2000PhR (). They are, however, not generic to the DSA (unlike, e.g., Bohm regime) and should be justified on a case-by-case basis.

###### Acknowledgements.
I’m indebted to R. Blandford, A. Bykov, and V. Dogiel, along with other participants at the second San Vito Conference “Cosmic Ray Origin beyond the Standard Model,” held in September 2016 in Italy, for interesting discussions of the present results. A useful exchange of views on the telegraph equation with Yu. Litvinenko is gratefully acknowledged as well. I also thank two anonymous referees for furnishing additional references. This work was supported by the NASA Astrophysics Theory Program under Grant No. NNX14AH36G.

## Appendix A Higher moments of particle distribution and summation formulae

Using Eq.(8), after some computer-assisted algebra, we obtain the following expressions for a few higher moments, , needed to compute the full solution given by Eq.(12)

 M0,4=1270e−6t−t+25e−2t+13t2−2645t+107270
 M0,6=131500e−12t−3t+21134e−6t+e−2t×(310t2+639350t+47431750)+59t3−3718t2+22663t−61432268
 M0,8=16945750e−20t−5t+2253125e−12t+(t2567+11t11907−5927783)e−6t−(1425t3+858125t2+1510425625t+18509371506250)e−2t+3527t4−22427t3+3554135t2−2811836075t+1234033375

This process can be continued using the recurrence relations in eqs.(8) and (13). Our purpose, however, is to derive a simplified, easy-to-use version of the complete solution. Considering the case , we first expand the moments up to the order This gives us a series in Eq.(14) for the moment-generating function in which we slightly rearrange the terms as follows

 fλ=1+λ2t23![(1−t3)2+2t29…]+λ4t45![(1−t3)4+2845t2+…]+λ6t67![(1−t3)6+65t2+…]+… (27)

This representation of the series suggests using a retarded time instead of . Separating then the terms in the brackets to form two individual series, of which the first one can be summed up immediately, we obtain

 fλ=12λt′(eλt′−e−λt′)+t245S(λt) (28)

The remaining series can be represented in the following compact form

 S(y)≡∞∑n=12n(2n+3)(2n+1)!y2n

and summed up by rewriting it as

 S=yddy1y2ddy∞∑n=1(y)2n+3(2n+1)!=yddy1y2ddyy2(sinhy−y)

Using this formula, we evaluate Eq.(28) to its final form

 fλ(t)=1λt′sinh(λt′)+t245[2cosh(λt)+(λt−2λt)sinh(λt)] (29)

It is important that can be arbitrarily large here, , even though our expansion technically requires . As we emphasized earlier, large values of are associated with the contribution of higher moments in the series in Eq.(27) which, in turn, is responsible for sharp fronts.

The leading term in Eq.(29) is the first one, while from the second one we extract a contribution proportional to . First, we assume which applies to sharp fronts and narrow initial CR distributions. Although the ballistic phase precedes the transdiffusive one, the value of is effectively larger namely during the ballistic phase. It decreases while the fronts become smoother with increasing time but with the condition held up.

To the same order of expansion in , Eq.(29) can then be written as ()

 fλ(t)≈1λt′sinh(λt′)(1+λ2t3t′45)≈1λt′sinh(λt′)eλ2t3t′/45 (30)

Turning to the opposite case, , instead of the last expression we have

 fλ(t)≈1λt′sinh(λt′)eλ2t3t′/27 (31)

Because of the approximation , the product of and in these formulae can be replaced by a different combination , as it enters a correction term. From the standpoint of asymptotic expansion in use, there is no difference between and in the highest order terms included. We used this aspect in obtaining the series expansion in Eq.(31). It is a remarkable feature of the above expansion of that it has a very similar form for small and large values of . Even more interesting, perhaps, is that the solution , obtained below by inverting the Fourier integral from the above relations, can be written using the same functional dependence of the solution on the relevant time-dependent variables in all three phases, Eq.(22).

Replacing with and substituting Eq.(30) into the inverse Fourier integral in Eq.(12), we obtain

where . By performing a straightforward integration the result given in Eq.(18) immediately follows. Note that the last relation is obtained for the case , while in the opposite case Eq.(31) should be used which simply results in the replacement in the exponent.

## Appendix B Numerical verification of the analysis

Regarding a full expansion of particle distribution in Legendre polynomials

 f(x,μ,t)=∞∑n=0fn(x,t)Pn(μ), (32)

so far we have focused on as it is an isotropic and only constituent of the particle distribution that contributes to their density. This does not mean, of course, that with are unimportant. During the ballistic and transdiffusive propagation phases, a dozen of the first are generally of the order of, or even larger than, especially near the two sharp fronts in Fig.1.

Substituting the above expansion into the FP equation (5) one obtains the following system for

 ∂fn∂t=−n(n+1)fn−n2n−1∂fn−1∂x−n+12n+3∂fn+1∂x+ϵ∂f2n∂x2 (33)

which we solve numerically for under the following conditions

 fn≡0,forn<0,n≥nmax
 f0(x,0)=1√πδe−x2/δ2,fn(x,0)=0,n>0

A small diffusive term is added to the r.h.s for wellposedness. Keeping makes the integration results practically insensitive to this parameter, even for rather steep initial conditions (small ). The optimum truncation parameter also depends on the scale of initial condition, for which we take a Gaussian with the width for and zero for . This setting is consistent with a Green function solution sought for . Under these conditions, the optimum , beyond which the integration results do not change noticeably. The integration domain in is taken large enough to ensure the conditions (at the maximum of ). The choice of depends, of course, on the integration time, as seen from Figs.1 and 2.

The system in Eq.(33) was integrated numerically for , in Eq.(32) using an adaptive mesh refinement (collocation) algorithm, described in, e.g., MuirBacoli2004 (). It is especially well suited for evolving sharp fronts that are formed during the ballistic and transdiffusive particle propagation regimes.

## References

• (1) J. R. Jokipii, Reviews of Geophysics and Space Physics 9, 27 (1971).
• (2) H. J. Völk, Astroph. and Space Sci.25, 471 (1973).
• (3) M. A. Malkov, Physics of Plasmas 22, 091505 (2015).
• (4) I. A. Grenier, J. H. Black, and A. W. Strong, Ann. Rev. Astron. Astroph. 53, 199 (2015).
• (5) J. R. Jokipii, Astrophys. J. 146, 480 (1966).
• (6) F. A. Aharonian, L. O. Drury, and H. J. Voelk, Astronomy and Astrophys.285, 645 (1994).
• (7) Note, however, that if the pressure of released particles is of the order or larger than the magnetic pressure of the ambient medium, the particles are self-confined by scattering off self-generated MHD waves. This kind of problems should be treated differently, [38-41] from what is described further in this paper, as the scattering frequency would strongly depend on , whereas here we consider the case .
• (8) M. A. Malkov and R. Z. Sagdeev, Astrophys. J. 808, 157 (2015).
• (9) N. A. Schwadron and T. I. Gombosi, J. Geophys. Res. 99, 19301 (1994).
• (10) V. Berezinsky and A. Z. Gazizov, Astrophys. J. 669, 684 (2007).
• (11) R. Aloisio, V. Berezinsky, and A. Gazizov, Astrophys. J. 693, 1275 (2009).
• (12) W. I. Axford, Plan. Space Sci. 13, 1301 (1965).
• (13) S. Goldstein, The Quarterly Journal of Mechanics and Applied Mathematics 4, 129 (1951).
• (14) J. A. Earl, Astrophys. J. 188, 379 (1974).
• (15) Y. E. Litvinenko, F. Effenberger, and R. Schlickeiser, Astrophys. J. 806, 217 (2015).
• (16) Y. E. Litvinenko and P. L. Noble, Physics of Plasmas 23, 062901 (2016).
• (17) Evidently, the above arguments do not apply to a telegraph equation derived from a Boltzmann equation with a simplified scattering operator of the form (so called BGK, or approximation). In this context, the telegraph equation has been studied in, e.g., [42-46].
• (18) G. I. Taylor, Proceedings of the London Mathematical Society s2-20, 196 (1922).
• (19) A. A. Vedenov, E. P. Velikhov, and R. Z. Sagdeev, NUCLEAR FUSION, Suppl.2, 465 (1962).
• (20) C. F. Kennel and F. Engelmann, Physics of Fluids 9, 2377 (1966).
• (21) J. Skilling, Mon. Not. R. Astron. Soc. 172, 557 (1975).
• (22) R. Blandford and D. Eichler, Phys. Rep. 154, 1 (1987).
• (23) P. Goldreich and S. Sridhar, Astrophys. J. 485, 680 (1997).
• (24) B. D. G. Chandran, Physical Review Letters 85, 4656 (2000).
• (25) A. Achterberg, Astronomy and Astrophys.98, 161 (1981).
• (26) M. A. Malkov, P. H. Diamond, L. O’C. Drury, and R. Z. Sagdeev, Astrophys. J. 721, 750 (2010).
• (27) A. M. Bykov, D. C. Ellison, S. M. Osipov, and A. E. Vladimirov, Astrophys. J. 789, 137 (2014).
• (28) A. Shalchi, T. Skoda, R. C. Tautz, and R. Schlickeiser, Astronomy and Astrophys.507, 589 (2009).
• (29) M. Reed and B. Simon, Methods of modern mathematical physics, Vol. II (Academic Press, New York, 1975).
• (30) Some of the earlier studies of FP equation also pursued moment approaches, e.g., [47]. However, a set of moments similar (to ) has been truncated typically at or no larger than two. Besides, no spatial profile for similar to Eq.(12) has been given in the above reference, so we make no comparison with our results.
• (31) Moreover, it successfully reproduces the true FP solution for all , when an obvious additional constraint is imposed in Eq.(26). On the other hand, all the three solutions merge into one Gaussian after , Fig2, so there is no need to discuss this range.
• (32) G. M. Simnett, E. C. Roelof, and D. K. Haggerty, Astrophys. J. 579, 854 (2002).
• (33) A. Klassen et al., Astronomy and Astrophys.593, A31 (2016).
• (34) D. Perrone et al., Space Sci. Rev. 178, 233 (2013).
• (35) M. A. Malkov and P. H. Diamond, Astrophys. J. 642, 244 (2006).
• (36) R. Metzler and J. Klafter, Phys. Rep. 339, 1 (2000).
• (37) R. Wang, P. Keast, and P. Muir, Journal of Computational and Applied Mathematics 169, 127 (2004).
• (38) M. A. Malkov, P. H. Diamond, and R. Z. Sagdeev, Plasma Physics and Controlled Fusion 52, 124006 (2010).
• (39) M. A. Malkov et al., Astrophys. J. 768, 73 (2013).
• (40) L. Nava et al., Mon. Not. R. Astron. Soc. 461, 3552 (2016).
• (41) M. D’Angelo, P. Blasi, and E. Amato, Phys. Rev. D 94, 083003 (2016).
• (42) T. I. Gombosi et al., Astrophys. J. 403, 377 (1993).
• (43) G. P. Zank, J. Y. Lu, W. K. M. Rice, and G. M. Webb, Journal of Plasma Physics 64, 507 (2000).
• (44) G. M. Webb, M. Pantazopoulou, and G. P. Zank, Journal of Physics A Mathematical General 33, 3137 (2000).
• (45) G. M. Webb, G. P. Zank, E. K. Kaghashvili, and J. A. le Roux, The Astrophysical Journal 651, 211 (2006).
• (46) Y. I. Fedorov et al., Kinematics and Physics of Celestial Bodies 32, 105 (2016).
• (47) A. Shalchi, Astronomy and Astrophys.448, 809 (2006).